计算128位整数模64位整数的最快方法

您可以使用的Division版本 Russian Peasant Multiplication .

要查找其余部分,请执行(伪代码):

X = B;

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

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

模量保留在A中。

您将需要实现移位、比较和减法,以便对由一对64位数字组成的值进行操作,但这非常简单。

这将最多循环255次(使用128位A)。当然,你需要对零除数做一个预检查。

也许你正在寻找一个完成的程序,但是多精度算法的基本算法可以在Knuth的 Art of Computer Programming ,第2卷。你可以在网上找到描述的分割算法 here . 这些算法处理任意多精度算法,因此比您需要的更通用,但是您应该能够将它们简化为64位或32位数字上的128位算法。为合理的工作量做好准备:(a)理解算法;(b)将其转换为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在32位x86上实现64位模块是一个毛茸茸的循环填充程序集。( 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;

至少可以再进行一次速度优化!在“大除数移位优化”之后,我们可以测试除数高位,如果是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

乘法和加法都可能溢出。我认为人们可以考虑到这一点,并且仍然使用一些修改后的一般概念,但是有一些东西告诉我这将变得非常可怕。

我很好奇64位模块操作是如何在MSVC中实现的,我试图找出一些东西。我真的不知道程序集,我所能得到的只是Express版本,没有vc\crt\src\intel\llrem.asm的源代码,但是我想我在玩了一段调试程序和反汇编输出之后,设法了解了正在发生的事情。我试图搞清楚,如果是正整数,除数=2^32,余数是如何计算的。当然,有一些代码处理负数,但我没有深入研究。

我是这样看的:

如果除数>=2^32,则被除数和除数都会根据需要右移,以使除数适合32位。换言之:如果用n位数字将除数写为二进制,n>32,除数和被除数的n-32最低有效位都将被丢弃。之后,使用硬件支持将64位整数除以32位整数来执行除法。结果可能不正确,但我认为可以证明,结果最多可能会偏离1。除法后,除数(原除数)乘以结果,再减去被除数的积。然后,如果需要的话,通过加或减除数来修正(如果除法的结果被除掉了一个除数)。

利用硬件支持64位32位除法,很容易将128位整数除以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 reduction 效率很高。当然,这假设您多次使用相同的模量,并且将环的元素转换为特殊的表示形式是值得的。

为了对蒙哥马利减少的速度给出一个非常粗略的估计:我有一个旧的基准测试,它在2.4GHz内核2上用64位模和1600纳秒的模幂进行模幂。这个求幂可以完成大约96个模块化乘法(和模块化约简),因此每个模块化乘法需要大约40个周期。

我做了两个版本的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;

@caf接受的答案确实不错,评价也很高,但它包含了多年未见的错误。

为了帮助测试这一点和其他解决方案,我发布了一个测试工具,并让它成为社区wiki。

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的有用构建块。 __umoddi3 (下面链接的源代码)给出了如何执行这类操作的概念,但它只在2n/n=>n除法之上实现2n%2n=>2n,而不是4n%2n=>2n。

提供更广泛的多精度库,例如 https://gmplib.org/manual/Integer-Division.html#Integer-Division .

64位机器上的GNU C 确实提供了 __int128 type 和libgcc函数在目标体系结构上尽可能高效地进行乘法和除法运算。

X86—64 div r/m64 指令执行128B/64B=>64B除法(也将余数作为第二个输出产生),但如果商溢出则会出错。所以如果 A/B > 2^64-1 但是您可以让gcc为您使用它(或者甚至内联libgcc使用的相同代码)。

编译: Godbolt compiler explorer )一两个 div 说明(发生在 libgcc 函数调用)。如果有更快的方法,libgcc可能会使用它。

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

这个 __umodti3 它调用的函数计算一个完整的128b/128b模,但是该函数的实现会检查除数的高半部分为0的特殊情况。 see in the libgcc source . (libgcc根据该代码构建函数的si/di/ti版本,以适合目标体系结构。 udiv_qrnnd 是一个内联asm宏,它对目标体系结构执行无符号2n/n=>n除法。

对于x86-64 (以及其他具有硬件划分指令的架构) 快速路径 (什么时候 high_half(A) < B ;保证 div 不会犯错) 只是两条不带树枝 ,一些错误的CPU需要仔细考虑, 单身 div r64 指令,大约需要50-100个周期 根据 Agner Fog's insn tables . 一些其他的工作可以同时进行 div 但是整数除法单元并不是很流水线 div 解码到许多UOP(不同于FP除法)。

回退路径仍然只使用两个64位 div 案例说明 B 只有64位,但是 A/B 不适合64位,所以 A/B 直接会出错。

注意libgcc的 γ-UMODTI3 只是内联 __udivmoddi4 放入只返回余数的包装器中。

对于相同的重复模 乙

可能值得考虑计算 fixed-point multiplicative inverse 对于 乙 ,如果存在。例如,对于编译时间常量,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 指令执行64b*64b=>128b(rdx:rax)乘法,并可用作构建块来构造128b*128b=>256b乘法以实现相同的算法。因为我们只需要完整256b结果的高一半,这就节省了一些乘数。

现代的英特尔CPU有很高的性能 mul :3c延迟,每时钟吞吐量一个。但是,所需的移位和加法的确切组合随常量的变化而变化,因此在运行时计算乘法逆的一般情况并不如每次使用JIT编译或静态编译版本(甚至在预计算开销之上)时的效率高。

IDK盈亏平衡点在哪里。对于JIT编译,它将高于约200次重用,除非您缓存生成的代码以供常用 乙 价值观。对于“正常”方式,它可能在200次重复使用的范围内,但是IDK为128位/64位除法找到一个模块化的乘法逆运算是多么昂贵。

libdivide 可以这样做,但只能用于32位和64位类型。不过,这可能是一个很好的起点。

一般来说,除法运算速度慢,乘法运算速度快,移位速度快。从我迄今为止看到的答案来看,大多数答案都使用了使用位移位的蛮力方法。还有另一种方式。它是否更快还有待观察(也就是说分析它)。

不用除法,用倒数乘以。因此,要发现a%b,首先计算b的倒数。1/b.这可以通过使用牛顿-拉斐逊收敛方法的几个循环来实现。要做到这一点,将取决于表中一组良好的初始值。

关于牛顿-拉斐逊方法的更多细节,请参考 http://en.wikipedia.org/wiki/Division_(digital)

一旦有了倒数,商q=a*1/b。

余数r=a-q*b。

为了确定这是否会比蛮力更快(因为我们将使用32位寄存器模拟64位和128位数字,所以会有更多的乘法),请对其进行分析。

如果代码中的b是常量,则可以预先计算倒数,并使用最后两个公式进行简单计算。这,我肯定会比移位更快。

希望这有帮助。

如果您最近使用的是x86计算机,那么SSE2+有128位寄存器。我从来没有尝试过为除基本x86以外的任何东西编写程序集,但我怀疑那里有一些指南。

如果128位无符号乘63位无符号足够好,那么它可以在一个循环中完成,最多完成63个循环。

通过将MSN的溢出问题限制为1位,将其视为一个建议的解决方案。我们通过将问题分解为2,模块化乘法,并在最后添加结果来实现。

在下面的示例中,upper对应于最高有效的64位,lower对应于最低有效的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位溢出(信息丢失),从而产生错误的结果。

我还没有想出一个处理溢流的好方法,这让我很恼火。

由于在C中没有预先定义的128位整数类型,a的位必须在数组中表示。虽然B(64位整数)可以存储在 无符号长长整型 变量,为了有效地处理A和B,需要将B的位放入另一个数组中。

之后,b递增为bx2、bx3、bx4,…直到它是最大的B小于A。然后(A-B)可以计算,使用一些基础2的减法知识。

这是您正在寻找的解决方案吗?