15.5.1 平方根のアルゴリズム

平方根は，Paul Zimmermann (see References)が提唱した“Karatsuba 平方根” アルゴリズムを使用します。 入力値nをk bit単位で4つに分割し，b=2^kを用いて，n = a3*b^3 + a2*b^2 + a1*b + a0とします。a3の部分は，最大bitもしくはその次のbitが1になるように“正規化”されていなければなりません。 GMPでは, kはリムの境界上に置かれ，入力値は偶数ビット数分左シフトして正規化します。

上位から二つ分の平方根は下記の反復を適用することで得られます(1リム用のNewton法で収束させる）。

s1,r1 = sqrtrem (a3*b + a2)

これによって平方根の近似値が得られるので，s,rを求めるための除算を行って平方根の精度を上げていきます。

q,u = divrem (r1*b + a1, 2*s1)
s = s1*b + q
r = u*b + a0 - q^2

a3を正規化する理由は，この時点においては，sは正しい値か1だけ大きい値になっているからです。後者の場合はrは負数になっているので，下記のように補正します。

if r < 0 then
  r = r + 2*s - 1
  s = s - 1

このアルゴリズムは，分割統治法のスタイルで表現されていますが，論文で指摘しているように，Newon法の離散バージョンとしてみることも，一度に2桁平方根を出していくスクールボーイ法(いちいち説明しねぇよ）の変種と見ることもできます。

剰余rが不要の場合は，r と uの上の数リムだけ求めて，sの調整量が必要かどうかを見るだけに使用できますが，このような最適化は今のところ実装していません。

Karatsuba 乗算の範囲内では，このアルゴリズムはO(1.5*M(N/2))となります。ここで M(n)はnリム長の2数の乗算の計算時間です。FFT乗算の範囲内では， O(6*M(N/2))で抑えられます。実際には，Karatsuba乗算と3方向Toom-Cookでは1.5～1.8倍ぐらいですが，FFT乗算の範囲ではこれが2～3倍に増えます。

このアルゴリズムは全て整数演算で行います。mpn_sqrtrem関数の結果はmpz_sqrt関数やmpf_sqrt関数で利用されています。 mpf_sqrt_ui関数では，ゼロリムをパディングすることで長い桁の平方根を導出しています。