アルゴリズム論特論(塩田)第10回 (2) 多倍長整数の平方根

多倍長整数の平方根

 多倍長整数の平方根は倍精度実数の sqrt では求められませんので、 「開平法」と呼ばれる方法を用いたり、 次のニュートン法を用いて計算します。
Algorithm 5(ニュートン法) 
  • 入力:自然数 $a$
  • 出力:$a$ の平方根の切り捨て
def mysqrt(a):
  if a < = 1:
    return a
  x = a
  while x * x > a:
    x = (x + a//x) // 2
  return x
 数値解析で実数のニュートン法を習いましたね。方程式 $f(x)=0$ の近似解は数列
$\displaystyle{x_{k+1}=x_k - \frac{x_k}{f'(x_k)}}$
を作って求めました。$\sqrt{a}$ は $f(x)=x^2-a=0$ の解で、そのとき上の数列は
$\displaystyle{x_{k+1}=x_k + \frac{x_k^2-a}{2x_k}=\frac{1}{2}\left(x_k+\frac{a}{x_k}\right)}$
となります。除算を整数商に代えたものが Alg.5 です。
Ex.6 
a = 9956894514697920852 ( 64 bit )

Step   0 : x = 9956894514697920852
Step   1 : x = 4978447257348960426
Step   2 : x = 2489223628674480214
Step   3 : x = 1244611814337240108
Step   4 : x = 622305907168620057
Step   5 : x = 311152953584310036
Step   6 : x = 155576476792155033
Step   7 : x = 77788238396077548
Step   8 : x = 38894119198038837
Step   9 : x = 19447059599019546
Step  10 : x = 9723529799510028
Step  11 : x = 4861764899755525
Step  12 : x = 2430882449878786
Step  13 : x = 1215441224941440
Step  14 : x = 607720612474815
Step  15 : x = 303860306245599
Step  16 : x = 151930153139183
Step  17 : x = 75965076602359
Step  18 : x = 37982538366715
Step  19 : x = 18991269314429
Step  20 : x = 9495634919358
Step  21 : x = 4747817983966
Step  22 : x = 2373910040558
Step  23 : x = 1186957117429
Step  24 : x = 593482753008
Step  25 : x = 296749765032
Step  26 : x = 148391659099
Step  27 : x = 74229378923
Step  28 : x = 37181757882
Step  29 : x = 18724773827
Step  30 : x = 9628261807
Step  31 : x = 5331196948
Step  32 : x = 3599431402
Step  33 : x = 3182836172
Step  34 : x = 3155572503
Step  35 : x = 3155454726
Step  36 : x = 3155454723

answer = 3155454723
Prop.7 Alg.5 の計算量は $O(\log^3 n)$.
証明 終盤までは $x$ がほとんど半分半分になっていきますので ( Ex.6 参照 )、ループ回数は $O(\log n)$ で、 ループの中の計算は乗算と除算1回ずつゆえ。(証明終)