54

128 ビットの符号なし整数 A と 64 ビットの符号なし整数 B があります。計算する最も速い方法は何A % Bですか? A を B で割った (64 ビット) 剰余ですか?

これを C またはアセンブリ言語で実行しようとしていますが、32 ビット x86 プラットフォームをターゲットにする必要があります。残念ながら、これは、128 ビット整数のコンパイラ サポートや、単一の命令で必要な操作を実行する x64 アーキテクチャの機能を利用できないことを意味します。

編集:

これまでの回答ありがとうございます。ただし、提案されたアルゴリズムは非常に遅いようです.128ビット×64ビットの除算を実行する最速の方法は、64ビット×32ビットの除算に対するプロセッサのネイティブサポートを活用することではないでしょうか? いくつかの小さな分割に関して、より大きな分割を実行する方法があるかどうかは誰にもわかりませんか?

Re: B はどのくらいの頻度で変わりますか?

主に一般解に興味があります。A と B が毎回異なる可能性がある場合、どのような計算を行いますか?

しかし、2 番目に考えられる状況は、B が A ほど頻繁に変化しないということです。各 B で割るために 200 もの A が存在する可能性があります。この場合、あなたの答えはどのように異なりますか?

4

15 に答える 15

34

ロシアの農民の掛け算の除算バージョンを使用できます。

残りを見つけるには、(疑似コードで) 次を実行します。

X = B;

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

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

係数は A のままです。

64 ビット数値のペアで構成される値を操作するには、シフト、比較、および減算を実装する必要がありますが、それはかなり簡単です (おそらく left-shift-by-1 as を実装する必要がありますX + X)。

これは最大で 255 回ループします (128 ビット A の場合)。もちろん、ゼロ除数の事前チェックを行う必要があります。

于 2010-04-02T12:15:19.590 に答える
13

完成したプログラムを探しているのかもしれませんが、多精度演算の基本的なアルゴリズムは、Knuth のArt of Computer Programming、Volume 2 にあります。オンラインで説明されている除算アルゴリズムは、こちら にあります。アルゴリズムは任意の多倍精度演算を処理するため、必要以上に一般的ですが、64 ビットまたは 32 ビットの数字で実行される 128 ビットの演算用に単純化できるはずです。(a) アルゴリズムを理解し、(b) それを C またはアセンブラーに変換する、妥当な量の作業に備えてください。

Hacker's Delightもチェックしてみてください。これには、非常に巧妙なアセンブラーや、多精度演算を含むその他の低レベルのハッカーが満載です。

于 2010-04-06T06:49:32.513 に答える
12

uint64_t +操作がラップしないほど B が小さい場合:

与えられた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 ビット整数をサポートしている場合、これがおそらく最も簡単な方法です。32 ビット x86 での 64 ビット モジュロの MSVC の実装は、毛むくじゃらのループで満たされたアセンブリ (VC\crt\src\intel\llrem.asm勇敢な人向け) であるため、個人的にはそれを使用します。

于 2010-04-05T03:31:23.117 に答える
8

これは、Mod128by64 'Russian peasant' アルゴリズム関数を部分的に速度変更した、ほとんどテストされていない機能です。残念ながら、私は 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;

少なくとももう 1 つの速度の最適化が可能です! 「Huge Divisor Numbers Shift Optimisation」の後、除数の上位ビットをテストできます。それが 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:
于 2010-04-06T04:37:15.533 に答える
6

質問が32ビットコードを指定していることは知っていますが、64ビットの回答は他の人にとって有用または興味深いかもしれません.

はい、64b/32b => 32b 除算は、128b % 64b => 64b の有用なビルディング ブロックを作成します。libgcc __umoddi3(以下にリンクされているソース) は、そのようなことを行う方法のアイデアを提供しますが、4N % 2N => 2N ではなく、2N / N => N 分割の上に 2N % 2N => 2N のみを実装します。

https://gmplib.org/manual/Integer-Division.html#Integer-Divisionなど、より広い多精度ライブラリが利用可能です。


64 ビット マシン上の GNU C は、ターゲット アーキテクチャ上で可能な限り効率的に乗算および除算するための__int128typeおよび libgcc 関数を提供します。

x86-64 のdiv r/m64命令は 128b/64b => 64b の除算を行いますが (2 番目の出力として剰余も生成します)、商がオーバーフローするとエラーになります。したがって、 if を直接使用することはできませんA/B > 2^64-1が、gcc に使用させることはできます (または、libgcc が使用する同じコードをインライン化することもできます)。


これにより、( Godbolt コンパイラ エクスプローラー) が 1 つまたは 2 つのdiv命令 ( libgcc関数呼び出し内で発生) にコンパイルされます。もっと速い方法があれば、libgcc はおそらくそれを代わりに使用するでしょう。

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

それ__umodti3が呼び出す関数は、完全な 128b/128b モジュロを計算しますが、その関数の実装は、libgcc ソースでわかるように、除数の上位半分が 0 である特別なケースをチェックします。(libgcc は、ターゲット アーキテクチャに応じて、そのコードから関数の si/di/ti バージョンを構築します。ターゲット アーキテクチャに対して udiv_qrnnd符号なし 2N/N => N 除算を行うインライン asm マクロです。

x86-64 (およびハードウェア除算命令を使用するその他のアーキテクチャ) の場合、高速パス( high_half(A) < B; が失敗しないことを保証divする場合)は、2 つの未処理の分岐であり、順序が狂った CPU が噛むための綿毛であり、div r64Agner Fog の insn テーブルによると、最新の x86 CPU で約 50 ~ 100 サイクル1かかる単一の命令。と並行して他の作業を行うことができますdivが、整数除算ユニットはあまりパイプライン化されておらずdiv、多くの uops にデコードされます (FP 除算とは異なります)。

フォールバック パスは、 が 64 ビットのみdivの場合でも 2 つの 64 ビット命令のみを使用しますが、64 ビットに適合しないため、直接エラーが発生します。BA/BA/B

libgccは、残りのみを返すラッパーに__umodti3インライン展開するだけであることに注意してください。__udivmoddi4

脚注 1: 32 ビットdivは、Intel CPU で 2 倍以上高速です。AMD CPU では、実際の入力値が 64 ビット レジスタの小さな値であっても、パフォーマンスはそのサイズにのみ依存します。小さい値が一般的である場合は、64 ビットまたは 128 ビットの除算を行う前に、単純な 32 ビットの除算バージョンへの分岐をベンチマークする価値があるかもしれません。


同じモジュロを繰り返す場合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命令は 64b*64b => 128b (rdx:rax) の乗算を行い、同じアルゴリズムを実装するために 128b * 128b => 256b の乗算を構築するビルディング ブロックとして使用できます。完全な 256b の結果の上位半分だけが必要なので、数倍の節約になります。

最新の Intel CPU は非常に高いパフォーマンスを発揮しmulます。3c のレイテンシ、1 クロックあたり 1 つのスループットです。ただし、必要なシフトと加算の正確な組み合わせは定数によって異なるため、実行時に乗法逆数を計算する一般的なケースは、JIT コンパイルまたは静的にコンパイルされたバージョンとして使用されるたびに効率的ではありません (計算前のオーバーヘッドに加えて)。

損益分岐点となる IDK。JIT コンパイルの場合、一般的に使用される値に対して生成されたコードをキャッシュしない限り、再利用は最大 200 回を超えBます。「通常の」方法では、200回の再利用の範囲にある可能性がありますが、IDKは、128ビット/ 64ビット除算の剰余乗法逆数を見つけるのにどれだけ費用がかかりますか.

libdivideはこれを行うことができますが、32 ビットと 64 ビットの型に対してのみです。それでも、それはおそらく良い出発点です。

于 2016-09-01T08:22:58.747 に答える
4

Mod128by64の「ロシアの農民」除算関数の両方のバージョンを作成しました。クラシックと速度が最適化されています。最適化された速度は、私の3Ghz PCで1秒あたり1000.000以上のランダム計算を実行でき、従来の機能より3倍以上高速です。128x64ビットと64x64ビットを法として計算する実行時間を比較すると、この関数の方が約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;
于 2010-04-10T12:03:04.243 に答える
4

解決策は、正確に何を解決しようとしているかによって異なります。

たとえば、リング モジュロ 64 ビット整数で算術演算を行っている場合、 モンゴメリの簡約を使用すると非常に効率的です。もちろん、これは同じモジュラスを何度も使用し、リングの要素を特別な表現に変換することで利益が得られることを前提としています。


このモンゴメリー削減の速度を非常に大まかに見積もると、2.4Ghz Core 2 で 64 ビットの係数と指数を使用して累乗剰余を 1600 ns で実行する古いベンチマークがあります。この累乗は、約 96 の剰余乗算を行います (およびモジュラ リダクション) であるため、モジュラ乗算ごとに約 40 サイクルが必要です。

于 2010-04-06T21:40:36.530 に答える
4

@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;
}
于 2016-08-31T19:00:43.187 に答える
4

いくつかの考えを共有したいと思います。

残念ながら、MSN が提案するほど単純ではありません。

式では:

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

乗算と加算の両方がオーバーフローする可能性があります。それを考慮して、いくつかの変更を加えた一般的な概念を使用することもできると思いますが、何かが本当に恐ろしくなると言っています.

MSVC で 64 ビット モジュロ演算がどのように実装されているのかに興味があり、何かを見つけようとしました。アセンブリについてはよくわかりません。VC\crt\src\intel\llrem.asm のソースがない Express エディションしかありませんでしたが、少し遊んだ後、何が起こっているのかがわかったと思います。デバッガーと逆アセンブル出力で。正の整数と除数 >=2^32 の場合の剰余の計算方法を理解しようとしました。もちろん、負の数を処理するコードはいくつかありますが、私はそれを掘り下げませんでした。

これが私がそれをどのように見るかです:

除数 >= 2^32 の場合、被除数と除数の両方が、除数が 32 ビットに収まるように必要なだけ右にシフトされます。つまり、除数を 2 進数で書き出すのに n 桁が必要で、n > 32 の場合、除数と被除数の両方の最下位 n-32 桁が破棄されます。その後、64 ビット整数を 32 ビット整数で除算するためのハードウェア サポートを使用して除算が実行されます。結果は間違っているかもしれませんが、結果が最大で 1 ずれている可能性があることは証明できると思います。次に、必要に応じて除数を加算または減算して修正します (除算の結果が 1 ずれている場合)。

64 ビット x 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 より大きい場合、良いニュースはありません。MSVC が使用していると思われる前述の手順では、シフト後の結果が 1 を超えないという完全な証拠はありません。ただし、破棄される部分が少なくとも除数の 2 ^ 31 倍少なく、被除数が 2 ^ 64 未満で、除数が 2 ^ 32-1 より大きいという事実と関係があると思います。であるため、結果は 2^32 未満です。

被除数が 128 ビットの場合、ビットを破棄するトリックは機能しません。したがって、一般的に、最良の解決策はおそらく GJ または caf によって提案されたものです。(ビットの破棄が機能したとしても、おそらくそれが最善でしょう。128 ビット整数の除算、乗算、減算、および修正は遅くなる可能性があります。)

浮動小数点ハードウェアの使用も考えていました。x87 浮動小数点ユニットは、小数部が 64 ビット長の 80 ビット精度形式を使用します。64ビット×64ビットの除算の正確な結果を得ることができると思います。(直接の剰余ではなく、「MSVC 手順」のように乗算と減算を使用した剰余も)。被除数 >=2^64 かつ < 2^128 の場合、それを浮動小数点形式で格納することは、「MSVC 手順」で最下位ビットを破棄することに似ています。おそらく誰かがその場合のエラーがバインドされていることを証明し、それが役立つことを発見できるでしょう。GJのソリューションよりも高速になる可能性があるかどうかはわかりませんが、試してみる価値があるかもしれません.

于 2010-04-10T23:00:11.487 に答える
2

符号なし 128 ビット×符号なし 63 ビットで十分な場合は、最大 63 サイクルのループで実行できます。

これを 1 ビットに制限することにより、MSN のオーバーフロー問題の提案された解決策と考えてください。問題を 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 ビットのオーバーフロー (情報の損失) が発生して、誤った結果が生じることです。

オーバーフローを処理するきちんとした方法を見つけられなかったことが私を悩ませています。

于 2016-11-03T10:49:50.050 に答える
2

原則として、除算は遅く、乗算は高速であり、ビット シフトはさらに高速です。これまでの回答から、ほとんどの回答はビットシフトを使用したブルート フォース アプローチを使用していました。別の方法があります。それがより速いかどうかはまだわかりません(別名、それをプロファイルします)。

割る代わりに、逆数を掛けます。したがって、A % B を見つけるには、まず B の逆数を計算します ... 1/B. これは、収束のニュートン ラフソン法を使用していくつかのループで行うことができます。これをうまく行うには、テーブル内の適切な初期値のセットに依存します。

逆数に収束するニュートン ラフソン法の詳細については、http://en.wikipedia.org/wiki/Division_(digital)を参照してください。

逆数を求めると、商は Q = A * 1/B になります。

剰余 R = A - Q*B。

これがブルート フォースよりも高速かどうかを判断するには (32 ビット レジスタを使用して 64 ビットと 128 ビットの数値をシミュレートするため、さらに多くの乗算が行われるため、プロファイリングします。

コード内で B が定数の場合、逆数を事前に計算し、最後の 2 つの式を使用して簡単に計算できます。これは、ビットシフトよりも高速になると確信しています。

お役に立てれば。

于 2010-04-08T11:03:55.900 に答える
0

最近の x86 マシンを使用している場合は、SSE2+ 用の 128 ビット レジスタがあります。基本的な x86 以外のアセンブリを作成しようとしたことはありませんが、いくつかのガイドがあると思います。

于 2010-04-02T20:08:47.633 に答える
-2

C には定義済みの 128 ビット整数型がないため、A のビットは配列で表現する必要があります。B (64 ビット整数) はunsigned long long int変数に格納できますが、A と B を効率的に処理するには、B のビットを別の配列に配置する必要があります。

その後、B は Bx2、Bx3、Bx4 のようにインクリメントされ、A よりも小さい最大の B になります。そして、基数 2 の減算の知識を使用して、(AB) を計算できます。

これはあなたが探している種類のソリューションですか?

于 2010-04-02T11:04:47.077 に答える