山本ワールド

多倍長整数の割り算

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