Получение старшей части 64-битного целочисленного умножения

Если вы используете gcc и ваша версия поддерживает 128-битные числа (попробуйте использовать __uint128_t), то выполнение 128-битного умножения и извлечение верхних 64-битных чисел, вероятно, будет наиболее эффективным способом получения результата.

Если ваш компилятор не поддерживает 128-битные числа, то ответ Yakk правильный. Однако он может быть слишком кратким для общего употребления. В частности, реальная реализация должна быть осторожна с переполнением 64-битных целых чисел.

Простое и портативное решение, которое он предлагает, состоит в том, чтобы разбить каждое из a и b на 2 32-битных числа, а затем умножить эти 32-битные числа с помощью операции 64-битного умножения. Если мы напишем:

uint64_t a_lo = (uint32_t)a;
uint64_t a_hi = a >> 32;
uint64_t b_lo = (uint32_t)b;
uint64_t b_hi = b >> 32;

то очевидно, что:

a = (a_hi << 32) + a_lo;
b = (b_hi << 32) + b_lo;

а также:

a * b = ((a_hi << 32) + a_lo) * ((b_hi << 32) + b_lo)
      = ((a_hi * b_hi) << 64) +
        ((a_hi * b_lo) << 32) +
        ((b_hi * a_lo) << 32) +
          a_lo * b_lo

при условии, что расчет выполняется с использованием 128-битной (или более высокой) арифметики.

Но эта проблема требует, чтобы мы выполняли все вычисления с использованием 64-битной арифметики, поэтому нам нужно беспокоиться о переполнении.

Поскольку a_hi, a_lo, b_hi и b_lo являются 32-битными числами без знака, их произведение будет соответствовать 64-битному числу без знака без переполнения. Однако промежуточных результатов приведенного выше расчета не будет.

Следующий код будет реализовывать mulhi (a, b), когда математика должна выполняться по модулю 2 ^ 64:

uint64_t    a_lo = (uint32_t)a;
uint64_t    a_hi = a >> 32;
uint64_t    b_lo = (uint32_t)b;
uint64_t    b_hi = b >> 32;

uint64_t    a_x_b_hi =  a_hi * b_hi;
uint64_t    a_x_b_mid = a_hi * b_lo;
uint64_t    b_x_a_mid = b_hi * a_lo;
uint64_t    a_x_b_lo =  a_lo * b_lo;

uint64_t    carry_bit = ((uint64_t)(uint32_t)a_x_b_mid +
                         (uint64_t)(uint32_t)b_x_a_mid +
                         (a_x_b_lo >> 32) ) >> 32;

uint64_t    multhi = a_x_b_hi +
                     (a_x_b_mid >> 32) + (b_x_a_mid >> 32) +
                     carry_bit;

return multhi;
                                              

Как указывает Якк, если вы не против отклонения на +1 в старших 64 битах, вы можете опустить вычисление бита переноса.

TL: DR с GCC для 64-битного ISA: (a * (unsigned __int128)b) >> 64 прекрасно компилируется в одну инструкцию полного умножения или умножения с высокой половиной. Не нужно возиться со встроенным asm.

К сожалению, текущие компиляторы не оптимизируют красивую портативную версию @craigster0, поэтому, если вы хотите использовать преимущества 64-битных процессоров, вы не можете использовать ее, кроме как в качестве запасного варианта. для целей, для которых у вас нет #ifdef. (Я не вижу универсального способа его оптимизации; вам нужен 128-битный тип или встроенный.)

GNU C (gcc, clang или ICC) имеет unsigned __int128 на большинстве 64-битных платформ. (Или в более старых версиях __uint128_t). Однако GCC не реализует этот тип на 32-битных платформах.

Это простой и эффективный способ заставить компилятор выдать 64-битную инструкцию полного умножения и сохранить старшую половину. (GCC знает, что приведение uint64_t к 128-битному целому числу по-прежнему имеет верхнюю половину, равную нулю, поэтому вы не получите 128-битное умножение с использованием трех 64-битных умножений.)

MSVC также имеет __umulh внутреннюю функцию для 64-битного умножения старшей половины, но, опять же, он доступен только на 64-битных платформах (и, в частности, x86-64 и AArch64. В документации также упоминается IPF (IA-64), имеющий _umul128, но у меня нет MSVC для Itanium (в любом случае, вероятно, это не актуально. )

#define HAVE_FAST_mul64 1

#ifdef __SIZEOF_INT128__     // GNU C
 static inline
 uint64_t mulhi64(uint64_t a, uint64_t b) {
     unsigned __int128 prod =  a * (unsigned __int128)b;
     return prod >> 64;
 }

#elif defined(_M_X64) || defined(_M_ARM64)     // MSVC
   // MSVC for x86-64 or AArch64
   // possibly also  || defined(_M_IA64) || defined(_WIN64)
   // but the docs only guarantee x86-64!  Don't use *just* _WIN64; it doesn't include AArch64 Android / Linux

  // https://docs.microsoft.com/en-gb/cpp/intrinsics/umulh
  #include <intrin.h>
  #define mulhi64 __umulh

#elif defined(_M_IA64) // || defined(_M_ARM)       // MSVC again
  // https://docs.microsoft.com/en-gb/cpp/intrinsics/umul128
  // incorrectly say that _umul128 is available for ARM
  // which would be weird because there's no single insn on AArch32
  #include <intrin.h>
  static inline
  uint64_t mulhi64(uint64_t a, uint64_t b) {
     unsigned __int64 HighProduct;
     (void)_umul128(a, b, &HighProduct);
     return HighProduct;
  }

#else

# undef HAVE_FAST_mul64
  uint64_t mulhi64(uint64_t a, uint64_t b);  // non-inline prototype
  // or you might want to define @craigster0's version here so it can inline.
#endif

(или с clang -march=haswell, чтобы включить BMI2: mov rdx, rsi / mulx rax, rcx, rdi, чтобы напрямую поместить верхнюю половину в RAX. gcc тупой и по-прежнему использует дополнительный mov.)

     # x86-64 gcc7.3.  clang and ICC are the same.  (x86-64 System V calling convention)
     # MSVC makes basically the same function, but with different regs for x64 __fastcall
    mov     rax, rsi
    mul     rdi              # RDX:RAX = RAX * RDI
    mov     rax, rdx
    ret

Для AArch64 (с gcc unsigned __int128 или MSVC с __umulh):

С постоянной степенью времени компиляции умножителя 2 мы обычно получаем ожидаемый сдвиг вправо, чтобы захватить несколько старших битов. Но gcc забавно использует shld (см. Ссылку Godbolt).

test_var:
    umulh   x0, x0, x1
    ret

К сожалению, нынешние компиляторы не оптимизируют красивую портативную версию @craigster0. Вы получаете 8x shr r64,32, 4x imul r64,r64 и кучу _21 _ / _ 22_ инструкций для x86-64. то есть он компилируется во множество 32x32 => 64-битных умножений и распаковок результатов. Так что, если вы хотите что-то, что использует преимущества 64-битных процессоров, вам понадобится #ifdefs.

Инструкция полного умножения mul 64 составляет 2 мопа на процессорах Intel, но все же задержка всего 3 цикла, такая же, как imul r64,r64, которая дает только 64-битный результат. Таким образом, __int128 / внутренняя версия в 5-10 раз дешевле по задержкам и пропускной способности (влияние на окружающий код) на современных x86-64, чем портативная версия, судя по быстрому предположению, основанному на http://agner.org/optimize/.

Проверьте это в проводнике компилятора Godbolt по указанной выше ссылке.

Однако gcc полностью оптимизирует эту функцию при умножении на 16: вы получаете один сдвиг вправо, что более эффективно, чем при умножении unsigned __int128.

Это проверенная мною версия, которую я придумал сегодня вечером, обеспечивает полный 128-битный продукт. При осмотре это кажется проще, чем большинство других решений в Интернете (например, в библиотеке Botan и других ответах здесь), потому что оно использует то, как СРЕДНЯЯ ЧАСТЬ не переполняется, как описано в комментариях к коду.

Для контекста я написал это для этого проекта github: https://github.com/catid/fp61

Длинное умножение должно быть нормальным исполнением.

//------------------------------------------------------------------------------
// Portability Macros

// Compiler-specific force inline keyword
#ifdef _MSC_VER
# define FP61_FORCE_INLINE inline __forceinline
#else
# define FP61_FORCE_INLINE inline __attribute__((always_inline))
#endif


//------------------------------------------------------------------------------
// Portable 64x64->128 Multiply
// CAT_MUL128: r{hi,lo} = x * y

// Returns low part of product, and high part is set in r_hi
FP61_FORCE_INLINE uint64_t Emulate64x64to128(
    uint64_t& r_hi,
    const uint64_t x,
    const uint64_t y)
{
    const uint64_t x0 = (uint32_t)x, x1 = x >> 32;
    const uint64_t y0 = (uint32_t)y, y1 = y >> 32;
    const uint64_t p11 = x1 * y1, p01 = x0 * y1;
    const uint64_t p10 = x1 * y0, p00 = x0 * y0;
    /*
        This is implementing schoolbook multiplication:

                x1 x0
        X       y1 y0
        -------------
                   00  LOW PART
        -------------
                00
             10 10     MIDDLE PART
        +       01
        -------------
             01 
        + 11 11        HIGH PART
        -------------
    */

    // 64-bit product + two 32-bit values
    const uint64_t middle = p10 + (p00 >> 32) + (uint32_t)p01;

    /*
        Proof that 64-bit products can accumulate two more 32-bit values
        without overflowing:

        Max 32-bit value is 2^32 - 1.
        PSum = (2^32-1) * (2^32-1) + (2^32-1) + (2^32-1)
             = 2^64 - 2^32 - 2^32 + 1 + 2^32 - 1 + 2^32 - 1
             = 2^64 - 1
        Therefore it cannot overflow regardless of input.
    */

    // 64-bit product + two 32-bit values
    r_hi = p11 + (middle >> 32) + (p01 >> 32);

    // Add LOW PART and lower half of MIDDLE PART
    return (middle << 32) | (uint32_t)p00;
}

#if defined(_MSC_VER) && defined(_WIN64)
// Visual Studio 64-bit

# include <intrin.h>
# pragma intrinsic(_umul128)
# define CAT_MUL128(r_hi, r_lo, x, y) \
    r_lo = _umul128(x, y, &(r_hi));

#elif defined(__SIZEOF_INT128__)
// Compiler supporting 128-bit values (GCC/Clang)

# define CAT_MUL128(r_hi, r_lo, x, y)                   \
    {                                                   \
        unsigned __int128 w = (unsigned __int128)x * y; \
        r_lo = (uint64_t)w;                             \
        r_hi = (uint64_t)(w >> 64);                     \
    }

#else
// Emulate 64x64->128-bit multiply with 64x64->64 operations

# define CAT_MUL128(r_hi, r_lo, x, y) \
    r_lo = Emulate64x64to128(r_hi, x, y);

#endif // End CAT_MUL128

Разделите a*b на (hia+loa)*(hib+lob). Это дает 4 32-битных умножения плюс некоторые сдвиги. Выполняйте их в 64-битном режиме и переносите вручную, и вы получите большую часть.

Обратите внимание, что аппроксимация верхней части может быть сделана с меньшим количеством умножений - с точностью до 2 ^ 33 или около того с 1 умножением и в пределах 1 с 3 умножениями.

Не думаю, что есть портативная альтернатива.

Вот asm для версии ARMv8 или Aarch64:

А вот asm для старых компиляторов DEC:

// High (p1) and low (p0) product
uint64_t p0, p1;
// multiplicand and multiplier
uint64_t a = ..., b = ...;

p0 = a*b; asm ("umulh %0,%1,%2" : "=r"(p1) : "r"(a), "r"(b));

Если у вас есть BMI2 x86 и вы хотите использовать mulxq:

p0 = a*b; p1 = asm("umulh %a0, %a1, %v0", a, b);

И общее умножение x86 с использованием mulq:

asm ("mulxq %3, %0, %1" : "=r"(p0), "=r"(p1) : "d"(a), "r"(b));

Ссылка: blogs.msdn.com/b/oldnewthing/ archive / 2014/12/08 / 10578956.aspx

asm ("mulq %3" : "=a"(p0), "=d"(p1) : "a"(a), "g"(b) : "cc");