Самый быстрый способ вычисления 128-битного целого по модулю 64-разрядного целого числа

Ответ 1

Вы можете использовать версию деления русского крестьянского умножения.

Чтобы найти остаток, выполните (в псевдокоде):

X = B;

while (X <= A/2)
{
    X <<= 1;
}

while (A >= B)
{
    if (A >= X)
        A -= X;
    X >>= 1;
}

Модуль оставлен в А.

Вам нужно будет реализовать сдвиги, сравнения и вычитания, чтобы работать со значениями, состоящими из пары 64-битных чисел, но это довольно тривиально (вероятно, вам следует реализовать сдвиг влево на 1 как X + X).

Это будет цикл не более 255 раз (со 128-битным A). Конечно, вам необходимо выполнить предварительную проверку делителя нуля.

Ответ 2

Возможно, вы ищете готовую программу, но основные алгоритмы для арифметики с несколькими точками можно найти в Knuth Искусство компьютерного программирования, Том 2. Вы можете найти алгоритм деления, описанный онлайн здесь. Алгоритмы имеют дело с произвольной арифметикой с несколькими точками, поэтому они более общие, чем вам нужно, но вы должны иметь возможность упростить их для 128-разрядной арифметики, выполненной с 64- или 32-разрядными цифрами. Будьте готовы к разумному объему работы (а) понимания алгоритма и (б) преобразования его в C или ассемблер.

Вы также можете проверить Hacker Delight, который полон очень умных ассемблеров и других хакерских атак низкого уровня, точность арифметики.

Ответ 3

Учитывая A = AH*2^64 + AL:

A % B == (((AH % B) * (2^64 % B)) + (AL % B)) % B
      == (((AH % B) * ((2^64 - B) % B)) + (AL % B)) % B

Если ваш компилятор поддерживает 64-битные целые числа, то это, вероятно, самый простой способ. MSVC-реализация 64-битного по модулю 32-битного x86 - это сборка, заполненная волосистой петлей (VC\crt\src\intel\llrem.asm для храбрых), поэтому я лично пошел бы с этим.

Ответ 4

Это почти непроверенная, частично измененная модификация Mod128by64 "Русская крестьянская" функция алгоритма. К сожалению, я пользователь Delphi, поэтому эта функция работает в Delphi.:) Но ассемблер почти такой же, как...

function Mod128by64(Dividend: PUInt128; Divisor: PUInt64): UInt64;
//In : eax = @Dividend
//   : edx = @Divisor
//Out: eax:edx as Remainder
asm
//Registers inside rutine
//Divisor = edx:ebp
//Dividend = bh:ebx:edx //We need 64 bits + 1 bit in bh
//Result = esi:edi
//ecx = Loop counter and Dividend index
  push    ebx                     //Store registers to stack
  push    esi
  push    edi
  push    ebp
  mov     ebp, [edx]              //Divisor = edx:ebp
  mov     edx, [edx + 4]
  mov     ecx, ebp                //Div by 0 test
  or      ecx, edx                
  jz      @DivByZero
  xor     edi, edi                //Clear result
  xor     esi, esi
//Start of 64 bit division Loop
  mov     ecx, 15                 //Load byte loop shift counter and Dividend index
@SkipShift8Bits:                  //Small Dividend numbers shift optimisation
  cmp     [eax + ecx], ch         //Zero test
  jnz     @EndSkipShiftDividend
  loop    @SkipShift8Bits         //Skip 8 bit loop
@EndSkipShiftDividend:
  test    edx, $FF000000          //Huge Divisor Numbers Shift Optimisation
  jz      @Shift8Bits             //This Divisor is > $00FFFFFF:FFFFFFFF
  mov     ecx, 8                  //Load byte shift counter
  mov     esi, [eax + 12]         //Do fast 56 bit (7 bytes) shift...
  shr     esi, cl                 //esi = $00XXXXXX
  mov     edi, [eax + 9]          //Load for one byte right shifted 32 bit value
@Shift8Bits:
  mov     bl, [eax + ecx]         //Load 8 bits of Dividend
//Here we can unrole partial loop 8 bit division to increase execution speed...
  mov     ch, 8                   //Set partial byte counter value
@Do65BitsShift:
  shl     bl, 1                   //Shift dividend left for one bit
  rcl     edi, 1
  rcl     esi, 1
  setc    bh                      //Save 65th bit
  sub     edi, ebp                //Compare dividend and  divisor
  sbb     esi, edx                //Subtract the divisor
  sbb     bh, 0                   //Use 65th bit in bh
  jnc     @NoCarryAtCmp           //Test...
  add     edi, ebp                //Return privius dividend state
  adc     esi, edx
@NoCarryAtCmp:
  dec     ch                      //Decrement counter
  jnz     @Do65BitsShift
//End of 8 bit (byte) partial division loop
  dec     cl                      //Decrement byte loop shift counter
  jns     @Shift8Bits             //Last jump at cl = 0!!!
//End of 64 bit division loop
  mov     eax, edi                //Load result to eax:edx
  mov     edx, esi
@RestoreRegisters:
  pop     ebp                     //Restore Registers
  pop     edi
  pop     esi
  pop     ebx
  ret
@DivByZero:
  xor     eax, eax                //Here you can raise Div by 0 exception, now function only return 0.
  xor     edx, edx
  jmp     @RestoreRegisters
end;

Возможно, будет еще одна оптимизация скорости! После "Огромной оптимизации сдвигов чисел Divisor" мы можем проверить высокий бит divisors, если это 0, нам не нужно использовать дополнительный регистр bh как 65-й бит для его хранения. Таким образом, развернутая часть цикла может выглядеть так:

  shl     bl,1                    //Shift dividend left for one bit
  rcl     edi,1
  rcl     esi,1
  sub     edi, ebp                //Compare dividend and  divisor
  sbb     esi, edx                //Subtract the divisor
  jnc     @NoCarryAtCmpX
  add     edi, ebp                //Return privius dividend state
  adc     esi, edx
@NoCarryAtCmpX:

Ответ 5

Решение зависит от того, что именно вы пытаетесь решить.

например. если вы выполняете арифметику в кольце по модулю 64-разрядного целого числа, то используя Уменьшение Montgomerys очень эффективно. Конечно, это предполагает, что вы один и тот же модуль много раз и что он платит, чтобы преобразовать элементы кольца в специальное представление.

Чтобы дать очень приблизительную оценку скорости этого сокращения Montgomerys: у меня есть старый тест, который выполняет модульное возведение в степень с 64-битным модулем и показателем в 1600 нс на 2.4Ghz Core 2. Это возведение в степень составляет около 96 модульные умножения (и модульные сокращения) и, следовательно, требуется около 40 циклов на модульное умножение.

Ответ 6

Я сделал обе версии функции Mod128by64 "русский крестьянин": классическая и оптимизированная по скорости. Оптимизация скорости может делать на моем 3Ghz PC более 1000 000 случайных вычислений в секунду и более чем в три раза быстрее, чем классическая функция. Если мы сравним время выполнения вычисления 128 на 64 и вычисляем 64 на 64 бит по модулю, чем эта функция только на 50% медленнее.

Классический русский крестьянин:

function Mod128by64Clasic(Dividend: PUInt128; Divisor: PUInt64): UInt64;
//In : eax = @Dividend
//   : edx = @Divisor
//Out: eax:edx as Remainder
asm
//Registers inside rutine
//edx:ebp = Divisor
//ecx = Loop counter
//Result = esi:edi
  push    ebx                     //Store registers to stack
  push    esi
  push    edi
  push    ebp
  mov     ebp, [edx]              //Load  divisor to edx:ebp
  mov     edx, [edx + 4]
  mov     ecx, ebp                //Div by 0 test
  or      ecx, edx
  jz      @DivByZero
  push    [eax]                   //Store Divisor to the stack
  push    [eax + 4]
  push    [eax + 8]
  push    [eax + 12]
  xor     edi, edi                //Clear result
  xor     esi, esi
  mov     ecx, 128                //Load shift counter
@Do128BitsShift:
  shl     [esp + 12], 1           //Shift dividend from stack left for one bit
  rcl     [esp + 8], 1
  rcl     [esp + 4], 1
  rcl     [esp], 1
  rcl     edi, 1
  rcl     esi, 1
  setc    bh                      //Save 65th bit
  sub     edi, ebp                //Compare dividend and  divisor
  sbb     esi, edx                //Subtract the divisor
  sbb     bh, 0                   //Use 65th bit in bh
  jnc     @NoCarryAtCmp           //Test...
  add     edi, ebp                //Return privius dividend state
  adc     esi, edx
@NoCarryAtCmp:
  loop    @Do128BitsShift
//End of 128 bit division loop
  mov     eax, edi                //Load result to eax:edx
  mov     edx, esi
@RestoreRegisters:
  lea     esp, esp + 16           //Restore Divisors space on stack
  pop     ebp                     //Restore Registers
  pop     edi                     
  pop     esi
  pop     ebx
  ret
@DivByZero:
  xor     eax, eax                //Here you can raise Div by 0 exception, now function only return 0.
  xor     edx, edx
  jmp     @RestoreRegisters
end;

Оптимизированный по скорости русский крестьянин:

function Mod128by64Oprimized(Dividend: PUInt128; Divisor: PUInt64): UInt64;
//In : eax = @Dividend
//   : edx = @Divisor
//Out: eax:edx as Remainder
asm
//Registers inside rutine
//Divisor = edx:ebp
//Dividend = ebx:edx //We need 64 bits
//Result = esi:edi
//ecx = Loop counter and Dividend index
  push    ebx                     //Store registers to stack
  push    esi
  push    edi
  push    ebp
  mov     ebp, [edx]              //Divisor = edx:ebp
  mov     edx, [edx + 4]
  mov     ecx, ebp                //Div by 0 test
  or      ecx, edx
  jz      @DivByZero
  xor     edi, edi                //Clear result
  xor     esi, esi
//Start of 64 bit division Loop
  mov     ecx, 15                 //Load byte loop shift counter and Dividend index
@SkipShift8Bits:                  //Small Dividend numbers shift optimisation
  cmp     [eax + ecx], ch         //Zero test
  jnz     @EndSkipShiftDividend
  loop    @SkipShift8Bits         //Skip Compute 8 Bits unroled loop ?
@EndSkipShiftDividend:
  test    edx, $FF000000          //Huge Divisor Numbers Shift Optimisation
  jz      @Shift8Bits             //This Divisor is > $00FFFFFF:FFFFFFFF
  mov     ecx, 8                  //Load byte shift counter
  mov     esi, [eax + 12]         //Do fast 56 bit (7 bytes) shift...
  shr     esi, cl                 //esi = $00XXXXXX
  mov     edi, [eax + 9]          //Load for one byte right shifted 32 bit value
@Shift8Bits:
  mov     bl, [eax + ecx]         //Load 8 bit part of Dividend
//Compute 8 Bits unroled loop
  shl     bl, 1                   //Shift dividend left for one bit
  rcl     edi, 1
  rcl     esi, 1
  jc      @DividentAbove0         //dividend hi bit set?
  cmp     esi, edx                //dividend hi part larger?
  jb      @DividentBelow0
  ja      @DividentAbove0
  cmp     edi, ebp                //dividend lo part larger?
  jb      @DividentBelow0
@DividentAbove0:
  sub     edi, ebp                //Return privius dividend state
  sbb     esi, edx
@DividentBelow0:
  shl     bl, 1                   //Shift dividend left for one bit
  rcl     edi, 1
  rcl     esi, 1
  jc      @DividentAbove1         //dividend hi bit set?
  cmp     esi, edx                //dividend hi part larger?
  jb      @DividentBelow1
  ja      @DividentAbove1
  cmp     edi, ebp                //dividend lo part larger?
  jb      @DividentBelow1
@DividentAbove1:
  sub     edi, ebp                //Return privius dividend state
  sbb     esi, edx
@DividentBelow1:
  shl     bl, 1                   //Shift dividend left for one bit
  rcl     edi, 1
  rcl     esi, 1
  jc      @DividentAbove2         //dividend hi bit set?
  cmp     esi, edx                //dividend hi part larger?
  jb      @DividentBelow2
  ja      @DividentAbove2
  cmp     edi, ebp                //dividend lo part larger?
  jb      @DividentBelow2
@DividentAbove2:
  sub     edi, ebp                //Return privius dividend state
  sbb     esi, edx
@DividentBelow2:
  shl     bl, 1                   //Shift dividend left for one bit
  rcl     edi, 1
  rcl     esi, 1
  jc      @DividentAbove3         //dividend hi bit set?
  cmp     esi, edx                //dividend hi part larger?
  jb      @DividentBelow3
  ja      @DividentAbove3
  cmp     edi, ebp                //dividend lo part larger?
  jb      @DividentBelow3
@DividentAbove3:
  sub     edi, ebp                //Return privius dividend state
  sbb     esi, edx
@DividentBelow3:
  shl     bl, 1                   //Shift dividend left for one bit
  rcl     edi, 1
  rcl     esi, 1
  jc      @DividentAbove4         //dividend hi bit set?
  cmp     esi, edx                //dividend hi part larger?
  jb      @DividentBelow4
  ja      @DividentAbove4
  cmp     edi, ebp                //dividend lo part larger?
  jb      @DividentBelow4
@DividentAbove4:
  sub     edi, ebp                //Return privius dividend state
  sbb     esi, edx
@DividentBelow4:
  shl     bl, 1                   //Shift dividend left for one bit
  rcl     edi, 1
  rcl     esi, 1
  jc      @DividentAbove5         //dividend hi bit set?
  cmp     esi, edx                //dividend hi part larger?
  jb      @DividentBelow5
  ja      @DividentAbove5
  cmp     edi, ebp                //dividend lo part larger?
  jb      @DividentBelow5
@DividentAbove5:
  sub     edi, ebp                //Return privius dividend state
  sbb     esi, edx
@DividentBelow5:
  shl     bl, 1                   //Shift dividend left for one bit
  rcl     edi, 1
  rcl     esi, 1
  jc      @DividentAbove6         //dividend hi bit set?
  cmp     esi, edx                //dividend hi part larger?
  jb      @DividentBelow6
  ja      @DividentAbove6
  cmp     edi, ebp                //dividend lo part larger?
  jb      @DividentBelow6
@DividentAbove6:
  sub     edi, ebp                //Return privius dividend state
  sbb     esi, edx
@DividentBelow6:
  shl     bl, 1                   //Shift dividend left for one bit
  rcl     edi, 1
  rcl     esi, 1
  jc      @DividentAbove7         //dividend hi bit set?
  cmp     esi, edx                //dividend hi part larger?
  jb      @DividentBelow7
  ja      @DividentAbove7
  cmp     edi, ebp                //dividend lo part larger?
  jb      @DividentBelow7
@DividentAbove7:
  sub     edi, ebp                //Return privius dividend state
  sbb     esi, edx
@DividentBelow7:
//End of Compute 8 Bits (unroled loop)
  dec     cl                      //Decrement byte loop shift counter
  jns     @Shift8Bits             //Last jump at cl = 0!!!
//End of division loop
  mov     eax, edi                //Load result to eax:edx
  mov     edx, esi
@RestoreRegisters:
  pop     ebp                     //Restore Registers
  pop     edi
  pop     esi
  pop     ebx
  ret
@DivByZero:
  xor     eax, eax                //Here you can raise Div by 0 exception, now function only return 0.
  xor     edx, edx
  jmp     @RestoreRegisters
end;

Ответ 7

Я хотел бы поделиться несколькими мыслями.

Это не так просто, как MSN предлагает, я боюсь.

В выражении:

(((AH % B) * ((2^64 - B) % B)) + (AL % B)) % B

и умножение и добавление могут переполняться. Я думаю, можно было учесть это и по-прежнему использовать общую концепцию с некоторыми изменениями, но что-то говорит мне, что это будет очень страшно.

Мне было любопытно, как в MSVC реализована 64-битная модульная операция, и я попытался найти что-то. Я действительно не знаю сборки, и все, что у меня было, было Express Edition, без источника VC\crt\src\intel\llrem.asm, но я думаю, мне удалось понять, что происходит, после небольшой игры с выходом отладчика и разборки. Я попытался выяснить, как рассчитывается остаток в случае положительных целых чисел и делителя >= 2 ^ 32. Конечно, есть код, который имеет дело с отрицательными номерами, но я не вникал в это.

Вот как я это вижу:

Если divisor >= 2 ^ 32, то и дивиденд, и делитель сдвигаются вправо настолько сильно, насколько это необходимо, чтобы соответствовать делителю на 32 бита. Другими словами: если n цифр для записи делителя вниз в двоичном и n > 32, n-32 младших значащих разрядов как делителя, так и дивиденда отбрасываются. После этого деление выполняется с помощью аппаратной поддержки для деления 64-битных целых чисел на 32-битные. Результат может быть неправильным, но я думаю, что можно доказать, что результат может быть не более 1. После деления делитель (исходный) умножается на результат и произведение вычитается из дивиденда. Затем он корректируется путем добавления или вычитания делителя, если необходимо (если результат разделения был отключен на единицу).

Легко разделить 128-битное целое число на 32-разрядное, используя аппаратную поддержку для 64-битного 32-разрядного деления. В случае, когда делитель < 2 ^ 32, можно вычислить остаток, выполняющий только 4 деления следующим образом:

Предположим, что дивиденд хранится в:

DWORD dividend[4] = ...

остаток войдет в:

DWORD remainder;

1) Divide dividend[3] by divisor. Store the remainder in remainder.
2) Divide QWORD (remainder:dividend[2]) by divisor. Store the remainder in remainder.
3) Divide QWORD (remainder:dividend[1]) by divisor. Store the remainder in remainder.
4) Divide QWORD (remainder:dividend[0]) by divisor. Store the remainder in remainder.

После этих 4 шагов остаток переменной будет содержать то, что вы ищете. (Пожалуйста, не убивайте меня, если я ошибаюсь в энтиансе, я даже не программист)

В случае, если делитель больше, чем 2 ^ 32-1, у меня нет хороших новостей. У меня нет полного доказательства того, что результат после сдвига отключен не более чем на 1 в описанной выше процедуре, которую я считаю MSVC. Я думаю, однако, что это имеет какое-то отношение к тому факту, что часть, которая отбрасывается, по крайней мере в 2 раза меньше, чем в делителе, в 31 раза меньше, дивиденд меньше 2 ^ 64, а делитель больше 2 ^ 32-1, поэтому результат меньше 2 ^ 32.

Если у дивиденда есть 128 бит, трюк с отбрасыванием битов не будет работать. Поэтому в общем случае наилучшим решением является, вероятно, тот, который предлагается GJ или caf. (Ну, было бы, наверное, лучше, даже если бы отбрасывали биты. Разделение, вычитание умножения и коррекция на 128-битное целое число может быть медленнее.)

Я также думал об использовании аппаратного обеспечения с плавающей запятой. Блок с плавающей точкой x87 использует 80-битный формат точности с долей 64 бит в длину. Я думаю, что можно получить точный результат от 64 бит до 64 бит. (Не остаток напрямую, но и остаток с использованием умножения и вычитания, как в "процедуре MSVC" ). Если дивиденд >= 2 ^ 64 и < 2 ^ 128, хранящее его в формате с плавающей запятой, похоже на отбрасывание младших значащих бит в "процедуре MSVC". Может быть, кто-то может доказать, что ошибка в этом случае связана и считает ее полезной. Я понятия не имею, есть ли у него шанс быть быстрее, чем решение GJ, но, возможно, стоит попробовать.

Ответ 8

Принятый ответ @caf был действительно приятным и высоко оцененным, но в нем содержится ошибка, которая не наблюдалась годами.

Чтобы протестировать это и другие решения, я отправляю тестовую ленту и создаю ее вики сообщества.

unsigned cafMod(unsigned A, unsigned B) {
  assert(B);
  unsigned X = B;
  // while (X < A / 2) {  Original code used <
  while (X <= A / 2) {
    X <<= 1;
  }
  while (A >= B) {
    if (A >= X) A -= X;
    X >>= 1;
  }
  return A;
}

void cafMod_test(unsigned num, unsigned den) {
  if (den == 0) return;
  unsigned y0 = num % den;
  unsigned y1 = mod(num, den);
  if (y0 != y1) {
    printf("FAIL num:%x den:%x %x %x\n", num, den, y0, y1);
    fflush(stdout);
    exit(-1);
  }
}

unsigned rand_unsigned() {
  unsigned x = (unsigned) rand();
  return x * 2 ^ (unsigned) rand();
}

void cafMod_tests(void) {
  const unsigned i[] = { 0, 1, 2, 3, 0x7FFFFFFF, 0x80000000, 
      UINT_MAX - 3, UINT_MAX - 2, UINT_MAX - 1, UINT_MAX };
  for (unsigned den = 0; den < sizeof i / sizeof i[0]; den++) {
    if (i[den] == 0) continue;
    for (unsigned num = 0; num < sizeof i / sizeof i[0]; num++) {
      cafMod_test(i[num], i[den]);
    }
  }
  cafMod_test(0x8711dd11, 0x4388ee88);
  cafMod_test(0xf64835a1, 0xf64835a);

  time_t t;
  time(&t);
  srand((unsigned) t);
  printf("%u\n", (unsigned) t);fflush(stdout);
  for (long long n = 10000LL * 1000LL * 1000LL; n > 0; n--) {
    cafMod_test(rand_unsigned(), rand_unsigned());
  }

  puts("Done");
}

int main(void) {
  cafMod_tests();
  return 0;
}

Ответ 9

Я знаю, в вопросе указан 32-битный код, но ответ на 64-битный может быть полезным или интересным для других.

И да, 64b/32b => 32b деление делает полезный строительный блок для 128b% 64b => 64b. libgcc __umoddi3 (источник связан с ниже) дает представление о том, как это сделать, но он реализует только 2N% 2N => 2N поверх деления 2N/N => N, а не 4N% 2N => 2N.

Доступны более широкие библиотеки с высокой точностью, например https://gmplib.org/manual/Integer-Division.html#Integer-Division.

GNU C на 64-битных машинах действительно обеспечивает __int128 тип и функции libgcc для умножения и деления настолько эффективно, насколько это возможно на целевой архитектуре.

div r/m64 x86-64 div r/m64 выполняет div r/m64 128b/64b => 64b (также производит остаток в качестве второго вывода), но не работает, если частное переполнение. Так что вы не можете напрямую использовать его, если A/B > 2^64-1, но вы можете заставить gcc использовать его для себя (или даже встроить тот же код, который использует libgcc).

Это компилирует ( //+godbolt!'s+"binary" mode (link and disassemble) apparently links+a shared libgcc, //because we get a call to __umodti3 that goes+through the PLT. //Most systems+link a static libgcc.a, so __umodti3 is+part of the binary. uint64_t AmodB(unsigned __int128 A, uint64_t B) { return A % B; } //gcc won!'t use 128x128 => high_half(256) multiplication to make this+optimization uint64_t modulo_by_constant(unsigned __int128 A) { return A % 0x12345678ABULL; } uint64_t modulo_by_constant64(uint64_t A) { return A % 0x12345678ABULL; } unsigned __int128 mul128(unsigned __int128 A, unsigned __int128 B) { return A*B%3B+}')),filterAsm:(commentOnly:!t,directives:!t,intel:!t,labels:!t),version:3(link and disassemble) apparently links+a shared libgcc, //because we get a call to __umodti3 that goes+through the PLT. //Most systems+link a static libgcc.a, so __umodti3 is+part of the binary. uint64_t AmodB(unsigned __int128 A, uint64_t B) { return A % B; } uint64_t modulo_by_constant(unsigned __int128 A) { return A % 0x12345678ABULL; } uint64_t modulo_by_constant64(uint64_t A) { return A % 0x12345678ABULL; }')),filterAsm:(commentOnly:!t,directives:!t,intel:!t,labels:!t),version:3 rel="nofollow noreferrer">проводник компилятора Godbolt) в одну или две инструкции div (которые происходят внутри вызова функции libgcc). Если бы был более быстрый путь, libgcc, вероятно, использовал бы это вместо этого.

#include <stdint.h>
uint64_t AmodB(unsigned __int128 A, uint64_t B) {
  return A % B;
}

Функция __umodti3 она вызывает, вычисляет полное 128b/128b по модулю, но реализация этой функции проверяет особый случай, когда старшая половина делителя равна 0, как вы можете видеть в источнике libgcc. (libgcc создает версию функции si/di/ti из этого кода в соответствии с целевой архитектурой. udiv_qrnnd - встроенный макрос asm, который выполняет деление без знака 2N/N => N для целевой архитектуры.

Для x86-64 (и других архитектур с инструкцией аппаратного деления) быстрый путь (когда high_half(A) < B; гарантия того, что div не сработает) - это просто две high_half(A) < B ветки, некоторые из-за пуха из-за По данным таблиц Agner Fog insn, заказываем процессорные процессоры и одну инструкцию div r64, которая занимает около 50-100 циклов на современных процессорах x86. Некоторая другая работа может выполняться параллельно с div, но целочисленная единица деления не очень конвейерна, и div декодирует много мопов (в отличие от деления FP).

Резервный путь все еще использует только две 64-битные инструкции div для случая, когда B является только 64-битным, но A/B не помещается в 64-битные, поэтому A/B напрямую может выйти из строя.

Обратите внимание, что libgcc __umodti3 просто вставляет __udivmoddi4 в оболочку, которая возвращает только остаток.

Для повторения по модулю того же B

Возможно, стоит подумать о вычислении мультипликативного обратного с фиксированной точкой для B, если таковой существует. Например, с константами времени компиляции, gcc выполняет оптимизацию для типов, меньших чем 128b.

uint64_t modulo_by_constant64(uint64_t A) { return A % 0x12345678ABULL; }

    movabs  rdx, -2233785418547900415
    mov     rax, rdi
    mul     rdx
    mov     rax, rdx             # wasted instruction, could have kept using RDX.
    movabs  rdx, 78187493547
    shr     rax, 36            # division result
    imul    rax, rdx           # multiply and subtract to get the modulo
    sub     rdi, rax
    mov     rax, rdi
    ret

Инструкция x86 mul r64 выполняет mul r64 64b * 64b => 128b (rdx: rax) и может использоваться как строительный блок для построения умножения 128b * 128b => 256b для реализации того же алгоритма. Поскольку нам нужна только верхняя половина полного результата 256b, это экономит несколько умножений.

Современные процессоры Intel имеют очень высокую производительность mul: 3c задержки, один за тактовый пропускную способность. Однако точная комбинация требуемых сдвигов и добавлений зависит от константы, поэтому общий случай вычисления мультипликативного обратного значения во время выполнения не столь эффективен каждый раз, когда он используется в качестве JIT-скомпилированной или статически скомпилированной версии (даже на вершине затрат на предварительные вычисления).

ИДК, где будет точка безубыточности. Для JIT-компиляции это будет больше, чем ~ 200 повторного использования, если только вы не кешируете сгенерированный код для часто используемых значений B Для "нормального" способа это может быть в диапазоне 200 повторных использований, но IDK, как дорого было бы найти модульное мультипликативное обратное для 128-битного /64-битного деления.

libdivide может сделать это за вас, но только для 32- и 64-битных типов. Тем не менее, это, вероятно, хорошая отправная точка.

Ответ 10

Как правило, деление медленное и умножение происходит быстрее, а смещение битов еще быстрее. Из того, что я видел в ответах до сих пор, большинство ответов использовало подход грубой силы, используя бит-сдвиги. Существует другой путь. Быстрее ли это еще не видно (AKA-профиль).

Вместо деления умножьте на обратное. Таким образом, чтобы обнаружить A% B, сначала вычислите обратную величину B... 1/B. Это можно сделать с помощью нескольких циклов с использованием метода конвергенции Ньютона-Рафсона. Для этого хорошо зависеть от хорошего набора начальных значений в таблице.

Для получения более подробной информации о методе схождения Ньютона-Рафсона на обратном обратитесь к http://en.wikipedia.org/wiki/Division_(digital)

Как только у вас есть обратное, фактор Q = A * 1/B.

Остаток R = A - Q * B.

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

Если B является постоянным в вашем коде, вы можете предварительно вычислить взаимный и просто вычислить, используя две последние формулы. Это, я уверен, будет быстрее, чем смещение бит.

Надеюсь, что это поможет.

Ответ 11

Если у вас есть последняя машина x86, для SSE2 + существуют 128-битные регистры. Я никогда не пытался писать сборку для чего-то другого, кроме базового x86, но я подозреваю, что там есть несколько руководств.

Ответ 12

Если 128-битная неподписанная по 63-битной без знака достаточно хороша, то это можно сделать в цикле, выполняющем не более 63 циклов.

Рассмотрите это предлагаемое решение проблемы переполнения MSNs, ограничив его 1-бит. Мы делаем это, разбивая задачу на 2, модулярное умножение и добавляя результаты в конце.

В следующем примере верхний соответствует самым значительным 64-битным, младшим и наименее значимым 64-битным, а div - делителем.

unsigned 128_mod(uint64_t upper, uint64_t lower, uint64_t div) {
  uint64_t result = 0;
  uint64_t a = (~0%div)+1;
  upper %= div; // the resulting bit-length determines number of cycles required

  // first we work out modular multiplication of (2^64*upper)%div
  while (upper != 0){
    if(upper&1 == 1){
      result += a;
      if(result >= div){result -= div;}
    }
    a <<= 1;
    if(a >= div){a -= div;}
    upper >>= 1;
  }

  // add up the 2 results and return the modulus
  if(lower>div){lower -= div;}
  return (lower+result)%div;
}

Единственная проблема заключается в том, что если делитель является 64-битным, то мы получаем переполнение 1-разрядной (потеря информации), дающее ошибочный результат.

Мне кажется, что я не понял аккуратного способа обработки переполнения.

Ответ 13

Поскольку в C не существует предопределенного 128-битного целочисленного типа, биты A должны быть представлены в массиве. Хотя B (64-разрядное целое число) может быть сохранено в переменной unsigned long long int, необходимо поместить биты B в другой массив, чтобы эффективно работать с A и B.

После этого B увеличивается в виде Bx2, Bx3, Bx4,... до тех пор, пока он не станет больше B, меньшим A. И тогда (A-B) можно рассчитать, используя некоторые знания о вычитании для базы 2.

Это то решение, которое вы ищете?