多倍長整数の割り算
1.定義
非除数(割られる数)
123454322の場合
除数
11111の場合
商
余り
基数(進数)
D=1桁の最大値+1 10進数の場合は10関係式
2.計算を効率化するために非除数・除数に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)式の右辺の数値の範囲は
したがって
3.3.さらに仮商の誤差を少なくする
3.3.1.手法
前ページでは仮商と真の商は2以上の誤差がないことが分かったが、さらに精度を上げ誤差を1以下とする方法を探る前ページの仮商を計算後、非除数上位3桁目および除数の上位2桁目を使用して誤差が2のケースを削除する。さらに誤差が1のケースを少なくする。
上記式が成立する場合q'から1を引くr'にb_{t}を加える
3.3.2.非除数3桁÷除数2桁がq<q'であることの証明
3.3.3.非除数3桁÷除数2桁がq≦q'+1の証明
3.4.仮商より1桁の商を算定。次の桁の非除数の上位桁つまり今回の余りを算定する。
余りは、次の桁の商を求めるときr'が上位桁となるので、uが現在の商の桁位置とすると、r'をa[u]に代入すればよい。
もっと簡単には、
商はq[u]=q'
3.5.次下位の桁へ移動する
商の位置uが最下位でない限り、一つ下位の桁へ移動する --u3.項へジャンプし処理を繰り返す。
最下位の商が定まれば、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,);
}