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

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

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

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

Модуль остается в A.

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

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

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

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

Учитывая, что 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 – это assembly, заполненная волосистой петлей ( VC\crt\src\intel\llrem.asm для храбрых), поэтому я лично поступил бы с этим.

Это почти непроверенная, частично измененная модификация 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: 

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

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

В выражении:

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

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

Мне было любопытно, как в MSVC реализована 64-битная модульная операция, и я попытался найти что-то. Я действительно не знаю сборки, и все, что у меня было, было Express Edition, без источника VC \ crt \ src \ intel \ llrem.asm, но я думаю, мне удалось понять, что происходит, после небольшой игры с выходом отладчика и parsingки. Я попытался выяснить, как рассчитывается остаток в случае положительных целых чисел и делителя> = 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, но, возможно, стоит попробовать.

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

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

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

Я сделал обе версии функции Mod128by64 «русский крестьянин»: classический и оптимизированный по скорости. Оптимизация скорости может делать на моем 3Ghz PC более 1000 000 случайных вычислений в секунду и более чем в три раза быстрее, чем classическая функция. Если мы сравним время выполнения вычисления 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; 

Принятый ответ от @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; } 

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

И да, 64b / 32b => 32b деление делает полезный строительный блок для 128b% 64b => 64b. libgcc’s __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 выполняет 128b / 64b => 64b-деление (также создавая остаток в качестве второго выхода), но он ошибочно, если фактор переполняется. Поэтому вы не можете напрямую использовать его, если A/B > 2^64-1 , но вы можете получить gcc, чтобы использовать его для вас (или даже встроить тот же код, который использует libgcc).

Это компилирует ( Godbolt compiler explorer ) одну или две инструкции div (которые происходят внутри вызова функции libgcc ). Если бы был более быстрый способ, libgcc, вероятно, использовал бы это вместо этого.

 #include  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 ветви , некоторый пух для вне- в соответствии с div r64 , и для одного div r64 , который занимает около 50-100 циклов на современных процессорах x86. Некоторая другая работа может происходить параллельно с div , но целочисленная единица деления не очень конвейерна, а div декодируется на множество uops (в отличие от разделения FP).

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

Обратите внимание, что libgcc's __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 

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

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

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

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

Как правило, деление происходит медленно, а умножение происходит быстрее, а смещение бит еще быстрее. Из того, что я видел в ответах до сих пор, большинство ответов использовало подход грубой силы, используя бит-сдвиги. Существует другой путь. Быстрее ли это видно (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 является постоянным в вашем коде, вы можете предварительно вычислить взаимный и просто вычислить, используя две последние формулы. This, I am sure will be faster than bit-shifting.

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

If you have a recent x86 machine, there are 128-bit registers for SSE2+. I’ve never tried to write assembly for anything other than basic x86, but I suspect there are some guides out there.

If 128-bit unsigned by 63-bit unsigned is good enough, then it can be done in a loop doing at most 63 cycles.

Consider this a proposed solution MSNs’ overflow problem by limiting it to 1-bit. We do so by splitting the problem in 2, modular multiplication and adding the results at the end.

In the following example upper corresponds to the most significant 64-bits, lower to the least significant 64-bits and div is the divisor.

 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; } 

The only problem is that, if the divisor is 64-bits then we get overflows of 1-bit (loss of information) giving a faulty result.

It bugs me that I haven’t figured out a neat way to handle the overflows.

Since there is no predefined 128-bit integer type in C, bits of A have to be represented in an array. Although B (64-bit integer) can be stored in an unsigned long long int variable, it is needed to put bits of B into another array in order to work on A and B efficiently.

After that, B is incremented as Bx2, Bx3, Bx4, … until it is the greatest B less than A. And then (AB) can be calculated, using some subtraction knowledge for base 2.

Is this the kind of solution that you are looking for?