無理数に見えるのに整数になる？
このページはマス旅の一部です。
ラマヌジャンというと、無限に続く式というイメージがありますね。
一見、無理数に見えるけれど、計算すると整数になる式。
そんな式からはじめてみよう。
１．数式とウォーミングアップ
sqrt(1+2 sqrt(1+ 3 sqrt(1+ 4 sqrt( 1+ 5 sqrt(1+ ...... ＝?
geogebraのTexでかくと、
？
無理数の計算をすれば無理数になるんじゃないか？
と思いがちですが整数になります。
では、ウォーミングアップから入りましょう。
＜平方根の連分数展開＞
入れ子の分数の極限が無理数になるというのは知ってますよね。
1+1/(2+1/(2+1/(2+1/......)))=√2
geogebraのTexでかくと、
√２
√2となる結果をわからないとしてｘとおきましょう。
すると、x=１＋1/□
という形に見えますよね。
その分母□のところをじっと見よう。
□＝２＋１/△になりますね。
ここで、□と△を見比べてみよう。無限があることから□と△が同じになってしまう。
ということは、
最初の式の分母□＝２＋1/□＝１＋(１＋1/□)=1+ｘ。
つまり、最初の式全体は
ｘ＝１＋1/(1+x)
となるね。
x-1=1/(1+x)
(x-1)(x+1)=1
x^2=2
x＝√2（正だから）
＜分数が√にかわると。。。＞
視界をずらす練習ができたところで、次にいこう。
？
これも同じ作戦でいこう。
式の結果をｘとすると、
ｘ＝√(2+□)と見える。
今度は□を見ると、□は無限性からｘに見えるね。
ｘ＝√(2+x)となるから、
x^2=2+x
x^2-x-2=(x-2)(x+1)=0
x=2（正だから）
２．数式の値を求めよう
さあいよいよ、
sqrt(1+2 sqrt(1+ 3 sqrt(1+ 4 sqrt( 1+ 5 sqrt(1+ ...... ＝?だ、
これまでの数式では一定数だったのが、sqrtの前が1,2,3,4,5と変化しているので、
とおいて、漸化式のようなものが作れたらいいですね。
そこで、計算結果をｘとおくのではなく、
平方根にかける数列をx,x+1,x+2,...とするときのｘの関数値になるとしてf(x)とおくのです。
たまたま、最初の問題ではｘ＝２にしていただけだということにするのですね。
もし、だから、計算結果はf(2)になるでしょう。
さて、
f(x)=sqrt(1+x sqrt(1+ (x+1) sqrt(1+ (x+2) sqrt( 1+ (x+3) sqrt(1+ ......
すると、
[f(x)^2-1]/x=sqrt(1+ (x+1) sqrt(1+ (x+2) sqrt( 1+ (x+3) sqrt(1+ ......
となるね。
右辺はf(x)のｘが１つ進んだだけなので、f(x+1)だ。
f(x)^2=1+ xf(x+1)
一見ややこしそうにみえるが、f(x)が1次式だとすると、左辺も右辺も2次式になる。
そして、左辺がf(x)の2乗の典型の定数項が、右辺で１になることから、f(x)=x+1が予想できるね。
まさに、(x+1)^2=1+x(x+1+1) が成り立つ。
もう解けたと同じです。
計算結果はf(２)＝２＋１＝３。
＜f(2)の収束の速さ＞
このf(2)の値をもし、推定するとしたら、どのくらいで３に近くなるでしょうか。
a1=√(1+2*0)=1
a2=√(1+2√(1+3*0)=√3=1.73205...
a3=√(1+2√(1+3√(1+4*0)))=√5=2.23606...
a4=√(1+2√(1+3√(1+4√(1+5*0)))=....=2.559...
結構、スピードが速そうですね。
akは、根号がｋ個あり、最後の根号の中は１+0になりますね。
根号の最後（内側）から前（外側）に向かってする1回の操作を関数にしましょう。
直前の値をprevとすると、次の結果は
ret=sqrt(1+k*prev)になります。この関数をnext(k,prev)としましょう。
たとえば、k倍はk=4の場合2,3,4では足りないので2,3,4,5と4回まわす形にしたいですね。
そこで、a4を出すには最初のprev=0にして、k=5,4,3,2とカウントダウンしながらnext(k,prev)を
まわせばよいですね。
そこで、
ループの変数をi として、range(n)とすると
result = next(n+1-i, result)にすると、
n=4のときは、iは0から3まで動き、n+1-iは4+1-0=5から4+1-3=2まで変化します。
これで、a4=√(1+2√(1+3√(1+4√(1+5*0)))と同じ計算ができることが、脳内デバッグできるでしょう。
実際まわしてみると、n=10でほぼ３になるという猛スピードで収束することがわかるね。
課題：どのくらいで３になるかをPythonで視覚化しよう。
from math import sqrt
import matplotlib.pyplot as plt
def next(k,prev):
ret =sqrt(1 + k * prev)
if ret<0:
ret=0
return ret
def nested(n):
result = 0
for i in range(n):
result = next(n+1-i, result)
return result
# ==========================================
max_depth = 15
# ==========================================
# 視覚化してみる：深さnに対する値の収束プロット
depths = range(1, max_depth + 1)
values = [nested(d) for d in depths]
plt.figure(figsize=(10, 6))
plt.plot(depths, values, marker='o', linestyle='-', color='#674ea7', markerfacecolor='#e06666')
plt.axhline(y=3, color='black', linestyle='--', label='convergence: 3')
plt.xlabel('depth(n)')
plt.ylabel('value')
plt.title('where to (√1+2√1+3√...)')
plt.legend()
plt.grid(True)
plt.show()
[OUT]
入れ子の根号式の収束
＜振り返り＞
ラマヌジャンの発想を想像してみよう。
以下は、zenの勝手な想像であることを承知して読んでください。
ラマヌジャンにとっての数式は全体としてある１つの生き物であり、
登場する数はその生き物の状態を変えられるパラメータだ。
しかもそれが、ちまちま変えるのでなく、
一斉にブロードキャストできるものだ。
つまり、こんなイメージ
f(2, 3,4,5,.....)＝sqrt(1+2 sqrt(1+ 3 sqrt(1+ 4 sqrt( 1+ 5 sqrt(1+ ......=???
f(1, 2,3,4,.....)＝sqrt(1+1 sqrt(1+ 2 sqrt(1+ 3 sqrt( 1+ 4 sqrt(1+ ......=??
f(0, 1,2,3,.....)＝sqrt(1+0 sqrt(1+ 1 sqrt(1+ 2 sqrt( 1+ 3 sqrt(1+ ......=sqrt(1+0 )=1
f(-1, 0,1,2,.....)＝sqrt(1-1 sqrt(1+ 0 sqrt(1+ 1 sqrt( 1+ 2 sqrt(1+ ......=sqrt(1-1 sqrt(1+0))=0
そこで、sqrtの前の数を一斉に変化させるために
f(2,3,4,5,.....)をf(2)とすると、
f(0)=1, f(-1)=0だから、f(x)=x+1が予測できる。
ラマヌジャン自身は証明しなくても自明だったのではないでしょうか。
証明をしない数学者は、
「証明するまでもなく自明な筋道」を見つけているということだと思います。
課題：geogebraでnextを回す感じで3に近づくことを確かめよう。
タイトルは「ラマヌジャンの入れ子の根号計算」
g(k,x)=sqrt(1+ k x)
a10=g(2,g(3,g(4,g(5,g(6,g(7,g(8,g(9,g(10,g(11,0))))))))))
#2.98992
a11=g(2,g(3,g(4,g(5,g(6,g(7,g(8,g(9,g(10,g(11,g(12,0)))))))))))
#2.9948
a12=g(2,g(3,g(4,g(5,g(6,g(7,g(8,g(9,g(10,g(11,g(12,g(13,0))))))))))))
#2.99733
IterationListを使わなくても、関数を入れ子にするだけで、済みますね。
1回の処理を１つの関数としておくことで、
関数の入れ子の式をそのまま再現していることもわかるし、
あっけなく計算できてしまうね。
ということは、Pythonでも同様な関数の入れ子ができるはずだね。
# 1回分のルートの皮を被せる関数
def g(k, x):
return sqrt(1 + k * x)
# 数式をそのままコードに翻訳する
a10 = g(2, g(3, g(4, g(5, g(6, g(7, g(8, g(9, g(10, g(11, 0))))))))))
print(a10) # 2.98992...