この記事では
前回の記事
に続いて
さて前回までの考察で以下のことまでわかっていたのでした。
が
が成り立つ。特に
前回の記事を書いてからというもの、どうにかここから評価を改善できないか試行錯誤したものですが何とも言えず仕舞いだったので、ここは一般の
見るからに
この予想の精密化と解決が今後の課題となりますが、私には思い付きそうもないのでこの記事の読者にでも委ねようと思います。
具体的に
手順としては
実際に使ったプログラムは以下のようである。
import decimal
d=750 #有効数字の桁数
decimal.getcontext().prec=d
m=42 #求めるkの最小値(やり直すとき用)
K=50 #求めるkの最大値
N=10000 #計算するnの上限値(n_kが膨大になったとき用)
def natti(A,k): #計算回数の少ない漸化式L_{n+1}=2L_n-L{n-k}の適用
n=len(A)
A+=[2*A[n-1]-A[n-k-1]]
return
def alpha(k): #α_kをα^{k+1}-2α^k+1=0からニュートン法で求める
a=decimal.Decimal('2')
while (a-2)*a**k+1)/(((k+1)*a-k)*(a**(k-1))>10**(-d):
a=a-((a-2)*a**k+1)/(((k+1)*a-k)*(a**(k-1)))
return a
for k in range(m,K+1):
a=alpha(k)
L=[k]+[2**(i+1)-1 for i in range(k)] #L_nの初期値
count=0 #L_nとround(a^n)が連続して一致した回数
notfind=True
for i in range(N):
natti(L,k) #L_{k+i}を計算する
R=L[k+i]-round(a**(k+i))
if R==0:
count+=1
else:
count=0
if count==k: #L_n-round(a^n)がk回連続して0になったらn_kを出力して抜ける
print(f'n_{k}={i+1}')
notfind=False
break
if R>1000: #L_n-round(a^n)が大きくなりすぎても抜ける
break
if notfind: #n_kがNよりも大きいときやa^nの誤差が大きくなったときに抜ける
print(f'n_{k} not find')
break
ニュートン法の精度を上げるには有効数字decimal.getcontext().precを大きくすればよい。無闇に上げすぎても計算時間が増えるだけなので難しいところ(パソコンの性能次第ではありますが)。