この記事では
前回の記事
に続いて$k$-リュカ数の四捨五入表示について考察していきます。
さて前回までの考察で以下のことまでわかっていたのでした。
$k$-リュカ数$L_n^{[k]}$の四捨五入表示
$$L_n^{[k]}=\lfloor\a_k^n\rceil$$
が$n_0\leq n$で成り立つような$n_0$のうち最小のものを$n_k$とおくと
$$n_k<\frac{k^2(2k+1)}{\pi^2(1-\frac{\pi^2}{48})}\quad\l(<\frac{k^3}{3}\r)$$
が成り立つ。特に$n_k=O(k^3)$と評価できる。
前回の記事を書いてからというもの、どうにかここから評価を改善できないか試行錯誤したものですが何とも言えず仕舞いだったので、ここは一般の$k$について考えてるばかりでなく$n_k$の特殊値を見てみなければと思いさっそく$n_k$を計算してみることにしました。一旦計算の過程は省くとして$n_k$と$n'_k=n_k/k^3$の値について以下の結果を得ることができました。
$\begin{array}{ll}
\\ n_{2}=2&n'_{2}=0.2500
\\ n_{3}=4&n'_{3}=0.1481\ldots
\\ n_{4}=6&n'_{4}=0.0937\ldots
\\ n_{5}=11&n'_{5}=0.0880
\\ n_{6}=14&n'_{6}=0.0648\ldots
\\ n_{7}=23&n'_{7}=0.0670\ldots
\\ n_{8}=26&n'_{8}=0.0507\ldots
\\ n_{9}=38&n'_{9}=0.0521\ldots
\\ n_{10}=42&n'_{10}=0.0420
\\ n_{11}=57&n'_{11}=0.0428\ldots
\\ n_{12}=63&n'_{12}=0.0364\ldots
\\ n_{13}=80&n'_{13}=0.0364\ldots
\\ n_{14}=100&n'_{14}=0.0364\ldots
\\ n_{15}=121&n'_{15}=0.0358\ldots
\\ n_{16}=145&n'_{16}=0.0354\ldots
\\ n_{17}=170&n'_{17}=0.0346\ldots
\\ n_{18}=198&n'_{18}=0.0339\ldots
\\ n_{19}=228&n'_{19}=0.0332\ldots
\\ n_{20}=278&n'_{20}=0.0347\ldots
\\ n_{21}=312&n'_{21}=0.0336\ldots
\\ n_{22}=360&n'_{22}=0.0338\ldots
\\ n_{23}=409&n'_{23}=0.0336\ldots
\\ n_{24}=463&n'_{24}=0.0334\ldots
\\ n_{25}=519&n'_{25}=0.0332\ldots
\\ n_{26}=590&n'_{26}=0.0335\ldots
\\ n_{27}=665&n'_{27}=0.0337\ldots
\\ n_{28}=732&n'_{28}=0.0333\ldots
\\ n_{29}=815&n'_{29}=0.0334\ldots
\\ n_{30}=902&n'_{30}=0.0334\ldots
\\ n_{31}=993&n'_{31}=0.0333\ldots
\\ n_{32}=1088&n'_{32}=0.0332\ldots
\\ n_{33}=1202&n'_{33}=0.0334\ldots
\\ n_{34}=1323&n'_{34}=0.0336\ldots
\\ n_{35}=1431&n'_{35}=0.0333\ldots
\\ n_{36}=1544&n'_{36}=0.0330\ldots
\\ n_{37}=1695&n'_{37}=0.0334\ldots
\\ n_{38}=1817&n'_{38}=0.0331\ldots
\\ n_{39}=1979&n'_{39}=0.0333\ldots
\\ n_{40}=2148&n'_{40}=0.0335\ldots
\\ n_{41}=2323&n'_{41}=0.0337\ldots
\\ n_{42}=2504&n'_{42}=0.0337\ldots
\\ n_{43}=2650&n'_{43}=0.0333\ldots
\\ n_{44}=2842&n'_{44}=0.0333\ldots
\\ n_{45}=3083&n'_{45}=0.0338\ldots
\end{array}$
見るからに$n'_k$は$0.0333\ldots\fallingdotseq\frac{1}{30}$くらいの値に収束しそうな感じがします。つまるところ以下の予想が立てられます。
$\dis n_k\fallingdotseq\frac{k^3}{30}$が成り立つ。
この予想の精密化と解決が今後の課題となりますが、私には思い付きそうもないのでこの記事の読者にでも委ねようと思います。
具体的に$n_k$を求めるにあたって、手計算で挑むには$L_n^{[k]}\fallingdotseq\a_k^n\fallingdotseq2^n$であるので人力では分が悪すぎるということでこのためにpythonで簡単なプログラムを書いた(プログラミングにも数値計算にも明るいわけではないので賢いプログラムであるかは不明)。
手順としては$\a_k$の近似値をニュートン法で求め、$L_n^{[k]}$とそれに対する$L_n^{[k]}-\lfloor\a_k^n\rceil$を順に計算していき、$L_n^{[k]}-\lfloor\a_k^n\rceil$が$k$回連続して$0$になったら
ランダムウォークの手法
で$n_k$が決定する。$\a_k$の精度が荒くて$L_n^{[k]}-\lfloor\a_k^n\rceil$の値が大きくなったときは停止して精度を良くしてからやり直す。といった具合である。
実際に使ったプログラムは以下のようである。
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を大きくすればよい。無闇に上げすぎても計算時間が増えるだけなので難しいところ(パソコンの性能次第ではありますが)。