Эффективный способ сделать основные 128-битные целочисленные вычисления в C++?
Несколько лет назад мне нужен был способ сделать некоторую базовую 128-битную целочисленную математику с Cuda: 128-битное целое число на cuda?, Теперь у меня та же проблема, но на этот раз мне нужно запустить базовую 128-битную арифметику (суммы, сдвиги и умножения) в 32-битной встроенной системе (Intel Edison), которая не поддерживает 128-битные типы данных. Однако напрямую поддерживаются 64-битные целые числа (unsigned long long int).
Я наивно пытался использовать код asm, который мне ответили в прошлый раз на процессоре, но я получил кучу ошибок. Я не очень разбираюсь в asm, так что: как наиболее эффективно, имея 64-битные целые числа, реализовывать сложения, умножения и сдвиг битов в 128 битах?
2 ответа
Обновление: поскольку ОП еще не принял ответ <подсказка> <подсказка>, я добавил немного больше кода.
Использование библиотек, обсужденных выше, вероятно, хорошая идея. Хотя сегодня вам может понадобиться всего несколько функций, в конечном итоге вы можете обнаружить, что вам нужна еще одна. Затем еще один после этого. До тех пор, пока вы не закончите писать, отлаживать и поддерживать свою собственную 128-битную математическую библиотеку. Что является пустой тратой вашего времени и усилий.
Это сказал. Если вы полны решимости бросить свой собственный:
1) Вопрос cuda, который вы задавали ранее, уже имеет код c для умножения. Была ли какая-то проблема с этим?
2) Сдвиг, вероятно, не принесет пользы от использования asm, поэтому и здесь решение для ac имеет смысл. Хотя, если производительность действительно является проблемой здесь, я бы посмотрел, поддерживает ли Edison SHLD/SHRD, что может сделать это немного быстрее. Иначе, м. Может быть, такой подход?
my_uint128_t lshift_uint128 (const my_uint128_t a, int b)
{
my_uint128_t res;
if (b < 32) {
res.x = a.x << b;
res.y = (a.y << b) | (a.x >> (32 - b));
res.z = (a.z << b) | (a.y >> (32 - b));
res.w = (a.w << b) | (a.z >> (32 - b));
} elseif (b < 64) {
...
}
return res;
}
Обновление: так как кажется, что Edison может поддерживать SHLD/SHRD, вот альтернатива, которая может быть более производительной, чем код 'c' выше. Как и весь код, предназначенный для ускорения, вы должны протестировать его.
inline
unsigned int __shld(unsigned int into, unsigned int from, unsigned int c)
{
unsigned int res;
if (__builtin_constant_p(into) &&
__builtin_constant_p(from) &&
__builtin_constant_p(c))
{
res = (into << c) | (from >> (32 - c));
}
else
{
asm("shld %b3, %2, %0"
: "=rm" (res)
: "0" (into), "r" (from), "ic" (c)
: "cc");
}
return res;
}
inline
unsigned int __shrd(unsigned int into, unsigned int from, unsigned int c)
{
unsigned int res;
if (__builtin_constant_p(into) &&
__builtin_constant_p(from) &&
__builtin_constant_p(c))
{
res = (into >> c) | (from << (32 - c));
}
else
{
asm("shrd %b3, %2, %0"
: "=rm" (res)
: "0" (into), "r" (from), "ic" (c)
: "cc");
}
return res;
}
my_uint128_t lshift_uint128 (const my_uint128_t a, unsigned int b)
{
my_uint128_t res;
if (b < 32) {
res.x = a.x << b;
res.y = __shld(a.y, a.x, b);
res.z = __shld(a.z, a.y, b);
res.w = __shld(a.w, a.z, b);
} else if (b < 64) {
res.x = 0;
res.y = a.x << (b - 32);
res.z = __shld(a.y, a.x, b - 32);
res.w = __shld(a.z, a.y, b - 32);
} else if (b < 96) {
res.x = 0;
res.y = 0;
res.z = a.x << (b - 64);
res.w = __shld(a.y, a.x, b - 64);
} else if (b < 128) {
res.x = 0;
res.y = 0;
res.z = 0;
res.w = a.x << (b - 96);
} else {
memset(&res, 0, sizeof(res));
}
return res;
}
my_uint128_t rshift_uint128 (const my_uint128_t a, unsigned int b)
{
my_uint128_t res;
if (b < 32) {
res.x = __shrd(a.x, a.y, b);
res.y = __shrd(a.y, a.z, b);
res.z = __shrd(a.z, a.w, b);
res.w = a.w >> b;
} else if (b < 64) {
res.x = __shrd(a.y, a.z, b - 32);
res.y = __shrd(a.z, a.w, b - 32);
res.z = a.w >> (b - 32);
res.w = 0;
} else if (b < 96) {
res.x = __shrd(a.z, a.w, b - 64);
res.y = a.w >> (b - 64);
res.z = 0;
res.w = 0;
} else if (b < 128) {
res.x = a.w >> (b - 96);
res.y = 0;
res.z = 0;
res.w = 0;
} else {
memset(&res, 0, sizeof(res));
}
return res;
}
3) дополнение может извлечь выгоду из асм. Вы можете попробовать это:
struct my_uint128_t
{
unsigned int x;
unsigned int y;
unsigned int z;
unsigned int w;
};
my_uint128_t add_uint128 (const my_uint128_t a, const my_uint128_t b)
{
my_uint128_t res;
asm ("addl %5, %[resx]\n\t"
"adcl %7, %[resy]\n\t"
"adcl %9, %[resz]\n\t"
"adcl %11, %[resw]\n\t"
: [resx] "=&r" (res.x), [resy] "=&r" (res.y),
[resz] "=&r" (res.z), [resw] "=&r" (res.w)
: "%0"(a.x), "irm"(b.x),
"%1"(a.y), "irm"(b.y),
"%2"(a.z), "irm"(b.z),
"%3"(a.w), "irm"(b.w)
: "cc");
return res;
}
Я просто накатал это, так что пользуйтесь на свой страх и риск. У меня нет Эдисона, но это работает с x86.
Обновление: если вы просто делаете накопления (подумайте to += from вместо кода выше, который c = a + b), этот код может послужить вам лучше:
inline
void addto_uint128 (my_uint128_t *to, const my_uint128_t from)
{
asm ("addl %[fromx], %[tox]\n\t"
"adcl %[fromy], %[toy]\n\t"
"adcl %[fromz], %[toz]\n\t"
"adcl %[fromw], %[tow]\n\t"
: [tox] "+&r"(to->x), [toy] "+&r"(to->y),
[toz] "+&r"(to->z), [tow] "+&r"(to->w)
: [fromx] "irm"(from.x), [fromy] "irm"(from.y),
[fromz] "irm"(from.z), [fromw] "irm"(from.w)
: "cc");
}
Если использование внешней библиотеки является опцией, взгляните на этот вопрос. Вы можете использовать TTMath, который является очень простым заголовком для математики с большой точностью. На 32-битных архитектурах ttmath:UInt<4> создаст 128-бит int тип с четырьмя 32-битными конечностями.
Если вы должны написать это самостоятельно, то на SO уже есть много решений, и я их кратко изложу здесь.
Для сложения и вычитания это очень легко и просто, просто добавьте / вычтите слова (которые большие int-библиотеки часто называют конечностями) от младшего значащего до более значимого, конечно, с переносом.
typedef struct INT128 {
uint64_t H, L;
} my_uint128_t;
inline my_uint128_t add(my_uint128_t a, my_uint128_t b)
{
my_uint128_t c;
c.L = a.L + b.L;
c.H = a.H + b.H + (c.L < a.L); // c = a + b
return c;
}
Вывод сборки можно проверить с помощью Compiler Explorer
Компиляторы уже могут сгенерировать эффективный код для операций с двумя словами, но многие недостаточно умны, чтобы использовать "добавить с переносом" при компиляции операций с несколькими словами из языков высокого уровня, как вы можете видеть в вопросе эффективного 128-битного сложения с использованием нести флаг. Поэтому, используя 2 long long Подобные действия сделают компилятор не только более читабельным, но и более легким для создания немного более эффективного кода.
Если это по-прежнему не соответствует вашим требованиям к производительности, вы должны использовать встроенную функцию или написать ее в сборке. Чтобы добавить 128-битное значение, хранящееся в bignum для 128-битного значения в {eax, ebx, ecx, edx} вы можете использовать следующий код
add edx, [bignum]
adc ecx, [bignum+4]
adc ebx, [bignum+8]
adc eax, [bignum+12]
Эквивалентная внутренняя будет такой же для Clang
unsigned *x, *y, *z, carryin=0, carryout;
z[0] = __builtin_addc(x[0], y[0], carryin, &carryout);
carryin = carryout;
z[1] = __builtin_addc(x[1], y[1], carryin, &carryout);
carryin = carryout;
z[2] = __builtin_addc(x[2], y[2], carryin, &carryout);
carryin = carryout;
z[3] = __builtin_addc(x[3], y[3], carryin, &carryout);
Вам нужно изменить внутреннее на то, которое поддерживается вашим компилятором, например __builtin_uadd_overflow в gcc или _addcarry_u32 для MSVC и ICC
Для получения дополнительной информации прочитайте эти
- Работа с большими числами с использованием инструкций x86
- Как я могу сложить и вычесть 128-битные целые числа в C или C++, если мой компилятор не поддерживает их?
- Создание хорошего дополнения с кодом переноса из Clang
- добавление нескольких слов с использованием флага переноса
Для сдвигов битов решение C можно найти в вопросе " Операция побитового сдвига на 128-битном числе". Это простой сдвиг влево, но вы можете развернуть рекурсивный вызов для большей производительности
void shiftl128 (
unsigned int& a,
unsigned int& b,
unsigned int& c,
unsigned int& d,
size_t k)
{
assert (k <= 128);
if (k >= 32) // shifting a 32-bit integer by more than 31 bits is "undefined"
{
a=b;
b=c;
c=d;
d=0;
shiftl128(a,b,c,d,k-32);
}
else
{
a = (a << k) | (b >> (32-k));
b = (b << k) | (c >> (32-k));
c = (c << k) | (d >> (32-k));
d = (d << k);
}
}
Сборка для менее чем 32-битных сдвигов может быть найдена в вопросе 128-битных сдвигов с использованием языка ассемблера?
shld edx, ecx, cl
shld ecx, ebx, cl
shld ebx, eax, cl
shl eax, cl
Сдвиг вправо может быть реализован аналогично или просто скопирован из приведенного выше связанного вопроса
Умножение и деление намного сложнее, и вы можете сослаться на решение в вопросе " Эффективное умножение / деление двух 128-битных целых чисел на x86 (без 64-битной)":
class int128_t
{
uint32_t dw3, dw2, dw1, dw0;
// Various constrctors, operators, etc...
int128_t& operator*=(const int128_t& rhs) __attribute__((always_inline))
{
int128_t Urhs(rhs);
uint32_t lhs_xor_mask = (int32_t(dw3) >> 31);
uint32_t rhs_xor_mask = (int32_t(Urhs.dw3) >> 31);
uint32_t result_xor_mask = (lhs_xor_mask ^ rhs_xor_mask);
dw0 ^= lhs_xor_mask;
dw1 ^= lhs_xor_mask;
dw2 ^= lhs_xor_mask;
dw3 ^= lhs_xor_mask;
Urhs.dw0 ^= rhs_xor_mask;
Urhs.dw1 ^= rhs_xor_mask;
Urhs.dw2 ^= rhs_xor_mask;
Urhs.dw3 ^= rhs_xor_mask;
*this += (lhs_xor_mask & 1);
Urhs += (rhs_xor_mask & 1);
struct mul128_t
{
int128_t dqw1, dqw0;
mul128_t(const int128_t& dqw1, const int128_t& dqw0): dqw1(dqw1), dqw0(dqw0){}
};
mul128_t data(Urhs,*this);
asm volatile(
"push %%ebp \n\
movl %%eax, %%ebp \n\
movl $0x00, %%ebx \n\
movl $0x00, %%ecx \n\
movl $0x00, %%esi \n\
movl $0x00, %%edi \n\
movl 28(%%ebp), %%eax #Calc: (dw0*dw0) \n\
mull 12(%%ebp) \n\
addl %%eax, %%ebx \n\
adcl %%edx, %%ecx \n\
adcl $0x00, %%esi \n\
adcl $0x00, %%edi \n\
movl 24(%%ebp), %%eax #Calc: (dw1*dw0) \n\
mull 12(%%ebp) \n\
addl %%eax, %%ecx \n\
adcl %%edx, %%esi \n\
adcl $0x00, %%edi \n\
movl 20(%%ebp), %%eax #Calc: (dw2*dw0) \n\
mull 12(%%ebp) \n\
addl %%eax, %%esi \n\
adcl %%edx, %%edi \n\
movl 16(%%ebp), %%eax #Calc: (dw3*dw0) \n\
mull 12(%%ebp) \n\
addl %%eax, %%edi \n\
movl 28(%%ebp), %%eax #Calc: (dw0*dw1) \n\
mull 8(%%ebp) \n\
addl %%eax, %%ecx \n\
adcl %%edx, %%esi \n\
adcl $0x00, %%edi \n\
movl 24(%%ebp), %%eax #Calc: (dw1*dw1) \n\
mull 8(%%ebp) \n\
addl %%eax, %%esi \n\
adcl %%edx, %%edi \n\
movl 20(%%ebp), %%eax #Calc: (dw2*dw1) \n\
mull 8(%%ebp) \n\
addl %%eax, %%edi \n\
movl 28(%%ebp), %%eax #Calc: (dw0*dw2) \n\
mull 4(%%ebp) \n\
addl %%eax, %%esi \n\
adcl %%edx, %%edi \n\
movl 24(%%ebp), %%eax #Calc: (dw1*dw2) \n\
mull 4(%%ebp) \n\
addl %%eax, %%edi \n\
movl 28(%%ebp), %%eax #Calc: (dw0*dw3) \n\
mull (%%ebp) \n\
addl %%eax, %%edi \n\
pop %%ebp \n"
:"=b"(this->dw0),"=c"(this->dw1),"=S"(this->dw2),"=D"(this->dw3)
:"a"(&data):"%ebp");
dw0 ^= result_xor_mask;
dw1 ^= result_xor_mask;
dw2 ^= result_xor_mask;
dw3 ^= result_xor_mask;
return (*this += (result_xor_mask & 1));
}
};
Вы также можете найти много похожих вопросов с помощью 128-битного тега.