Сборка — эффективный способ сделать основные 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, поэтому решение c также имеет смысл для меня и здесь. Хотя, если производительность действительно проблема здесь, я бы посмотрел, поддерживает ли 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, вот альтернатива, которая может быть быть более производительным, чем код «с» выше. Как и весь код, предназначенный для ускорения, вы должны проверить его.

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");
}
4

Другие решения

При использовании внешняя библиотека вариант тогда взглянуть на этот вопрос. Ты можешь использовать 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;
}

Выход сборки можно проверить с помощью Проводник компилятора

Компиляторы уже могут генерировать эффективный код для операций с двумя словами, но многие недостаточно умны, чтобы использовать «добавить с переносом» при компиляции операций из нескольких слов из языков высокого уровня, как вы можете видеть в вопросе. эффективное 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

Для получения дополнительной информации прочитайте эти


За сдвиги бит Вы можете найти решение 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));
}
};

Вы также можете найти много связанных вопросов с тег

2

По вопросам рекламы [email protected]