Hundredsoft

前回、除算 x = p / q　をNewton法を使って解くプログラムを紹介しましたが、
単独の初期値を使った方法では、あまり早くありません。
数万bitを超える多桁演算の場合、64bit精度の初期値を使ったとしても、あまり速度に差が出ません。
そこで大きな桁(数万bitを超える)で計算を開始するのではなく、除数(q)の先頭部分を使った小さな桁で計算を開始して、これをNewton法の初期値とし、もう少し大きな桁でNewton法を解く、
これを繰り返して最終的な答えを得る方法を紹介します。
少しコード修正しました(2011.2.11)

Bigint newtondiv(const Bigint& q, const int& n, const Bigint& init){
  Bigint x(init), m(0), c2(2);
  c2 <<= n;
  while(m != x){
    m = x;
    x *= c2 - q * x;
    x >>= n;
  }
  return x;
}

Bigint iDiv(Bigint p, Bigint q){
  int pb = p.length(); // 被除数のbit数
  int qb = q.length(); // 除数のbit数
  int n;
  int b = 64;
  int step = qb / 4;

  Bigint d = 1;
  d <<= b;      // Newton法の初期値として, 1 << b を与える。

  d = newtondiv(q >> (qb-b), b*2, d);  // 1回目
  b += step;
  d <<= step;

  d = newtondiv(q >> (qb-b), b*2, d);  // 2回目
  b += step;
  d <<= step;

  d = newtondiv(q >> (qb-b), b*2, d);  // 3回目
  d <<= (pb-b);
  n = pb + qb;

  d = newtondiv(q, n, d);  // 4回目

  d *= p;        // d = p / q
  d >>= n;

  // 誤差修正
  if (p >= (d+1)*q){
      d += 1;
  }
  return d;
}

newtondivは、Newton法で （2ⁿ / q） を求める関数です。
初回は先頭64bit程度で計算します。
2回目は、初回の計算結果を初期値として、64+(除数のBit数/4)での計算。
3回目は、２回目の計算結果を初期値として、64+(除数のBit数/2)での計算。
4回目は、３回目の計算結果を初期値として、全桁での計算。
と計算桁を増やしていきます。

除数・被除数のとり方にもよりますが、計算時間は明らかに短くなります。
この方法は、Bigintの演算回数としては、効果はないのですが、
可変長のFFT乗算を用いる場合、演算量は２の冪数桁単位に増加するので、
Totalの計算量としては、こちらの方が有利です。

4分割としていますが、利用するBigint(多倍長整数)クラスで最速となるよう調節する必要があります。
これについては、「高橋 大介, 金田 康正. "多倍長平方根の高速計算法". 情報処理学会研究報告 95-HPC-58, pp.51-56.」が詳しいです。
文献に従えば、初期値に倍精度結果を使う・３〜４分割程度・FFT再利用、を行えば最大限の高速化を望めますが、このサンプルコードではそこまで対応できていません。

ここでは説明を簡単にするため初期値に（2^64）を与えています。
より正確な初期値のとり方や、詳細アルゴリズムはこちら
整数除算の詳細解説

除数・非除数の桁差が小さければ、このような方法もあります。
回復法を使った２進多桁の除算

同様の計算方法を行っている、整数平方根解説はこちら
整数平方根での詳細解説

このアルゴリズムを利用した16進・２進・10進電卓です。
整数だけですがおよそ100万桁まで計算できます。


(注)ここでのBigintクラスは、多倍長整数を格納する仮想のクラスであり、実在するものではありません。

---

Tags: プログラムメモ