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: プログラムメモ

前回、平方根 √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: プログラムメモ