Получение старшей части 64-битного целочисленного умножения
В C++ скажем, что:
uint64_t i;
uint64_t j;
затем i * j даст uint64_t что имеет значение как нижнюю часть умножения между i а также jт.е. (i * j) mod 2^64, Теперь, что если бы я хотел высшую часть умножения? Я знаю, что существует инструкция по сборке, которая делает что-то подобное при использовании 32-битных целых чисел, но я совсем не знаком со сборкой, поэтому я надеялся на помощь.
Какой самый эффективный способ сделать что-то вроде:
uint64_t k = mulhi(i, j);
4 ответа
Если вы используете gcc и ваша версия поддерживает 128-битные числа (попробуйте использовать __uint128_t), то выполнение умножения 128 и извлечения старших 64-битных, вероятно, будет наиболее эффективным способом получения результата.
Если ваш компилятор не поддерживает 128-битные числа, то ответ Якка правильный. Тем не менее, он может быть слишком коротким для общего потребления. В частности, фактическая реализация должна быть осторожна с переполнением 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 битах, вы можете опустить расчет бита переноса.
К сожалению, современные компиляторы не оптимизируют красивую портативную версию @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
Для x86-64, AArch64 и PowerPC64 (и других) это компилируется в один mul инструкция и пара movs, чтобы иметь дело с соглашением о вызовах (которое должно оптимизироваться после этого inline). Из проводника компилятора Godbolt (с источником + asm для x86-64, PowerPC64 и AArch64):
# 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
(или с clang -march=haswell включить BMI2: mov rdx, rsi / mulx rax, rcx, rdi положить верхнюю половину в RAX напрямую. GCC тупой и все еще использует дополнительный mov.)
Для AArch64 (с gcc unsigned __int128 или MSVC с __umulh):
test_var:
umulh x0, x0, x1
ret
С постоянной мощностью 2 во время компиляции мы обычно получаем ожидаемое смещение вправо, чтобы получить несколько старших бит. Но gcc забавно использует shld (см. ссылку Godbolt).
К сожалению, современные компиляторы не оптимизируют красивую портативную версию @craigster0. Вы получаете 8x shr r64,32, 4x imul r64,r64 и куча add / mov инструкция для x86-64. т.е. он компилируется с большим количеством 32х32 => 64-битных умножений и распаковывает результаты. Так что если вы хотите что-то, что использует преимущества 64-битных процессоров, вам нужно немного #ifdef s.
Полное умножение 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 и других ответах здесь), потому что оно использует преимущество того, что MIDDLE PART не переполняется, как объяснено в комментариях к коду.
Для контекста я написал это для этого проекта 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:
// 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));
А вот asm для старых компиляторов DEC:
p0 = a*b; p1 = asm("umulh %a0, %a1, %v0", a, b);
Если у вас есть BMI2 x86 и вы хотите использовать mulxq:
asm ("mulxq %3, %0, %1" : "=r"(p0), "=r"(p1) : "d"(a), "r"(b));
И общий x86 умножается с помощью mulq:
asm ("mulq %3" : "=a"(p0), "=d"(p1) : "a"(a), "g"(b) : "cc");