世界は2乗でできている 自然にひそむ平方数の不思議 (ブルーバックス)
- 作者: 小島寛之
- 出版社/メーカー: 講談社
- 発売日: 2013/08/20
- メディア: 新書
- この商品を含むブログ (13件) を見る
この本に、オイラーのバーゼル問題に対する執念が書かれていました。バーゼル問題とは、
S = 1 + 1/22 + 1/32 + …
を求めよ、というもので、π2/6となります。これを求めるために、オイラーは小数点以下20桁まで求めたというのです。
これを何の気なしに読むと、二乗の逆数を延々と計算し続けたのかと思ってしまうかもしれませんが、そんなことはありません。なぜなら、途中の第n項で打ち切ると、真値に対する差は、
と近似できて、これはおよそ1/nです。すなわち、1020程度まで計算を行わなければ小数点以下20桁までは求められないということになります。
ではEulerはどのように20桁まで求めたのでしょうか。それには上の近似を使うことが考えられます。
これももちろん誤差があります。これを評価しましょう。誤差は、
2000000くらいまで計算すれば20桁の精度が得られるようです。まだまだですね。
このように次々と近似をよくすることができます。この近似を一般化しましょう。
を近似します。
また、
と置くと、
と近似を進めることができます。本当は一般的には近似が進むわけではないのですが、その辺には目をつむります。
これを繰り返し適用すると、
となります。でもう一度計算すると、
ただ、f2を求めるときにf1を積分するとlogが出てきてややこしくなるので、Taylor展開をすることにします。
この計算を次々に行うのはオイラーならやるかもしれませんが、我々には計算機代数という武器があるのでそれを使いましょう。
を、Pythonで[ 0, a1, a2, … ]と表します。Taylor展開は適当な項で打ち切ります。結果は、
これを使ってSの近似式を書くと、
n=100を代入してDecimalモジュールを使って高精度に計算すると、
3.1415926535 8979323846 2643383279 5028841971 699258655
と真値
3.1415926535 8979323846 2643383279 5028841971 693993751
に対して小数点以下42桁まで一致しました。
実際にオイラーがどのように計算したかわかりませんが、オイラーならこれくらい手で計算できそうな気がします。
from itertools import * from fractions import Fraction from decimal import Decimal, getcontext import time getcontext().prec = 50 def C(n, m): return reduce(lambda x, y: x * (n - y + 1) / y, range(1, m + 1), 1) def integral(f): return [ 0 ] + [ -c / (n - 1) for n, c in enumerate(f[2:], 2) ] + [ 0 ] def subtract(f, g): return [ a - b for a, b in izip(f, g) ] def add(f, g): return [ a + b for a, b in izip(f, g) ] # 1/(x-a)^n -> 1/x^n + na/x^{n+1} + n(n+1)/2a^2/x^{n+2} + ... def expand(a, n): return [ 0 ] * n + [ C(k - 1, k - n) * a ** (k - n) for k in range(n, M + 1) ] # f(x) -> f(x + a) def shift(f, a): g = [ 0 ] * (M + 1) for k in range(1, M + 1): for l, c in enumerate(expand(-a, k)): g[l] += c * f[k] return g def next_f(f): g = integral(f) half = Fraction(1, 2) return subtract(f, subtract(shift(g, half), shift(g, -half))) def compute_fs(N): f0 = [ 0, 0, Fraction(1) ] + [ 0 ] * (M - 2) fs = [ f0 ] f = f0 for _ in range(N): f = next_f(f) fs.append(f) return fs def min_dim(f): for d, c in enumerate(f): if c != 0: return d def TeX_term(c, dim, var): d = c.denominator n = c.numerator xd = (var + "^{%d}") % dim if dim > 1 else var if d == 1: return r"\frac{%d}{%s}" % (n, xd) elif n > 0: return r"\frac{%d}{%d%s}" % (n, d, xd) elif n < 0: return r"-\frac{%d}{%d%s}" % (-n, d, xd) def TeX(f, var = "x"): min_d = min_dim(f) s = "" head = True for d, c in enumerate(f[min_d:], min_d): if c > 0 and not head: s += "+" + TeX_term(c, d, var) elif c != 0: s += TeX_term(c, d, var) head = False return s def frac_to_dec(f): return Decimal(f.numerator) / f.denominator def value(f, x): return reduce(lambda p, q: p / x + q, map(frac_to_dec, reversed(f)), 0) def approx_S(N): L = 100 fs = compute_fs(N) for k, f in enumerate(fs): print "[tex: f_{%d}(x) = %s]" % (k, TeX(f)) g = reduce(lambda x, y: subtract(x, integral(y)), fs, [ 0 ] * (M + 1)) print TeX(g, "(n+1/2)") res = value(g, L + 1 / Decimal(2)) return sum(1 / Decimal(n * n) for n in range(1, L + 1)) + res t0 = time.clock() N = 10 M = N * 2 + 2 S = approx_S(N) print (S * 6).sqrt() print time.clock() - t0