アルゴリズム論特論（塩田）第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)$ で、
ループの中の計算は乗算と除算１回ずつゆえ。(証明終)