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

多倍長整数の平方根

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

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

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