多倍長整数の割り算

1.定義

非除数(割られる数)

a=a_{s}D^s+a_{s-1}D^{s-1}+a_{s-2}D^{s-2} \, \cdot \cdot \cdot \, +a_{0}
123454322の場合a=123454322,s=8,a_{8}=1,a_{7}=2,a_{6}=3,a_{5}=4,a_{4}=5,a_{3}=4,a_{2}=3,a_{1}=2,a_{0}=2となる

除数

b=b_{t}D^t+b_{t-1}D^{t-1}+b_{t-2}D^{t-2} \, \cdot \cdot \cdot \, +b_{0}
11111の場合b=11111,t=4,b_{4}=1,b_{3}=1,b_{2}=1,b_{1}=1,b_{0}=1となる

商

q=q_{u}D^u+q_{u-1}D^{u-1}+q_{u-2}D^{u-2} \, \cdot \cdot \cdot \, +q_{0}

余り

r=r_{v}D^v+r_{v-1}D^{v-1}+r_{v-2}D^{v-2} \, \cdot \cdot \cdot \, +r_{0}

基数(進数)

D=1桁の最大値+1 10進数の場合は10

関係式

a=qb+r \, (0 \leq r<b)

2.計算を効率化するために非除数・除数にkを掛ける

全桁の割り算を実行するのは困難であるため、上位桁同士の割り算を実行し商を推測する
除数にkを掛けた結果、除数の最上位の値が\frac{D}{2} \leq b_{t} \cdot k \lt D になるようにkを決める
この結果btが非常に小さい場合、商の計算の効率が悪いのでkを掛けることにより商の候補が絞りこみやすくなる
k=\frac{D}{b_{t}+1}
k=\frac{10}{1+1}=5
a,bにkを掛ける
a=a*k=123454322*5=617271610,s=8,a_{8}=6,a_{7}=1,a_{6}=7,a_{5}=2,a_{4}=7,a_{3}=1,a_{2}=6,a_{1}=1,a_{0}=5 b=b*k=11111*5=55555,t=4,b_{4}=5,b_{3}=5,b_{2}=5,b_{1}=5,b_{0}=5

3.商の最上位桁位置及びu及び商の桁値quを求める

3.1.as≧bt 非除数の最上位桁の値≧除数の最上位桁の値 である場合

a≧bD^(s-t)の時 非除数≧除数の最上位をs-t桁
kを掛けていることにより
\frac{D}{2}\leq b_{t} , a_{s} < D
であることから、Dを消去すると
a_{s} < 2b_{t}
\frac{a_{s}}{b_{t}}<2
\frac{a_{s}}{b_{t}}\geq1
となるので最上位の商は1となる。
商の最上位桁の位置は
u=s-t
u=8-4=4
最上位桁の商の値は
q_{u}=1

3.2.as≧bt 非除数の最上位桁の値≧除数の最上位桁位の値である場合

a < bD^(s-t)の時 非除数の最上位桁<除数の最上位桁の場合(非除数2桁/除数1桁で計算)
商の最上位の位置はu=s-t-1となる。
\frac{ (a_{s}D+a_{s-1})D^{s-1} }{b} \ \leq \ q \ < \frac{ (a_{s}D+a_{s-1})D^{s-1}+D^{s-1}}{b} \ \cdot \cdot \cdot (1)
\frac{ (a_{s}D+a_{s-1})D^{s-1} }{b} \ \leq \ q' \ \leq \ \frac{a}{b_{t}D^{s-1}} \ \cdot \cdot \cdot (2)

(a_{s}D+a_{s-1})D^{s-1} \ \leq \ a \ < \ \{ a_{s}D+(a_{s-1}+1) \}D^{s-1} \ \cdot \cdot \cdot (1')
q'=\frac{ (a_{s}D+a_{s-1})D^{s-1} }{b_{t}D^{s-1}} \ \cdot \cdot \cdot (2')
分子 \leq u, 分母 \leq v
a=qb+r \ , \ r < b
a < qb + b \ , \ q > \frac{a}{b}-1
(1)(2)式より不等式を作成すると
q'-q \ >\ \{ \frac{ (a_{s}D+a_{s-1})D^{s-1} }{b} - \frac{ (a_{s}D+a_{s-1})D^{s-1}+D^{s-1} }{b} \}
q'-q \ >\ -\frac{D^{s-1}}{b} \ \geq \ -1
q'-q \ \leq \ \frac{a}{b_{t}D^{s-1}} - (\frac{a}{b}-1)
=\frac{ ab - ab_{t}D^{s-1} }{bb_{t}D^{s-1}}+1
=\frac{a}{b} \cdot \frac{b-b_{t}D^{s-1} }{b_{t}D^{s-1}}+1
< \frac{a}{b} \cdot \frac{D^{s-1} }{b_{t}D^{s-1}}+1
=\frac{a}{b} \cdot \frac{1}{b_{t}}+1
\frac{a}{b} \ < \ D
-1 \ <\ q'-q \ <\ \frac{D}{b_{t}}+1 \ \cdot \cdot \cdot (3)
\frac{D}{2} \leq b_{t} < Dならば
(3)式の右辺の数値の範囲は\frac{2D}{D}=2 \ , \ \frac{D}{D}=1
0 \ \leq \ q'-q \ \leq 2
したがってq=\frac{a}{b}\ ,\ q'=aの上位2桁 / bの上位1桁 とすると qはq'と等しいか(q'+1)または(q'+2)のいずれかである。

3.3.さらに仮商の誤差を少なくする

3.3.1.手法
前ページでは仮商と真の商は2以上の誤差がないことが分かったが、さらに精度を上げ誤差を1以下とする方法を探る
前ページの仮商を計算後、非除数上位3桁目および除数の上位2桁目を使用して誤差が2のケースを削除する。さらに誤差が1のケースを少なくする。
q'=D \ \ またはq' b_{t-1} > D r'+a_{s-2}であるかチェックする
上記式が成立する場合q'から1を引くr'にb_{t}を加える
r' < Dならチェックを繰り返す
3.3.2.非除数3桁÷除数2桁がq<q'であることの証明

a=bq+r
a-q'b \leq a-q'b_{t} D^{s-1}-q'b_{t-1}D^{s-2}
=a_{s}D^{s}+a_{s-1}D^{s-1} +a_{s-2}D^{s-2} \ \cdot \cdot \cdot \ +a_{0} -q'b_{t} D^{s-1}-q'b_{t-1}D^{s-2} \ \cdot \cdot \cdot \ (1)
a_{s}D^s+a_{s-1}D^{s-1}-q'b_{t}D^{s-1}=r'D^{s-1}(1)式に左式の左辺を代入
<a_{s-2}D^{s-2} \ \cdot \cdot \cdot \ +a_{0} + r'D^{s-1} -q'b_{t-1}D^{s-2}
<D^{s-2} \{ a_{s-2} + r'D -q'b_{t-1} \} \leq 0
a-q'b < 0 \ \, \ q<q'
3.3.3.非除数3桁÷除数2桁がq≦q'+1の証明
q \leq q' -2 なら
a < (q'-1)b < q' \{ b_{t} D^{s-1} +(b_{t-1}+1) D^{s-2} \}-b
< b_{t}D^{s-1}+ \{ Dr'+b_{t-1} \} D^{s-1}+D^{s-2}-b
=a_{s}D^{s} + a_{s-1}D^{s-1}+a_{s-2}D^{s-2}+D^{s-1}-b
\leq a_{s}D^{s} + a_{s-1} D^{s-1} +a_{s-2} D^{s-2} \leq a
a<aとなることから矛盾であるq=q'+2のケースはない
3.4.仮商より1桁の商を算定。次の桁の非除数の上位桁つまり今回の余りを算定する。
余りは、r'=a-q'bとなる。
次の桁の商を求めるときr'が上位桁となるので、uが現在の商の桁位置とすると、r'をa[u]に代入すればよい。
もっと簡単には、a=a-q'bD^{u-1}を算定すればよい。ただしaはs～jの桁を対象とすればよい。もし計算の結果aがマイナスとなった場合は、q'が大きすぎるのでq'から1を減ずる。またaにを加える
商はq[u]=q'
3.5.次下位の桁へ移動する
商の位置uが最下位でない限り、一つ下位の桁へ移動する --u
3.項へジャンプし処理を繰り返す。
最下位の商が定まれば、aに余りが得られる
ただしこの余りは、k倍されているのでaをkで割ると正規化前の余りが求められる
123454322/11111
450000000/59999
400,000,000,000/500,999,999(1000進数の場合)

アルゴリズム

// 多倍長整数同士の割算
// a/b=q r:余り
void div(int* a,int* b,int* q,int* r,int s,int t){
  D=10;
  // 正規化
  k=D/(b[t]);
  a=lmul(a,k); // a=a*k
  b=lmul(b,k); // b=b*k
  int qq; // 仮商
  int rr; // 仮余
  for(;;){
    // 仮商qqを求める
    if(a[s]==b[t])
      qq=D-1;
    }
    qq=(a[s]*D+a[s-1])/b[t];
    rr=(a[s]*D+a[s-1])%b[t];
    while(b[t-1]*qq > D*rr+a[s-2]){
      --qq;
      rr+=b[t-1];
      if(D<=rr)
        break;
    }
    // a-仮商(qq)*bを求める 次の桁の仮商を求めるときの上位桁となる
    z=lmul(b,qq); // b*qq;
    y=lsub(a,z,);
  }