√2を高精度で計算する

昨日はcos20°を求めましたが、今日は√2を求めます。100万桁を目標にしましょう。
√nの循環連分数になります。√2は、

√2 - 1 = 1 / (2 + 1 / (2 + ...)) = [ 0; 2, 2, ... ]

となります。連分数を途中で打ち切ったときの分数をa_n / b_nと置くと、

a_n+2 = 2a_n+1 + a_n
b_n+2 = 2b_n+1 + b_n
a₀ = 1, a₁ = 0, b₀ = 0, b₁ = 1

さて、漸化式は次のように式を加えることで行列に変えることができます。

a_n+2 = 2a_n+1 + a_n
a_n+1 = a_n+1

$A = \left(\begin{array}2&1\\1&0\end{array}\right)$
$u_n = \left(\begin{array}a_{n+1}\\a_n\end{array}\right)$

とおくと、

u₁ = Au₀

行列になると結合則が成り立つから分割統治法が使えますが、すべて同じ係数なので単にべき乗で表すことができますね。

u_n = Aⁿu₀

さて、100万桁の精度を出すにはどれだけのべき乗を求めなければならないでしょうか。
Aの固有値を求めると、最大の固有値は1 + √2となります。ゆえに、Aⁿの各項は(1 + √2)ⁿ程度となります。また、a_n / b_nが√2 - 1の近似になるのですが、(a_n + b_n)² - b_n² = ±1となるので、

(a_n + b_n) / b_n = √2(1 - 1 / 2b_n²)

で精度が1 / 2b_n²程度となります。E乗だとするとより正確には、

4 / (1 + √2)^2E = 1 / 10^N

ここから何乗すればいいかわかります。
100万桁を求めるのに82秒かかりました。もうちょっと速いイメージだったのですが。これを10進表示しようとするとさらに10倍とかかかりそうです。

from math import sqrt, log
import time

def square(A):
    b = A[0][1] * A[0][1]
    a11 = A[0][0] * A[0][0] + b
    a12 =  A[0][1] * (A[0][0] + A[1][1])
    a22 = A[1][1] * A[1][1] + b
    return [ [ a11, a12 ], [ a12, a22 ] ]

def mul_mat(A, B):
    return [ [ sum(A[i][k] * B[k][j] for k in xrange(2))
                for j in xrange(2) ] for i in xrange(2) ]

def pow_mat(A, e):
    if e == 1:
        return A
    else:
        B = pow_mat(A, e / 2)
        if e % 2 != 0:
            return mul_mat(square(B), A)
        else:
            return square(B)

t0 = time.clock()
N = 10 ** 6
E = int((log(4) + N * log(10) / (2 * log(1 + sqrt(2))))) + 1
A = [ [ 2, 1 ], [ 1, 0 ] ]
B = pow_mat(A, E)
(B[0][1] + B[0][0]) * 10 ** N / B[0][0]
print time.clock() - t0

さて、昨日の問題も同様に速く、と考えたくなりますが、しかし、平方根のようにうまく連分数で表すこともできず、またいつものニュートン法でもあまりうまく分数で近似をすることができないようなのです。