Zenn
🚄

多倍長整数の除算: Burnikel-Ziegler 法

2024/02/20に公開
C#
Go
Java
JavaScript
多倍長整数
tech

多倍長整数の除算: Burnikel-Ziegler 法

この記事では基数を \beta と表します。2 進数ならば \beta=2, 10 進数ならば \beta=10 です。

r 桁の整数を s 桁の整数で割るときの計算量は、ナイーブな実装の場合 O(s(r-s)) です。

r=n, s=\dfrac{n}{2} の場合は s(r-s)= \dfrac{n^2}{4} なので O(n^2) となってしまいます。

乗算の計算量を K(n) としたときに O(K(n) + n \log n) で除算できるアルゴリズムはいくつかあります。

競技プログラミングでの実装はニュートン法で実装している人が多い印象ですが、精度の管理などややこしいので実用には実装が簡潔でわかりやすい Burnikel-Ziegler 法 がよく用いられているようです。

参考: 競プロでの提出例

実用されている例としては Go, Java(OpenJDK), JavaScript(v8) のランタイム実装で使われています。

また、.NET の BigInteger は愚直計算による実装で遅かったので Burnikel-Ziegler 法での実装での PR を出しました。

この記事の説明は Burnikel と Ziegler の論文 Fast Recursive Division の説明に基づきますが、変数名は変更しています。

導入

r 桁の整数 A を s 桁の整数 B で割る除算を考えます。商を Q、あまりを R とします。

ここで r \ge s とします。r \lt s ならば自明に Q=0, R=A です。

また、下記のような表記を導入します。

\begin{aligned} A &= [A_2, A_1, A_0] = A_2 \beta^{2n/3} + A_1 \beta^{n/3} + A_0 \\ A &= [A_1, A_0] = A_1 \beta^{n/2} + A_0 \end{aligned}
  • [A_3, A_2, A_1] と書いたときは A の桁数 n が 3 の倍数であり、下位 n/3 桁を A_0、その次の n/3 桁を A_1、上位 n/3 桁を A_2 とする。
  • [A_2, A_1] と書いたときは A の桁数 n が 2 の倍数であり、下位 n/2 桁を A_0、上位 n/2 桁を A_1 とする。

という具合です。

桁が小さいとき・桁の差が小さいときの除算(愚直な計算)

B の桁が小さいとき、または A, B の桁の差が小さいときは上の方の桁から順に商を確定させていく愚直な計算でやります。

ひとまず B が DIV\_LIMIT 桁未満ならば愚直計算するということにします。 DIV\_LIMIT の決め方は ベンチマークに基づくのが良いでしょう。自分の計測だと 32 になりました。

論文の RecursiveDivision 関数では DIV\_LIMIT 桁以下という条件になっていますが、後述の補正を考えると「未満」とするのが妥当です。

2n 桁 / n 桁の除算

r \le 2n, s=n となるときの除算です。といっても n は偶数でなければならないので l=4p, m=2p 桁と考える方が良いかもしれません。また、A \lt B \beta^n 、言い換えると 商 Q が n 桁以下である必要があります。

n が奇数だったり DIV\_LIMIT 以下だったりするときは愚直に計算します。

そうでなければ、A=[A_3,A_2,A_1,A_0], B=[B_1,B_0], Q=[Q_1,Q_0] として次のように考えます。

  1. 3n 桁 / 2n 桁の除算により [A_3,A_2,A_1] = Q_1 B + S となる Q_1, S を求める。
  2. 3n 桁 / 2n 桁の除算により [S_1,S_0,A_0] = Q_0 B + R となる Q_0, R を求める。

このようにすると 2 回の 3n 桁 / 2n 桁の除算に落とし込むことで計算できます。

2n 桁 / n 桁の除算の図解

A/B を考える
|       A_3       |       A_2       |       A_1       |       A_0       | / |       B_1       |       B_0       |

Aは以下の和で表される
|                        Q_1 B                        |                 |
|                 |       S_1       |       S_0       |       A_0       |
=
|                        Q_1 B                        |                 |
|                 |                        Q_0 B                        |
|                 |                 |       R_1       |       R_0       |

3n 桁 / 2n 桁の除算

r \le 3n, s=2n となるときの除算です。

以下、A=[A_2,A_1,A_0], B=[B_1,B_0] として考えます。

A \lt B \beta^n すなわち 商 Q が n 桁以下である必要があります。分割形式で表すと [A_2,A_1] \lt [B_1,B_0] です。

  1. 商の概算 \widehat{Q} を求める

Q \le \widehat{Q} \le Q+2 であるような \widehat{Q} を求めます。

  • A_2 \lt B_1 ならば \widehat{Q} = [A_2, A_1]/B_1 (2n 桁 / n 桁の除算)
  • A_2 \ge B_1 ならば \widehat{Q} = \beta^n - 1
    • A \lt B \beta^n なのだから A_2 \gt B_1 にはならない気がする

また、S = [A_2, A_1] - \widehat{Q} B_1 となるような S を求めます。言い換えると、S は [A_2, A_1]/B_1 のあまりです。

Q \le \widehat{Q} \le Q+2 をみたすことの証明は論文に記載されています。

  1. あまりの概算 \widehat{R} を求める

まず

D = \widehat{Q} B_2

を導出します。A_2 \ge B_1 のときは D = B_2 \beta^n - B_2 とすると良いです。

すると

\widehat{R} = S \beta^n + A_0 - D

として、あまりの概算が求められます。

  1. 実際の商とあまりの導出

Q \le \widehat{Q} \le Q+2 なので \widehat{Q} は実際の商より少し大きい可能性があります。

Q \lt \widehat{Q} のときは \widehat{R} が負の値となっているので

// 擬似コード
// 実際には配列で保持されている
R_hat = S * pow(β, n) + A_0 - D;
while(R_hat < 0)
{
    R_hat += B;
    Q_hat++;
}

というようなロジックを記述して、\widehat{Q}, \widehat{R} を実際の Q, R にします。

  1. 補足

実際の値は配列で保持しているため、値は符号なしの非負整数として扱いたいです。そのため、「\widehat{R} が負」という条件でループするのは避けたいです。

実際には下記のようにするのが良いでしょう。

// 擬似コード
// 実際には配列で保持されている
R_hat = S * pow(β, n) + A_0;
while(R_hat < D)
{
    R_hat += B;
    Q_hat++;
}
R_hat -= D;

3n 桁 / 2n 桁の除算の図解

A/B を考える
|       A_2       |       A_1       |       A_0       | / |       B_1       |       B_0       |

Aは以下の和で表される
|               Q B_1               |                 |
|                 |               Q B_2               |
|                 |                 R                 |


[A_2, A_1]/B_1 を考える
|       A_2       |       A_1       | / |       B_1       |

[A_2, A_1] は以下の和
|             Q_hat B_1             |
|                 |        S        |

Q_hat, S を用いて A を表す
|             Q_hat B_1             |                 |
|                 |        S        |       A_0       |

Q_hat, R_hat が Q, R と同じ形をみたすように R_hat を定義したい。R_hat が負になる可能性は一旦無視する。
|             Q_hat B_1             |                 |
|                 |             Q_hat B_2             |
|                 |               R_hat               |

以下の方程式がたてられる
|        S        |       A_0       |
=
|             Q_hat B_2             |
|               R_hat               |

R_hat について解くと
|        S        |       A_0       |
|            -Q_hat B_2             |

2n 桁 / n 桁の除算の実践

2n 桁 / n 桁の除算 と 3n 桁 / 2n 桁の除算 ができるようになりました。10 進数で実践してみましょう。

\beta=10, DIV\_LIMIT = 2, A=14926421, B=7894 として計算します。

  • 14926421 / 7894
    • 上半分の計算: 149264 / 7894
      • 1492 / 78 → \widehat{Q}=19, S=10
      • \widehat{R} = 1064 - 19 \times 94 = 1064 - 1786 = -722
      • \widehat{Q} = \widehat{Q}-1, \widehat{R} += 7894
      • Q = 18, R=7172
    • 下半分の計算: 717221 / 7894
      • 7172 / 78 → \widehat{Q}=91, S=74
      • \widehat{R} = 7421 - 91 \times 94 = 7421 - 8554 = -1133
      • \widehat{Q} = \widehat{Q}-1, \widehat{R} += 7894
      • Q = 90, R=6761
  • Q=1890, R=6761

任意桁の除算

2n 桁 / n 桁の除算と 3n 桁 / 2n 桁の除算が再帰的に実行できるようになったので、任意の桁の除算を 2n 桁 / n 桁の除算に落とし込むことで計算できるようにします。

  1. A, B を正規化する

B の桁数を計算しやすいように増やします。計算しやすい桁数 n は下記のように求めます。

k を 2^k \times DIV\_LIMIT \gt s をみたす最小の整数とします。また、m=2^kとします。

s 以上で最小の m の倍数を n とします。

\begin{aligned} m &= min\{2^k|2^k \times DIV\_LIMIT \gt s\} \\ j &= \lceil s/m \rceil \\ n &= j \times m \end{aligned}

B を \sigma = n - s 桁だけ大きくするため、B を B\times \beta^\sigma で置き換えます。基数 \beta が \beta=2 ならばビットシフトで補正可能です。

同様に A も A\times \beta^\sigma で置き換えます。

※各変数の解説

2^k \times DIV\_LIMIT \gt s は「愚直計算に移行するまでに k 回分割できる長さ」を表します。

Burnikel-Ziegler 法では除数の桁数を半分ずつに分割していくので、桁数の素因数に 2 を多く含むと都合が良いです。

DIV\_LIMIT = 5 での例を挙げてみます。

除数の桁数が 24, 26, 28, 32 の場合をそれぞれ考えると

  • 24 桁で除算
    • 分割して 12 桁で除算
    • 分割して 6 桁で除算
    • 分割して 3 桁で除算: 奇数なので/ DIV\_LIMIT を下回るので愚直計算
  • 26 桁で除算
    • 分割して 13 桁で除算: 奇数なので愚直計算
  • 28 桁で除算
    • 分割して 14 桁で除算
    • 分割して 7 桁で除算: 奇数なので愚直計算
  • 32 桁で除算
    • 分割して 16 桁で除算
    • 分割して 8 桁で除算
    • 分割して 4 桁で除算: DIV\_LIMIT を下回るので愚直計算

と桁数が 8 の倍数のときには、愚直計算する桁数が DIV\_LIMIT を下回ります。

つまり m=8 とすると、愚直計算する桁数が DIV\_LIMIT を下回るということになります。

下記の表を見ると、5 未満のとき 1　の倍数、10 未満のとき 2 の倍数、20 未満のとき 4 の倍数、40 未満のとき 8 の倍数、...　というように m が推移していくのが読み取れます。

除数の桁数 愚直計算に移行する桁数
2 2
4 4
6 3
8 4
10 5
12 3
14 7
16 4
18 9
20 5
22 11
24 3
26 13
28 7
30 15
32 4
34 17
36 9
38 19
40 5
42 21
44 22
46 23
48 3
50 25
  1. A を分割

以下のような t を求めます。

  • A \lt \dfrac{\beta^{n}}{2} ならば t=2 とします
  • そうでなければ、A \lt \dfrac{\beta^{l n}}{2} をみたす最小の整数 l について t=l とします

A を n 桁ごとに分割して [A_{t-1}, ..., A_0] と置きます。A_{t-1} については桁数が n 未満となります。

A_{t-1} の桁数が n 未満となるように t を置くのは 2n 桁 / n 桁の除算で A\lt B \beta^n という条件を満たすためです。

  1. 計算する

上位の桁から順にあまりを求めていきます。

// 擬似コード
// 実際には配列で保持されている
Z[t-2] = [A[t-1], A[t-2]];
Q[t-2] = Z[t-2] / B;
R[t-2] = Z[t-2] % B;
for(int i = t-3; i>=0; i--)
{
    Z[i] = [R[i+1], A[i]];
    Q[i] = Z[i] / B;
    R[i] = Z[i] % B;
}

A/B を考える
|     A_{t-1}     |     A_{t-2}     |     A_{t-3}     |       ...       |       A_0       | / |        B        |

A を複数の和で表すと下記のようになる
|     A_{t-1}     |     A_{t-2}     |     A_{t-3}     |       ...       |       A_0       |

=
|             B Q_{t-2}             |                                                     |
|                 |     R_{t-2}     |     A_{t-3}     |       ...       |       A_0       |

=
|             B Q_{t-2}             |                                                     |
|                 |             B Q_{t-3}                                                 |
|                                                     |       R_1       |       A_0       |

=
|             B Q_{t-2}             |                                                     |
|                 |             B Q_{t-3}                                                 |
|                                                     |               B Q_0               |
|                                                                       |       R_0       |

求める商は Q=[Q_{t-2}, ..., Q_0]、あまりはR_0 になります。

ただし、A, B を \beta^\sigma だけ補正していることから、真に求めるあまり R は R_0/\beta^\sigma です。

任意桁の除算の実践

任意桁除算を 10 進数で実践してみましょう。

\beta=10, DIV\_LIMIT = 2, A=57543907443, B=532 として計算します。

  • 57543907443 / 532
  • B の桁数を補正: A=575439074430, B=5320
  • A を分割: A=[0000, 5754, 3907, 4430] 最上位の桁数が B と同じにならないようにすることに注意。
    • [0000, 5754] / 5320 を計算
      • 上半分の計算: 000057 / 5320
        • 0000 / 53 → \widehat{Q}=0, S=0
        • \widehat{R} = 0057 - 0 \times 20 = 57 - 0 = 57
        • Q = 0, R=57
      • 下半分の計算: 005754 / 5320
        • 0057 / 53 → \widehat{Q}=1, S=4
        • \widehat{R} = 0454 - 1 \times 20 = 454 - 20 = 434
        • Q = 1, R=434
      • Q_2 = [00, 01] = 1, R_2 = 434
    • [434, 3907] / 5320 を計算
      • 上半分の計算: 043439 / 5320
        • 0434 / 53 → \widehat{Q}=8, S=10
        • \widehat{R} = 1039 - 8 \times 20 = 1039 - 160 = 879
        • Q = 8, R=879
      • 下半分の計算: 087907 / 5320
        • 0879 / 53 → \widehat{Q}=16, S=31
        • \widehat{R} = 3107 - 16 \times 20 = 3107 - 320 = 2787
        • Q = 16, R=2787
      • Q_1 = [08, 16] = 816, R_1 = 2787
    • [2787, 4430] / 5320 を計算
      • 上半分の計算: 278744 / 5320
        • 2787 / 53 → \widehat{Q}=52, S=31
        • \widehat{R} = 3144 - 52 \times 20 = 3144 - 1040 = 2104
        • Q = 52, R=2104
      • 下半分の計算: 210430 / 5320
        • 2104 / 53 → \widehat{Q}=39, S=37
        • \widehat{R} = 3730 - 39 \times 20 = 3730 - 780 = 2950
        • Q = 39, R=2950
      • Q_0 = [52, 39] = 5239, R_0 = 2950
  • Q=[Q_2, Q_1, Q_0] = [1,0816,5239]=108165239, R=R_0=2950
  • R の補正を外す: Q=108165239, R=295
GitHubで編集を提案

Discussion