Gobel sequence の謎
ゲーベル数列は色々な部分によって異なる定義で示されていますね。
A003504 では
a(0)=a(1)=1; thereafter a(n+1) = (1/n)*Sum_{k=0..n} a(k)^2 (a(n) is not always integral!).
1,1,2,3,5,10,28,154,･･･
が並んでいき,a(43)まで整数でa(44)で有理数
より一般化(2乗をk乗にしてk-Gobel sequence)したもの
A108394 では
Least k for which f(k) = (1 + f(0)^n + f(1)^n + ... + f(k-1)^n)/k, f(0) = 1
n=2-->1,2,3,5,10,28,154,･･･
が並びf(42)までが整数でf(43)で有理数
n=3-->1,2,5,45,22815,2375152056927,･･･
が並びf(88)までが整数でf(89)で有理数
n=4-->1,2,9,2193,5782218987645,･･･
が並びf(96)までが整数でf(97)で有理数
n=5-->1,2,17,473297,5937570334133678310135701537,･･･
が並びf(213)までが整数でf(214)で有理数
これらは共通してf(1)=2となるタイプとして並ぶことになるので
さらのここをf(1)>=2 (即ちf(1)=3,4,5,6,･･･と変化させたことも含め)
A097398では
Matrix T(m,x(1)), m>=1, x(1)>=2, read by antidiagonals,
where T(m,x(1)) gives the position of the first noninteger term in the sequence
defined by x(n)=(x(n-1)*(x(n-1)^m+n-1))/n for n>=2
with exponent m and the given starting value x(1)
として
m\x1:--,2 ,3 ,4 ,5 ,6 ,7 ,8 ,9 ,10 ,11
1;--,43 ,7 ,17 ,34 ,17 ,17 ,51 ,17 ,7 ,34
2;--,89 ,89 ,89 ,89 ,31 ,151 ,79 ,89 ,79 ,601
3;--,97 ,17 ,23 ,97 ,149 ,13 ,13 ,83 ,23 ,13
4;--,214 ,43 ,139 ,107 ,269 ,107 ,214 ,139 ,251 ,107
5;--,19 ,83 ,13 ,19 ,13 ,37 ,13 ,37 ,347 ,19
6;--,239 ,191 ,359 ,419 ,127 ,127 ,239 ,191 ,239 ,461
7;--,37 ,7 ,23 ,37 ,23 ,37 ,17 ,23 ,7 ,37
8;--,79 ,127 ,158 ,79 ,103 ,103 ,163 ,103 ,163 ,79
9;--,83 ,31 ,41 ,83 ,71 ,83 ,71 ,23 ,41 ,31
10;--,239 ,389 ,169 ,137 ,239 ,239 ,239 ,239 ,239 ,389
の表が作られている。
確認のため
gp > gobel_kl(k, l, N) =
{
my(v = vector(N));
v[1] = l;
for(n = 1, N-1,
v[n+1] = v[n] * (n + v[n]^(k-1)) / (n+1);
);
return(v);
}
とプログラムを組み
上記表で最も少ない数で(m,x1)=(1,3),(1,10)での7を確認してみると
gp > gobel_kl(2, 3, 7)
%21 = [3, 6, 16, 76, 1216, 247456, 61235956672/7]
gp > gobel_kl(2, 10, 7)
%22 = [10, 55, 1045, 273790, 14992411852, 37462068869158688194, 1403406603957588515490448060329867110800/7]
となり確かに第7項目が初めて有理数となってしまう。
これ位の大きさなら十分通常のコンピュータの計算で確認は取れるが
上記の
n=5-->1,2,17,473297,5937570334133678310135701537,･･･
が並びf(213)までが整数でf(214)で有理数
などの様子など望むべくもありません。
一体どの様にして調べたのか教えてほしい。
f(213)を214で割った余りがわかればよい
→f(212)を213*214で割った余りがわかればよい
→f(211)を212*213*214で割った余りがわかればよい
→・・・
ということですので、例えば
gobel_kl(k,l,N)=
{
x=l;
for(n=1,N-1,
x=x%(N!/(n-1)!);
x=x*(n+x^(k-1))/(n+1);
);
return(x);
}
for(i=3,300,print(i," ",gobel_kl(5,2,i)))
のようにするとf(214)が非整数であることがわかりますね。
# 理由がよくわからないのですが、非整数が出てくると
# 次のx=x%(N!/(n-1)!)でエラーになりますので、
# 最初の非整数まで表示してエラー終了します。
# （非整数が出現したら終了というロジックを省略できて便利ではあります）
このプログラムでどこで有理数となる場所が判明できるのですね。
そこでこのプログラムを使って
A288641に並んでいる整数を算出していたら
forprime(i=3,900,print(i," ",gobel_kl(5,2,i)))で
･･･････････････
193 132351082600315838701
197 3340985130460343217
199 841924060846945
211 11068272793562666699
*** at top-level: ...rime(i=3,900,print(i," ",gobel_kl(5,2,i)))
*** ^-----------------
*** in function gobel_kl: x=l;for(n=1,N-1,x=x%(N!/(n-1)!);x=x*(n+x^(k-1)
*** ^---------------------------
*** _%_: impossible inverse in Fl_inv: Mod(0, 107).
の返しが起こりA288641でのn=5での値251にひっからないのです。
n＝2~20まではn=5以外ではちゃんと正解の値を与えてくれます。
この部分の謎は何なのでしょうか？
原因の予想はつきますが、対策はどうすればよいかわかりません。
n=5のとき最初に非整数になるのが214であり、214が素数でないから
という理由だと思いますが、非整数になってしまうとmod演算が
できなくなって今の方式では不都合です。
n=2～20はn=5を除きすべてたまたま素数項で非整数になりますが、
少し先を計算するとn=22のときも素数でない94で非整数になってしまい
179が得られません。
非整数を5乗とか22乗とかするとあっという間に分子分母が
巨大になってしまうのでうまくないですね。
さて、どうしましょう？
(追記)
上のプログラムで求まる43, 89, 97, 214, 19,・・・という数列はA108394ですね。
これを使うとてっきりA288641も見つかると思ってやっていました。
