Hundredsoft

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

Bigint newtonsqrt(const Bigint& q, const int& n, const Bigint& init){
  Bigint x(init), m(0), c3(3);
  c3 <<= n;
  while(m != x){
    m = x;
    x *= c3 - ((m * m * q) >> n);
    x >>= n+1;
  }
  return x;
}

Bigint iSqrt(Bigint s){
  int n = s.length(); // bit数
  int b;
  if (n & 1)    b = 65; // 初期演算Bit数(n:奇数)
  else          b = 64; // 初期演算Bit数(n:偶数)

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

  Bigint w;
  int m = n / 8;  // 4分割程度

  while(1){
    w = s >> (n - b);        // (s)先頭から指定bit数(n-b)切り取り
    x = newtonsqrt(w, b, x);

    if (b + 2 * m > n)
      break;

    b += 2 * m;
    x <<= m;
  }
  x <<= (n - b) / 2;         // 必ず割り切れる
  x = newtonsqrt(s, n, x);

  x *= s;         // x = s * (2^n) / √s = (2^n) * √s
  x >>= n;        // x = √s

  // 誤差修正
  w = x + 1;
  if (s >= w * w){
     x = w;
  }
  return x;
}

newtonsqrtは、Newton法で （2ⁿ / √q） を求める関数です。
初回は先頭64bit程度で計算します。
2回目は、初回の計算結果を初期値として、(n/8)桁増やして計算。
3回目は、２回目の計算結果を初期値として、更に(n/8)桁増やして計算。
と次々と桁を増やしていきます。

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

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

また、この方法は平方根に限らず、Newton法を使った様々な計算に応用できます。

整数平方根での詳細解説

整数除算での詳細解説

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


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

---

Tags: プログラムメモ