多倍長整数の割り算

1.定義

非除数(割られる数)


123454322の場合となる

除数


11111の場合となる


余り


基数(進数)

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

関係式


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

全桁の割り算を実行するのは困難であるため、上位桁同士の割り算を実行し商を推測する
除数にkを掛けた結果、除数の最上位の値が になるようにkを決める
この結果btが非常に小さい場合、商の計算の効率が悪いのでkを掛けることにより商の候補が絞りこみやすくなる


a,bにkを掛ける

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

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

a≧bD^(s-t)の時 非除数≧除数の最上位をs-t桁
kを掛けていることにより

であることから、Dを消去すると



となるので最上位の商は1となる。
商の最上位桁の位置は


最上位桁の商の値は

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

a < bD^(s-t)の時 非除数の最上位桁<除数の最上位桁の場合(非除数2桁/除数1桁で計算)
商の最上位の位置はu=s-t-1となる。





分子 , 分母


(1)(2)式より不等式を作成すると









ならば
(3)式の右辺の数値の範囲は

したがっての上位2桁 / bの上位1桁 とすると qはq'と等しいか(q'+1)または(q'+2)のいずれかである。

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

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




(1)式に左式の左辺を代入



3.3.3.非除数3桁÷除数2桁がq≦q'+1の証明
なら




となることから矛盾であるq=q'+2のケースはない
3.4.仮商より1桁の商を算定。次の桁の非除数の上位桁つまり今回の余りを算定する。
余りは、となる。
次の桁の商を求めるときr'が上位桁となるので、uが現在の商の桁位置とすると、r'をa[u]に代入すればよい。
もっと簡単には、を算定すればよい。ただし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,);
  }

実行ファイルとソースファイルのダウンロード