チラ裏バッファ: 野蛮な解法
野蛮な解法
野蛮な解法とは…（このゲームの遊び方）
x2 − 421y2 = 1 のような式を満たす正の整数 x, y を見つける。
ただの数学の問題じゃん？ いやいやゲームなので「縛り」があるの。x なり y なりに 1, 2, 3, … を順に入れてみる方法（ブルートフォース＝「野蛮な力」）でこれを解きなさい、という縛りが…。暗算や筆算でやるのではなく、コンピューターに検索させる。1, 2, 3, … を順に入れて試してみるだけなので、計算機を使えるなら、誰にでもできるゲーム。
だったら総当たりのループを書くだけじゃん？ それで解けるならそれもありだけど、ゲームなので制限時間がある。理想は一瞬、目標は1秒未満、秒単位（10秒とか）も可だが、何時間もかかるようでは駄目。マシンや環境に左右されるのも不公平なので、直接検索していいのは15桁（999兆）までとしよう。
具体的には、x2 − dy2 = 1 の形の d が 2, 3, 5, 6, 7, 8, 10, 11, 12, … とだんだん増えていくとして（平方数の場合は「解なし」なのでスキップ）、どこまで解けるか？ d = 60 まではたわいもない。d = 61 は有名な「最初の壁」。フェルマーは、他の数学者たちに対して「このケースを解いてみろ！」と挑戦状を突きつけたという（インドでは何世紀も前に、とっくに解かれていたのだが…。インドすげえ）。現代のコンピューターをもってすれば、フェルマー・チャレンジは少しの工夫で突破できる。次の壁 d = 109 は、もっと難易度が高い。
基本技として、例えば x が奇数か偶数か分かれば、1, 3, 5, … のように一つおきに試せばいいので、2倍の高速化になる。同様に「3で割って1余る数」とか「4で割って2余る数」とか対象を絞れれば、それぞれ3倍・4倍の高速化になり、このような「フィルター」を組み合わせることで、何万倍というような高速化が可能なこともある（ただし「フィルターを作って組み合わせる」という準備作業にも、時間とメモリーが必要であり、やり過ぎると逆効果）。
別のアプローチとして、x2 − dy2 = 1 の代わりに
x2 − dy2 = −1
のような似た形の方程式の解を検索して、そこから間接的に問題を解いてもいい。うまくやると劇的な高速化になる。d = 181 を突破するには、この技が必須（最小の解が「直接検索していい範囲」の15桁を超えるので）。これで d = 330 まで進めるようになる。d = 331 が次の壁だが…
「野蛮な解法」というメモの大部分は、このゲームのプレイ記録である。
「野蛮な解法」第1部 インド人もびっくり
2020-04-22 野蛮な解法 インド人もびっくり
佐藤郁郎先生の「ｘ^2－９２ｙ^2＝１（その１）」に「数学者とは
ｘ^2－９２ｙ^2＝１
を１年以内に解ける人のことである」という言葉が書かれていた。ただし7世紀の話らしい。自明解 x=1, y=0 は除外しよう。
この手の不定方程式の整数解は、極端に難しいことがある。良い例は
x^3 + y^3 + z^3 = 33
だろう（2019年になってようやく解が得られた）。
しかし解が小さければ、全数検索という野蛮な方法が使える。さっそく PARI を起動して
test(x,y) = { if(x^2-92*y^2==1,print("x=",x,", y=",y)); }
brute(x) = { for(y=1,10000,test(x,y)); }
とタイプして
for(x=1,10000,brute(x))
とすると:
x=1151, y=120
1万以下の範囲には、これしかない。タイプするのはせいぜい1分、実行は10秒だが…もちろん、こんな解法、数学的には美しくない。21世紀のわれわれには7世紀の人にはない「おもちゃ」があるので、ゲームバランスが変わり、数学者の定義も変わってしまう。
やみくもに全数検索せず、少しは数学的に考えると、x が奇数、y が偶数なので（理由は後日説明）:
brute2(x) = { forstep(y=2,10000,2,test(x,y)); }
forstep(x=1,10000,2,brute2(x))
とすれば、とりあえず4倍の高速化が得られる。同様に mod 3, mod 5 等々で考えると、10倍くらい高速化するかもしれない。x と y は大ざっぱに同じ大きさ（例えば x が10くらい、y が1000くらいのとき、x^2−92y^2 が 1 になるわけない）という事実も、捜索範囲を絞るのに役立ち得る…
2020-04-23 野蛮な解法（その2） 素朴に検索
ｘ^2－９２ｙ^2＝１
の整数解を知りたいということは、要するに
９２ｙ^2+1
が平方数になるような自然数 y を見つけたい、ということ（言い換えれば、曲線が通る整数格子点）。ある数 n が平方数かどうか、どうやって判定しましょう。とりあえず「平方根の整数部分を2乗して元に戻るかチェック」というベタな実装を試してみる:
is_square(n) = { r=floor(sqrt(n)); r*r==n; }
（注: 数学的には r=sqrt(n); r==floor(r) の方がシンプルだが、計算機上では、その方法は必ずしも良いアイデアではない。）
その上で
forstep(y=0,10^9,2,if(is_square(y*y*92+1),print(y)))
のような検索を行うと、数分かかるものの
P0 = [1, 0]
P1 = [1151, 120]
P2 = [2649601, 276240]
P3 = [6099380351, 635904360]
が見つかり、x は末尾が1であること、y は120の倍数で1000を法として120ずつ増加することが予想される。さらに、P2以降について、x座標、y座標の値が、前の点の対応する値と比べて約2302倍ずつ増加すると予想される。こうした観察に基づき検索範囲を絞ると、与式の整数解たちを高速にどんどん生成できる:
P4 = [14040770918401, 1463851560480]
P5 = [32321848554778751, 3369785656320600]
P6 = [74404881332329766401, 7757245116998460720]
P7 = [171280004505174567476351, 17857174889544800256840]
...
x は末尾が1というだけでなく 01 と 51 が交互に現れている。比 x/y は92の平方根に近いが、
９２ｙ^2+1
の +1 の部分があるので、92の平方根ぴったりにはならない。
どうしてそうなるのか論理的に考えず、数値的に答えを見つける…という一種「野蛮」なアプローチ。エレガントの正反対のじゅうたん爆撃、物量攻撃。上記の is_square には「平方数かどうか知りたいだけなのに、全例について、いちいち平方根を精密計算している」という無駄（速度面の重大な欠陥）もある…。とはいえ、それなりに面白い結果が得られた。92の平方根は、23の平方根の2倍だが、この「23」という数と座標値の倍率「約2302」には関連があるのだろうか。
x, y 両方について片っ端から全数検索をかける方法は、非常に効率が悪い（P1 は容易に見つかるものの、P2 の先が見えてこない）。
2020-04-24 野蛮な解法（その3） あえて x を検索
ｘ^2－９２ｙ^2＝１
を満たす最初の3組の自然数解 P1, P2, P3 をブルート・フォースで得る方法として、前回は
９２ｙ^2+1
が平方数になるような偶数 y を検索した。同じ野蛮な方法でも、y ではなく、逆に
(ｘ^2－１)/92
が平方数になるような奇数 x を検索すると、多少、効率が上がる（約20%の高速化）。x は y より √92倍くらい大きいので、そのままでは検索範囲が10倍近くに広がって、ざっと10倍遅くなるし、割り算の手間も増える。けれど、別の角度から見ると、
ｘ^2－１
は92の倍数。簡単な考察によると x ≡ ±1, ±45 (mod 92) が必要条件。そのため、検索範囲内で実際にチェックする数は 4/92 つまり約4%で十分。検索範囲が10倍になっても、実際にはそのうちの4%だけを見ればいいのだから、実質40%の試行回数。前回のやり方は、範囲内の全偶数をチェックするのだから50%であり、結局、今回の方が2割お得。実装テストは:
T(x) = { n=(x*x-1)/92; y=floor(sqrt(n)); if(y*y==n, print([x,y])); }
を使い、次のように、歩幅92の大きなステップの中に、赤ちゃんステップを4個ずつ入れればいい。
forstep( a=0, floor(sqrt(92)*1e9), 92, T(a+1); T(a+45); T(a+47); T(a+91) )
これで「その2」と同じ範囲の総当たり検索になる。定量的に言えば、13ループ50試行で最初の解 P1 が得られる（前回の方法では60試行）。
2020-04-25 野蛮な解法（その4） 偶奇の判定・高速化
整数 x, y が
x^2 − 92y^2 = 1
を満たすとき、x は奇数、y は偶数。…今まで理由を記さず、この事実を使っていた。
x が奇数なのは一目瞭然だろう（もし偶数なら、左辺は偶数になってしまう）。従って、x は 8k±1 または 8k±3 の形を持ち、どちらにしても x^2 ≡ 1 (mod 8)。同じ理由から、もし仮に y が奇数なら y^2 ≡ 1 (mod 8) だが、そのとき与式の左辺は ≡ 1−92*1 ≡ 5 (mod 8) となり、題意を満たさない（左辺が整数の 1 に等しくなるためには、左辺を8で割ると1余ること、つまり左辺 ≡ 1 (mod 8) が必要）。ゆえに y は奇数になり得ない。
一方、y ≡ 0, 2, 4, 6 (mod 8) のいずれでも、92y^2 ≡ 0 (mod 8) であり、左辺 ≡ 1 (mod 8) となり得る。
この考え方は、いろいろな方向に拡張可能。例えば mod 16 では x ≡ ±1, ±7 が必要。「その3」からも分かるように x ≡ ±1 (mod 23) も必要。この二つを組み合わせると
x ≡ ±1, ±47, ±137, ±183 (mod 368)
となり、同一区間内のテスト対象の個数が、前回の半分に。従って、同じ検索が半分の時間でできる:
forstep( a=0, floor(sqrt(92)*1e9), 368,\
T(a+1); T(a+47); T(a+137); T(a+183); T(a+185); T(a+231); T(a+321); T(a+367) )
頑張ってもっと大きな法で考えれば、さらに10～100倍くらいは高速化できそうな感触だが…
2020-04-26 野蛮な解法（その5） もっと高速化
「その4」では、16*23 を法とすることで「その3」の検索を2倍、高速化した。「整数368個につき、8個の候補だけを調べればいい」というものだった。この手法を推し進めてみる。試しに mod 64 で考えると、x ≡ ±1, ±7, ±25, ±31 が必要条件。mod 9 で考えると x ≡ 0, ±1, ±3 が必要。これらを mod 23 と組み合わせると「整数13248個につき、80個の候補」。候補が3.6分の1に減り、3.6倍の高速化が見込まれる。
「その4」同様、80個のテストを T(a+1); T(a+231); T(a+321); T(a+415);... とベタで書いてもいいのだが、法をさらに大きくすることに備え
v = [ 1, 231, 321, 415, 505, 735, 737, 921, 1151, 1241, （中略） ];
という80要素のベクトル（配列）を別に作っておいて
tests(a) = { for( i=1, 80, T(a+v[i]) ) }
forstep( a=0, 10^10, 13248, tests(a) )
とする方が扱いやすい。ベタで書くより、多少速度が低下するかもしれないが、それでも「その4」の3倍は速い。（自明解 x=1 を別にすれば）たった8個の候補を調べるだけで、題意を満たす解 x=1151 が得られ、さらに第2解・第3解も容易に得られる。
ほんの少し工夫すれば、配列の長さを前半の40個だけにして、次のようにすることもできる。
tests1(a) = { for( i=1, 40, T(a+v[i]) ) }
tests2(a) = { for( i=1, 40, T(a+v[i]); T(a-v[i]) ) }
tests1(0); forstep( a=13248, 10^10, 13248, tests2(a) )
2020-04-27 野蛮な解法（その6） x^2−61y^2=1 の攻略
ペル方程式 x^2 − d*y^2 = 1 の正の整数解を 10^10 のようなオーダーまで全数検索することは、d=92 の場合には必要ない（2番目の解 P2 以降は、最小解 P1 から機械的に生成可能）。しかし
x^2 − 61*y^2 = 1 … ☆
とか
x^2 − 109*y^2 = 1 … ★
のような場合、一番小さい解 P1 が既に巨大。全数検索でそれを発見したいなら、検索の高速化が望まれる。
☆ の解は、条件 x ≡ ±1 (mod 61) を満たさなければならない。これだけでも候補をざっと30分の1に絞れるが、例えば x ≡ ±1, ±9, ±23, ±31 (mod 64) かつ x ≡ ±1, ±8 (mod 27) かつ x ≡ 0, ±1, ±5, ±10 (mod 25) という条件を併用すると、64*27*25*61個の整数につき、候補は448個。これは約6000分の1の絞り方。素朴に全数検索するより数千倍、速い。前回同様、
vec61 = [ 1, 3295, 13175, 18665, ..., 2635199 ];
のようなものを用意して、
T61(x) = { n=(x*x-1)/61; y=floor(sqrt(n)); if(y*y==n, print([x,y])); }
test61(a) = { for( i=1, 448, T61(a+vec61[i]) ) }
forstep( a=0, 10^10, 64*27*25*61, test61(a) )
とすると、1秒ほどで解
[1766319049, 226153980]
が得られる。これは速い！
…ように思えるが、10^10 の検索に1秒かかるなら、10^15 まで行くのに24時間以上かかる。★ の最小解は、そのくらい大きい。1万倍程度の高速化では、まだ ★ を素早く解くことができない。
☆ については、単に x ≡ ±1 (mod 61) という条件だけでも…つまり、448要素のベクトルを準備する手間を省いて、
forstep( a=0, 10^10, 61, T61(a+1); T61(a+60) )
とタイプするだけでも…割とすぐ解が得られる。素朴に考えると、y座標（比較的小さい）を検索対象にした方が、x座標を検索対象にするより効率的なはず。けれど「その3」で見たように（上記 ☆ の例からも分かるように）、x に関する条件 x^2 ≡ 1 (mod d) を利用するメリットは、検索範囲が √d 倍に広がるデメリットを補って余りある。
高速化が奏功するためには、一定の前提条件がある。「フェルマーのクリスマス～§5」の例8（x^2 − 79*y^2 = 5）のように「解なし」の場合、やみくもに広範囲をチェックしても、らちが明かない。
2020-04-28 野蛮な解法（その7） 平方数判定の高速化
ここまでの方法では、平方数かどうかの判定がベタで、非効率だった（「その2」参照）。「平方数かどうか」だけを知りたいのに、いちいち浮動小数点数で実際に平方根を計算するという…。
PARI の場合、組み込み関数として「整数の平方根の整数部分を求める」 sqrtint や「平方数かどうかを判定する」 issquare が用意されている。例えば
T(x) = { n=(x*x-1)/92; y=floor(sqrt(n)); if(y*y==n, print([x,y])); }
の場合、次のようにすると、機能は同じで、順に効率が上がる。
T(x) = { n=(x*x-1)/92; y=sqrtint(n); if(y*y==n, print([x,y])); }
T(x) = { n=(x*x-1)/92; if( issquare(n,&y), print([x,y]) ); }
T(x) = { n=(x*x-1)/92; if( issquare(n), print([ x, sqrtint(n) ]) ); }
最後の方法では、最初の方法に比べ時間コストが約半分で済む。PARI の組み込み関数には、上記の ☆ や ★ が一瞬にして解けるようなものもあるのだが…それは、ここで考えている全数検索的アプローチとは全く別の話。
2020-04-29 野蛮な解法（その8） 作戦会議
次のようなベタな平方数判定
T(x) = { n=(x*x-1)/92; y=floor(sqrt(n)); if(y*y==n, print([x,y])); }
に比べ、
T(x) = { n=(x*x-1)/92; if( issquare(n), print([ x, sqrtint(n) ]) ); }
とすると時間コストが半分で済む。これを
T(x) = { if( issquare((x*x-1)/92), print([ x, sqrtint((x*x-1)/92) ]) ); }
と書き直すと、そこからさらに 20～25% の時間短縮効果が得られた。PARI は単に「一時変数を作って代入する」というだけで、結構コストがかかるらしい。（プログラマーなら「n がグルーバルなのが最適化の妨げでは」と疑うだろう。ところが n をローカル宣言すると、むしろわずかに遅くなった。）
今までのやり方でも、
x^2 − 109*y^2 = 1 … ★
をブルート・フォースで10時間以内に解けそうだが、10時間も待ちたくない。候補絞りのベクトルを桁違いに巨大にして速くできるか、試してみたい（やり過ぎると空間コストの壁でかえって遅くなるだろうが、まだ行けそう）。
野蛮な方法（ブルートフォース）では当然そのうち限界がくるので、気が向いたら、教科書的な普通の方法（連分数）や、もうちょっと現代的な方法（二次体）も実装してみたい。ペル方程式は古典的な話題のようだが、「最小の正の解は必ずこの範囲にある」と上から押さえる問題は、現在でも未解決らしい。数学のロマンが絡んでる！
野蛮な方法であっても「小さい y より大きい x 側を検索した方が速い」というのは素朴な直観に反していて、ある意味、未知の領域かもしれない。「工夫すれば、ほとんどの場合、教科書の解法より幼稚なブルートフォースの方が速い」と結論できれば、それは一つの知見だろう。
2020-04-30 野蛮な解法（その9） x^2−109y^2=1 が解けた！
ペル方程式
x^2 − 109*y^2 = 1 … ★
の最小の正の解
x =158070671986249, y = 15140424455100
を次のようにして、ブルートフォースで得た。所要時間は、準備10分、検索10分程度だった。
5種類の法 512、81、25、49、109 について、★ の左辺が ≡1 になり得る x を調べた（0以上、法未満）。その結果を基に、解 x を M = 512*81*25*49*109 で割ったときの余りを考え、余りとなる可能性のある値の配列 vec109 を生成。vec109 の要素は50176個だが、M は約55億なので、素朴な全数検索に比べ11万倍の高速化が見込まれる。
この準備の後、次の3行のテストコードを実行し、上記の解を得た。
T109(x) = { if( issquare((x*x-1)/109), print([ x, sqrtint((x*x-1)/109) ]) ); }
test109(a) = { for( i=1, 50176, T109(a+vec109[i]) ) }
forstep( a=0, 2*10^14, M, test109(a) )
これは「その4」～「その6」と同様の高速化を大規模化したもの。
x^2 − d*y^2 = 1 の d が小さい場合、全数検索は実用上、非常に効果的。例外として、d = 61, d = 109 などについては、素朴な全数検索では時間がかかり過ぎる。前者については、適当に高速化して対処できたが（「その6」参照）、後者については、時間コストだけでなく空間コストが問題となり始めた。実際、M をあと2桁くらい大きくして ★ の解を数分で検索することにも成功したが、その場合、vec109 が巨大になり気軽にロードできないし、何より vec109 を準備する部分で時間がかかるので、トータルでは速くならなかった。
2020-05-03 野蛮な解法（その9の補足） 「その9」では M = 512*81*25*49*109 という数が天下り的に使われていますが、M の選び方に決まりはなく、でたらめに設定しても（所要時間は変わるとしても）解は得られます。大ざっぱに言って、単純に検索だけを考えた場合、「M の大きさと vec109 の要素数の比」（高速化率）が大きければ大きいほど効率が良く、「その9」では高速化率が10万倍を超えるように M が選択されています。幾つかの小さい素数べき（2^9 や 3^4 など）と素数 d のそれぞれを法として x の取り得る範囲を絞り、中国剰余定理のような感じで「ふるいを合成」（何重にもフィルターをかける）したわけです。
「野蛮な解法」第2部 勝ち方にスタイルは必要ないぜ！
2020-05-10 野蛮な解法（その10） いろいろな d についての x^2−d*y^2=1
d = 1, 2, 3, ... などと変わるとき、x^2−d*y^2=1 の最小の正の整数解 [x, y] は、どういう挙動を見せるでしょうか。次の（ばか正直な）コードを使って、PARIで検索してみた。
pell_test(d,y) = {
my( x2 = d*y^2+1 );
if( !issquare(x2), return(0) );
print( "d=", d, " : ", [ sqrtint(x2), y ] );
1;
}
brute_pell( d, yMax ) = {
if( !yMax, yMax=10^6 );
for( y=1, yMax, if( pell_test(d,y), return(1) ) );
print( "*** Failed with d=", d, " : yMax=", yMax );
0;
}
for(d=1,50,brute_pell(d)) とタイプすると:
*** Failed with d=1 : yMax=1000000 d=2 : [3, 2] d=3 : [2, 1] *** Failed with d=4 : yMax=1000000 d=5 : [9, 4] d=6 : [5, 2] d=7 : [8, 3] d=8 : [3, 1] *** Failed with d=9 : yMax=1000000 d=10 : [19, 6] d=11 : [10, 3] d=12 : [7, 2] d=13 : [649, 180] d=14 : [15, 4] d=15 : [4, 1] *** Failed with d=16 : yMax=1000000 d=17 : [33, 8] d=18 : [17, 4] d=19 : [170, 39] d=20 : [9, 2] d=21 : [55, 12] d=22 : [197, 42] d=23 : [24, 5] d=24 : [5, 1] *** Failed with d=25 : yMax=1000000 d=26 : [51, 10] d=27 : [26, 5] d=28 : [127, 24] d=29 : [9801, 1820] d=30 : [11, 2] d=31 : [1520, 273] d=32 : [17, 3] d=33 : [23, 4] d=34 : [35, 6] d=35 : [6, 1] *** Failed with d=36 : yMax=1000000 d=37 : [73, 12] d=38 : [37, 6] d=39 : [25, 4] d=40 : [19, 3] d=41 : [2049, 320] d=42 : [13, 2] d=43 : [3482, 531] d=44 : [199, 30] d=45 : [161, 24] d=46 : [24335, 3588] d=47 : [48, 7] d=48 : [7, 1] *** Failed with d=49 : yMax=1000000 d=50 : [99, 14]
d=1, 4, 9, 16, 25, 49 つまり平方数だと失敗するのが分かる。これらの場合、自然数 n があって d=n^2 (n≠0) と書けるので、与式は
x^2 − (n^2)*y^2 = 1
となる。これを成り立たせるような正の整数 y があったとして、m = n*y と置くと、
x^2 − (n*y)^2 = (x + n*y)(x − n*y) = 1
となるが、これは (x + n*y) = (x − n*y) = 1 または (x + n*y) = (x − n*y) = −1 を意味し、どちらの場合も y=0 となる。すなわち d が平方数のとき、正の整数解は存在しない（以下では考察対象外とする）。
最初にデフォルト検索（y が10^6以下の範囲を探す）が失敗するのは d=61 のケース。この場合の解は「その6」で既に得ている（y は 10^8 のオーダー）。次に d=97 でも失敗するが、これは捜索範囲を10倍に広げて
brute_pell(97,10^7)
とすれば、それだけで解決。d = 101...150 の範囲では 106, 109, 139, 149 が難しい。106 と 139 はちょっと捜索範囲を広げれば済むが、d=109 はしんどい。「その9」で20分くらいかけて無理やり解いたのだった…。d=149 の場合の解はそこまで大きくはないが、d=61 の場合よりは大きな数になる。
ばか正直に全数検索しないで解だけ見たい場合、PARI では次のチートコードを使える。
cheat_pell(d) = { my(z=quadunit(4*d)); if( norm(z)==1, z, z^2 ); }
2020-05-11 野蛮な解法（その11） どんな d に対して x^2−d*y^2=1 は難しい？
一番小さい例として d=2 を考えてみましょう（前回見たように d は平方数では駄目なので d=1=1^2 は考察対象外）。
x^2−2y^2=1 の最小の正の整数解は、次のようにして、簡単に見つかります。まず y=1 だと x^2−2=1, x^2=3 ですが、これは無理。次に y=2 だと x^2−8=1, x^2=9 ですが、これはもちろん x=3 という解を持ちます。
同様に d=3 や d=5 の場合の解…つまり x^2−3y^2=1 や x^2−5y^2=1 の解…は、y=1, 2, 3, ... などを順に試せば、暗算でもすぐ見つかるでしょう。前回は、この「順に試す」という部分をコンピューターにやらせたのでした。「野蛮な解法」と呼ぶわけは、こんな「片っ端から調べる」ことをしなくても、もっとスマートな解法があるから。でも d が小さい場合、「エレガントな解法」など考えるより、直接全数検索した方が手っ取り早いのも事実。勝ち方に美学は要らないぜ！
とはいえ d=13 すなわち x^2−13y^2=1 になると、暗算ではきつい。d が大きいと必ず難しい…というわけではなく、x^2−20y^2=1 のように、掛け算の九九の範囲で、簡単に暗算できる場合もあります。一体 d の値と、最小の正の解の難しさ（どのくらいでかくなるか）には、どういう関連性があるのでしょうか。前回の PARI コード検索によると、d が大きくなるにつれ「これまでなかった新記録のでかい解」が出るケースとして、次の例があります。
d=13 のとき y=180; d=29 のとき y=1820; d=46 のとき y=3588; d=53 のとき y=9100; d=61 のとき y=226153980; d=109 のとき y=15140424455100; ...
このうち d=46 を除くと、d が「4で割って1余る素数」のとき、新記録が出やすいようです（さらによく見ると「8で割って5余る素数」がヤバそう）。素数に限らず d が4で割って1余るとき、題意を満たす y は偶数になるようです。 d が4で割って3余るときには奇数になりやすいですが、これは例外あり。「4で割って3余る素数」に限れば、y は奇数になるようです。さらに、素数に限らず、d が4で割って2余るとき、y は偶数になるようです。d=92 のケースについて、与式を満たす y は必ず偶数になることを既に証明しましたが（「その4」参照）、このような証明は一般化できるでしょうか…？
2020-05-14 野蛮な解法（その12） 真面目に d=4k+1型
x2 − 5y2 = 1 の解は x = 9, y = 4 など。
x2 − 13y2 = 1 の解は x = 649, y = 180 など。
x2 − 17y2 = 1 の解は x = 33, y = 8 など。
一般に d ≡ 1 (mod 4) つまり d が「4で割って1余る整数」のとき、x2 − dy2 = 1 の整数解は、x が奇数、y が偶数になると予想されます。それを示すため、この条件で y が奇数になり得たと仮定してみましょう。つまり y を4で割ると1または3余ると仮定。これは y ≡ ±1 (mod 4) と同じ意味なので、この仮定上の世界においては y2 ≡ 1 (mod 4) であり、d は奇数なので x は偶数、その平方は4の倍数…すなわち x2 ≡ 0 (mod 4)。ということは
x2 − dy2 ≡ 0 − 1×1 = −1 ≡ 3 (mod 4)
ですが、これでは条件に合いません（方程式の右辺が1なので、x2 − dy2 = 1 ≡ 1 (mod 4) となる必要がある）。ゆえに y が奇数になり得るという仮定は誤りであり、y は偶数でなければ駄目、ということが分かります。同様に mod 8 で考えると、d が4で割って1余る場合、y は（単に偶数というだけなく）4の倍数になる…ということを確認できます。よって、この場合、全数検索するにしても、4の倍数だけを試せば十分。「ベタの全数検索」と「4の倍数限定の全数検索」を比較すると…
brute1( d ) = {
if( issquare(d), return );
for( y=1, 10^7, if( pell_test(d,y), return ) );
print( "*** Failed with d=", d );
}
brute4( d ) = {
if( issquare(d), return );
forstep( y=4, 10^7, 4, if( pell_test(d,y), return ) );
print( "*** Failed with d=", d );
}
forstep( d=1, 120, 4, brute1(d) ) と比べ forstep( d=1, 120, 4, brute4(d) ) は、約4倍速くて結果は同じ。
2020-05-15 野蛮な解法（その13） 奇数の平方を4で割ると1余る（8で割っても1余る）
ペル方程式 x2 − dy2 = 1 において、今度は d が 4k + 2 型の自然数（4で割ると2余る数）である場合を考える。
x2 − 2y2 = 1 の解は x = 3, y = 2 など。
x2 − 6y2 = 1 の解は x = 5, y = 2 など。
x2 − 10y2 = 1 の解は x = 19, y = 6 など。
x が奇数で y が偶数になりそうだ…。ここでは d が偶数（ただし4の倍数以外）と仮定しているので、x が奇数になるのは明白だろう。さもなければ、方程式の左辺は偶数ひく偶数で偶数になってしまい、右辺の1と等しくなり得ないから。ゆえに x ≡ ±1 (mod 4), x2 ≡ 1 (mod 4) であり、このタイプのペル方程式が成り立つ条件は 1 − dy2 ≡ 1 (mod 4) すなわち dy2 ≡ 0 (mod 4)。これは y が偶数であることと同値（d は4の倍数ではないので）。「その4」では、d = 92 のケースについて、同様のことを mod 8 で考えた。別の狙いがあって、「その4」以降では mod 16 や mod 64 などでも考えたが、偶奇の分類だけなら mod 4 で考えれば足りるし、d ≡ 2 (mod 4) 型を考えているのだから、それが自然でもある。
「野蛮な解法」第3部 小鳥のホテルは不思議がいっぱい
2020-05-28 野蛮な解法（その14） ふしぎなふしぎなハト・ホテル
次の短いPARIコードの意味は…
pigeon(x,N) = {
my( h = Map(), f, r, j );
for( k=0, N, f=frac(k*x); r=floor(f/(1/N));
if( mapisdefined(h,r,&j), return(k-j), mapput(h,r,k) ) );
}
N個の客室があるハトのホテル（ハト小屋）を作ります（Nは自然数）。k=0番のハトから、k=N番のハトまで、最大 N+1羽のハトさんをチェックインさせようとしますが、客室はN個しかないので、絶対どこかで問題（ダブルブッキング）が生じます（笑）。k番目のハトには、k*x の小数部分 f に応じて r号室を割り当てますが、ダブルブッキングの問題があるので、ハトを入れる前に、そこが今空室かどうか確認。既に先客がその部屋を使っている場合、割り当てマップを見ると「j番のハトがその部屋を使用中」と分かるので、その場合、ホテルを終了し k−j を返します。そういう問題がなければ、ハトさんを r号室（現在空室）に入れ、ホテルの営業を続けます。例えば…
gp > pigeon(Pi,10) k=0 : k*x=0 : r=0 k=1 : k*x=3.141592653589793238462643383 : r=1 k=2 : k*x=6.283185307179586476925286767 : r=2 k=3 : k*x=9.424777960769379715387930150 : r=4 k=4 : k*x=12.56637061435917295385057353 : r=5 k=5 : k*x=15.70796326794896619231321692 : r=7 k=6 : k*x=18.84955592153875943077586030 : r=8 k=7 : k*x=21.99114857512855266923850368 : r=9 k=8 : k*x=25.13274122871834590770114707 : r=1
この例では、円周率について、部屋が10個のハト・ホテルを作りました。円周率の2倍は小数部分が0.2台なので、ハト2を2号室へ入れる。円周率の3倍は小数部分が0.4台なので、ハト3を4号室に。…同様に続けると、ハト8は1号室に行くはずですが、そこは既にハト1が使用中なので、そこでストップ。8−1=7 を返します。これは「円周率の7倍は整数に近い」ということを意味しています（なぜでしょう？）。
もっと部屋の数が大きい（細かく部屋割りをできる）ホテルの例。pigeon(Pi,1000) とすると 113 が返ります。pigeon(Pi,10^6) なら 66317 が…。実際、
113π = 354.9999698556466359462787023
66317π = 208341.0000081143181951271213
は、非常に整数に近い値を持っています。これらの数値に見覚えのある方も多いでしょう。特に π ≈ 22/7 や π ≈ 355/113 は有名なので…。
円周率に限らず、任意の正の無理数 x について、同様のことを実行でき、この原理を出発点として「d が平方数の場合を除き、ペル方程式には、必ず正の整数解が存在する」ということを証明できるのです。
注: このコードは必ずしも最適近似を与えません。「一番際どいニアミスが起きる場所」を全数検索せず、単にダブルブッキングが判明したら、緩めのニアミスでも検索を打ち切るから。とはいえ、1/N オーダーのニアミスを誘発するので、最適でないとしても、N を大きくすればするほど、いくらでも高精度の近似分数が得られます。
2020-05-31 野蛮な解法（その15） 「小鳥のホテル」は不思議がいっぱい
pigeonhole（ピジンホール）を直訳すると「ハト穴」、数学用語としては「鳩の巣」と呼ばれますが、この用語が最初に使われたとき、もともと「ハトの宿舎のように、仕切りのついた区画がいっぱい並んでいる整理ボックス」のような意味だったとか…。でも、かわいいので、ここでは「小鳥のホテル」と呼んじゃいましょう！
「小鳥のホテル」の原理（鳩の巣原理）は、要するに「客室数より多い小鳥さんをチェックインさせたら、最低2羽は相部屋になるよ」…どう考えても当たり前ですが、この当たり前のことが、いろんな不思議と関わってくる…。例えば、
53√89 =
500.0009999990000019999950000139999580001319995…
は、前回の pigeon で遊んでいて見つけた不思議な数。無理数なのに数字の並び方が興味津々。計算自体は簡単に試せるので、Windows ユーザーの方は電卓を起動して [Alt]+[2] で科学計算モードにしてから、
53*89@=
とタイプしてみてください。
さて、この原理を使って x2 − dy2 = 1 には必ず正の整数解があることを証明するのが本来の目標ですが（それは、パズルとして面白いだけでなく、数学的にも深い意味を持つ）、ここでは寄り道して「円周率にこの計算法を適用したらどうなるか？」ということを少し丁寧に考えてみましょう。
わたしたちは、客室が10個だけある「小鳥のホテル」を作りました。次の11個の数を11羽の小鳥だと思うと…
0π, 1π, 2π, 3π, 4π, 5π, 6π, 7π, 8π, 9π, 10π
…部屋の数より鳥の方が多いので、少なくともどれか2羽が相部屋になることは明白。部屋の割り当て方は、小数部分が 0.0台なら「0号室」、0.1台なら「1号室」、0.2台なら「2号室」…0.9台なら「9号室」として、例えば、8π = 25.13… は「1号室」に入れます。でも 1π = 3.14… も「1号室」なので、1π, 8π は相部屋に。これは、jπ, kπ の小数部分の差が 1/N 未満になるような、0以上N以下の相異なる整数 j, k が存在する…ということ（この場合 N=10, j=1, k=8）。ホテルの構造上、客室数より客の方が多いのだから、問題にする数が π でなくても、あるいは、部屋の数 N が10でなくても、常に同様のことが成立。
かくして 1π, 8π の小数部分はだいたい同じなので（両方とも0.1台）、この二つの数の差は、整数に近い:
7π = 8π − 1π = 25.1327412… − 3.1415926… ≈ 21.9911486
= 22 − 0.0088514
この 0.0088514 の逆数 を計算すると、112.976…。つまり 7π は、整数22と、距離1/112 程度のニアミスを起こす。幅1/10の半開区間に2羽の小鳥を入れるので、小鳥同士が距離1/10未満まで接近するのは確実ですが、実際には1/112以内まで接近しています。この観察を基に、今度は小数部分を1/113単位で区切れば（※）、さらに際どい1/113未満のニアミスが確実に起きる（※小数部分が0/113以上1/113未満なら0号室、1/113以上2/113未満なら1号室…というふうに、113室のホテルに114羽の小鳥たちを割り当てる）。
やってみると…0番の小鳥 0π = 0.0 と、106番の小鳥 106π = 333.008821… の小数部分の差は、実際 1/113 = 0.008849… 未満（N=113）。106π は距離 0.008821 ≈ 1/113 まで整数に接近しました。前回が1/112程度なので、1/113では大して違いがないようですが、
22/7 = 3.142857…
333/106 = 3.141509…
なので、近似精度は確実にアップ！ 調子に乗って、今度は部屋の幅を1/114にすると（N=114）、1π = 3.14159… と 114π = 358.14156… の小数部分が、1/33173 程度のすごいニアミス。両者の差を考え分数を作ると:
355/113 = 3.14159292…
小数6桁目まで正しい（有名な）近似分数が得られました。さらに調子に乗って N=33174 にすると、ニアミス度は約1/52275。9桁目まで正しい近似分数が得られます:
103993/33102 = 3.14159265301…
以下同様に、ニアミス度の分母が n 台のとき、n+1 を次の N にして新しいホテルを作ればニアミス度が上がり、そこから、より優秀な近似分数が得られます。教科書的にはこれは「円周率の連分数展開」というお決まりの話題ですが、このように「小鳥のホテル」を考えると…
必ずニアミス（相部屋）が起きるのは当たり前（整数に近い数が現れるので、それを使って簡単に近似分数を作れる）
Nの選択によって、必ず前回より際どいニアミスを引き起こせる＝無限に（いくらでも高精度の）近似分数を作れる
…ということが実感でき、天下り的で無味乾燥な連分数の話より面白いのでは…。といっても、実用上の計算は、直接、連分数として計算した方がはるかに効率的ですが。ときめく方がいいよね！
ボーナスクイズ: プログラム上で「少なくとも10万羽の小鳥を収容できるホテル」を実装し、あと3歩先まで近似分数を作ってください。分母が10万までの範囲での、円周率の最良近似分数を見つけてください。実行する前に、その近似分数の値が円周率の真の値とどのくらい一致するか、勘で予想してください。小数20桁くらい？ 10桁くらい？ 10桁いかない？
2020-06-01 野蛮な解法（その16） サインカーブの波打ち際
「小鳥のホテル」で遊んでいたら、今度は、こんなもんが見つかりました:
1 + sin(344) = 0.00000966049382909903…
1 + cos(355) = 0.00000000045434101983…
前回の系: π ≈ 333/106, 333 ≈ 106π、ゆえに sin(333) ≈ 0, cos2(333) = 1 − sin2(333) ≈ 1。π ≈ 355/113, 355 ≈ 113π についても同様。
きちんと文献を調べたわけではありませんが、鳩の巣原理を使って 1/N 未満のニアミスの存在を導くアプローチは、1842年に Dirichlet が発見したらしい（？）。だとすると、Lagrange が、ペル方程式の解の存在を最初に証明したときは（1760年代？）、別のアプローチを使ったのかな…。
意味ということは、定義上、心の中で起きているマッピングであり、解釈がなければ、意味があるということもないし、意味がないということもないわけです。そうすると、実用上・便宜上の観点からは、起きることには全て意味があると解釈した方が得でしょう。例えば…
ホウレンソウを買うつもりだったのに、間違えて小松菜を買ってしまった
…というとき、「意味のないばかな行為」と解釈して悲しい気持ちになることもできるのですが、楽天的に「今日は小松菜を食べるとラッキーというお告げかもしれないな。小松菜で何を作ろうかな～」と考えてもいいわけで。もっとすごくひどいことが起きた場合でも、何という不条理だと悲嘆に暮れることもできるけれど、「何だか分からないけど、これにはきっと深い意味がある」と考えてみてもいいわけです。注意点として、解釈はあくまで自分の心の中にとどめておくのが無難。悲しんでいる第三者に対して安易に「これは意味のある出来事。元気を出して」などと言うと、往々にして相手が気分を害し、関係がこじれる原因になるからです。ホロコースト犠牲者の遺族に、第三者が「ホロコーストには意味があったんだよ」などと言わない方がいいのは明白でしょう。自分自身が遺族の場合にも「意味がある出来事でした」などというと「頭がおかしい」と思われかねないので、やはり「解釈」はプライベートな心の中の問題にとどめておくのが吉かと…。選択の余地がない場合…例えば、裁判のような公的な場で「遺族の意見」を問われた場合には、控えめに「なるべく軽い罰にしてあげて」くらいは言うかもしれませんが…自分に無関係な論争には、どの陣営にも加わらないのが得策だと思うのです。
しかし sin(344) の値を考えることに何の意味があるのか…と問われた場合、ぶっちゃけ意味がないといえばない。「連分数の循環の周期を数えれば済む話なのに、なぜ全数検索するのか」「なぜ教科書通りに効率的にやらないのか」と問われた場合、何と答えればいいのでしょう…。例えばの話、世の中には「チェスに興味がない人」もいれば「チェスが好きな人」もいれば「チェスなんて、どうせコンピューターにかなわないのに人間がやっても意味がないと考える人」もいるわけで、どの立場であるにせよ、それは心の中の解釈の問題。他人を説得して同意させる必要があることではないでしょう。
2020-06-12 野蛮な解法（その17） 53√89再訪
53√89 =
500.0009999990000019999950000139999580001319995…
は、数字の並び方が面白い。
(53√89)2 = 532 × 89 = 250001 = 5002 + 1
なので、53√89 は 500 よりわずかに大きな数。同様に、例えば、
13 × 13 × 147929 = 25000001 = 50002 + 1
なので、
13√147929 =
5000.0000999999990000000199999995000000139999995800000131999995710000142999995…
これだけでは、数字の並び方の面白さを詳細に説明できないが、とりあえず、だいたいこういう大きさの「整数に近い値」になることは分かる。そして
532 × 89 = 5002 + 1
は、負のペル方程式
x2 − 89y2 = −1
の基本解 (x, y) = (500, 53) と関係あることも分かる。
2020-06-13 野蛮な解法（その18） 筆算で√250001
53√89 = 500.0009999990 0000199999 5000013999 9580001319 9957100142 9995138016…
の数字パターンは、最初のうち「きれい」だけど、だんだん崩れて、結局はごちゃ混ぜになる。これを感覚的に把握するため、小数第21位の 5 まで直接、筆算してみる。
53√89 = √250001 の根号内は500の2乗プラス1なので、筆算で計算すると500が立って、かすの1が余る。下記の [A] 地点。結果は [B]。攻撃対象（割られる数＝分子）は、3桁レベル最弱の100。こちらの攻撃力（割る数＝分母）は1万台の5桁レベルなので、相手が弱過ぎて戦闘にならず、経験値が稼げない。戦果 0 が続く。
5 0 0. 0 0 0 9 9 9 9 9 9 00000 1 99999 5
----------
25 00 01 5
25
00 10_0
00
01 100_0 [A]
0
100 1000_0 [B] 3d vs. 5d
0
10000 10000_0
0
1000000 100000_0
0
100000000 1000000_ [C] 9d vs. 8d
1ターンごとに敵は2桁ずつ回復する一方、こちらの攻撃力は1桁ずつしかアップしないので、3ターン後、再び戦闘状態に突入。それが [C] 地点。敵は9桁レベル最弱だが、言い換えれば8桁レベル上限 99999999 プラス1。一方、当方のレベルは8桁レベル最軽量級の 1000000x なので、ここからはおいしい。最大戦果9を挙げ、さらにお釣りから9を挙げ…という状態がしばらく続く。具体的には、次のように、6連続「9」ゲット。小数第4位以下の部分に当たる。
100000000 1000000_9 [C]
90000081
999991900 10000018_9
900001701
9999019900 100000198_9
9000017901
99900199900 1000001998_9
90000179901
990001999900 10000019998_9
900001799901
9000019999900 100000199998_9
9000017999901
199999900 1000001999998_ [D] 9d vs. 14d
[D] 地点において、敵の勢力は、1から始まる9桁レベルにまで落ちてしまう。当方は14桁（「相手から14桁の数を1回引くごとに1ポイント」というレベルの戦闘中）なので、これ以上の戦果は挙げられない。敵が2桁ずつ回復する間、こちらは地味に1桁ずつ大きくなって、戦闘再開を待つ。
199999900 1000001999998_0 [D] 19999990000 10000019999980_0 1999999000000 100000199999800_0 199999900000000 1000001999998000_0 19999990000000000 10000019999980000_0 1999999000000000000 100000199999800000_1 [E] 19d vs. 19d 1000001999998000001 99999700000199999900 1000001999998000002_9 [F] 20d vs. 20d
[D] 地点から戦果「0」が5連続。小数点以下0が5個続く部分に相当する。[E] 地点で19桁同士にまで回復。「1」が立つだけだが、1倍を引くと、[F] の状態になり、敵はオール9に近い20桁。当方は1から始まる20桁なので、また荒稼ぎができそう…。しかし冷静に分析すると、前回ほどは引きまくれない。[C] においては、敵はオール9より一つ上という燃料満載状態、当方は最下位桁以外 100... という軽量状態だったので、相手から引きまくることができたが、[F] の敵はオール9よりは劣化しているし、当方も桁の真ん中あたりに「ノイズ」が生じている。敵の最上位桁には9が5個しかないので、9のダメージを最大5回与えたら、そこで打ち止め。具体的には…
99999700000199999900 1000001999998000002_9 [F]
90000179999820000261
999952000037999963900 10000019999980000038_9
900001799998200003501
9995020003979996039900 100000199999800000398_9
9000017999982000035901
99500200399799600399900 1000001999998000003998_9
90000179999820000359901
950002039997960003999900 10000019999980000039998_9
900001799998200003599901
5000023999976000039999900 100000199999800000399998_5 [G]
5000009999990000019999925
1399998600001999997500 1000001999998000003999990_ [H] 22d vs. 26d
5連続で「9」が立ったものの、敵の劣化・当方のノイズ（微妙な肥大化）の影響を受け [G] では「5」倍しか引けない。[H] では22桁 vs. 26桁となり、[D] に似たシチュになるものの、当方のノイズも増加してしまった。以下、敵の回復を待って最上位桁の「13」を消せたとしても、その先に9が4個しかないので、次の連続「9」ゲットは最大でも4回と予想される。小数31桁までは
53√89 = 500.0009999990 0000199999 5000013999 9…
となったとして、その先ではさらに劣化が進行し「ごちゃ混ぜ」度がアップするだろう。しょせん無理数なので、循環小数にはならず、結局「ランダム」になってしまう。
パターンの本質が解明されたわけではないが、パターンが「だんだん崩れる」仕組みは何となく分かった。
注: 平方根の筆算については、簡略表記。
2020-06-16 野蛮な解法（その19） 500+53√89 と黄金比
53√89 を小数で書くと面白いパターンが見られること…。それは「鳩の巣原理」の具体例をいろいろ考えていて、たまたま気付いたことだったが、
532 × 89 = 5002 + 1
であり、このことは
x2 − 89y2 = −1
の解 (x, y) = (500, 53) とも関係しているのだった。この場合の「鳩の巣原理の例」とは「a√89 と b√89 の小数部分がほとんど等しいような（※）、そんな相異なる自然数 a, b が存在する」ということだった。（※）の意味＝「小数部分の差を1万分の1以内にしろ」「1兆分の1以内にしろ」…等々、どんなに際どいニアミスを要求されたとしても、必ずその条件を満たすような a, b がある。
ところで x = 500 + 53√89 は、二次方程式
x2 − 1000x − 1 = 0
の解。これを変形すると:
x − 1000 − 1/x = 0
x = 1000 + 1/x ← この関係を(♪)と呼ぶことにする。
(♪) の右辺の分数の分母には x があるが、左辺を見ると、(♪) は x を定義する式…と解釈可能。そこで右辺の x に（♪）の右辺自身を代入すると:
x = 1000 + 1/(1000 + 1/x)
この手順は、いくらでも繰り返すことができる。左辺の x を定義する右辺に、定義したい x 自身が含まれているので、代入すると一種の無限ループ（「再帰的」と呼ばれる状態）が発生する:
x = 1000 + 1/(1000 + 1/(1000 + 1/x))
x = 1000 + 1/(1000 + 1/(1000 + 1/(1000 + 1/x)))
…
x = 500 + 53√89 と書けば済むのに、なぜわざわざこんな書き方をするのか。第一に、
500 + 53√89 = 1000.0009999990 0000199999 5000013999 9580001319 9957100142 9995138016…
…という一見不思議な数字の並び方が（少なくともある程度）説明される。連分数表示を見ると、小数表示がこんな感じになってもおかしくないような気がしてくる。第二に、役に立つとか分かりやすいとかそういうことと無関係に、ワクワクする。詳しい人から見れば自明のことだろうし、興味ない人にとっては「小難しい数学の話」なのだろうけど、偶然 53√89 の数字の並び方が目に留まったら、素朴に「不思議なパターンだなぁ」と思わないだろうか。第三に、連分数表示は理論的にも重要な意味を持つ。
上記は黄金比 x = (1 + √5)/2 ≈ 1.618 を定める二次方程式（参考記事）
x2 − x − 1 = 0
x = 1 + 1/x
と似ている。（♪）と同様の原理から、黄金比も
x = 1 + 1/(1 + 1/(1 + 1/(1 + 1/x)))
のような形で書くことができる。
「野蛮な解法」第4部 解の存在証明
2020-06-18 野蛮な解法（その20） ペル方程式は非自明な解を持つ（前編）
鳩の巣原理の話から脱線しまくってしまったが、本題に戻って、次のことを証明したい。「平方数ではない任意の自然数 d について、
x2 − dy2 = 1
は無限個の整数解を持つ」。これだけなら、どの教科書にも書いてあるようなことだが、そこで止まらず、最終的には「フェルマーのクリスマス定理で遊ばせて！」§5、例8で説明をうやむやにした次の事実を示すのが目標である。
事実:「x2 − 79y2 = 5 は整数解を持たない」
この事実は「剰余演算や平方剰余では証明できない」という意味において、初等的な方法では歯が立たないように思われるのだが、どうも高校数学のような範囲の「初等的な道」がありそうなのだ！ 簡単ではないが、何回かに分けてその道をたどってみたい。
補題1 無理数 √d をうまく整数倍すれば、結果はいくらでも整数に近くなる。具体的に、
| x − y√d | < 1/y （☆）
を満たす正の整数のペア (x, y) が無限個存在する。
〔証〕 任意の整数 N ≥ 2 を選ぶ（例えば N = 10）。鳩の巣原理により、j√d, k√d の小数部分の差の絶対値が 1/N 未満になるような、整数 j, k が存在する（「その15」参照）。ただし 0 ≤ j < k ≤ N。このとき k√d − j√d は整数に近いので（最寄りの整数との差が 1/2 未満）、その整数を x とする。k − j の最大値は N であり、その場合 y = N と置けば（☆）成立。そうでなければ、y = k − j と置くと、（☆）左辺は 1/N 未満だから 1/y 未満であり、やはり（☆）成立。どちらにしても（☆）左辺は 1/y 未満。これをもっと厳密に評価して「1/s 未満だが 1/(s + 1) 未満ではない」としよう（s は y 以上の整数）。今 s + 1 をあらためて N と置けば、上記と同じ議論から、N に対応する新しい (x, y) があって（☆）左辺はさらに小さくなる。
このとき、古い (x, y) では（☆）左辺は 1/(s + 1) 未満ではなかったが、新しい (x, y) では（☆）左辺が 1/(s + 1) 未満になるので、古い組と新しい組は同一ではない。
この手順を反復して「さらに小さく」「さらに小さく」…と（☆）左辺をいくらでもゼロに近づけることができる。反復のたびに新しいペア (x, y) が見つかり、それは既に見つかっているどのペアとも同一ではない。√d は無理数なので、有限回の反復では（☆）左辺はゼロにはならない。つまり（☆）を満たす (x, y) は無限個存在。□
〔補足〕 j < k なので y = k − j ≥ 1。y = N は、そのうち j = 0, k = N の場合に当たる。k√d − j√d = y√d は、最寄りの整数 x との差の絶対値が最悪（N = 2 の場合）でも 1/2 未満: y√d − 1/2 < x < y√d + 1/2。y は 1 以上なので y√d − 1/2 は正であり、従って x も 1 以上。√d は無理数なので、その整数倍がちょうど 1/2 の端数を持つことはない。0 倍すれば 0 にはなるが、（☆）左辺の y は 0 ではなく、（☆）左辺が 0 になる可能性はない。
2020-06-19 野蛮な解法（その21） ペル方程式は非自明な解を持つ（中編）
補題2 平方数ではない任意の正の整数 d に対して、
x2 − dy2 = c （△）
が無限個の正の整数解を持つような整数 c が存在する。
例1 d = 10 のとき、c = 6 は条件を満たす。
x2 − 10y2 = 6 の解は
x = 4, y = 1 ⇒ x2 − 10y2 = 6
x = 16, y = 5 ⇒ x2 − 10y2 = 256 − 10 × 25 = 6
x = 136, y = 43 ⇒ x2 − 10y2 = 18496 − 10 × 1849 = 6
など。（△）において c = 1 にできればペル方程式の解の存在証明が完了するのだが、まずは（△）まで。
補題2の証明 補題1（☆）を利用して、（△）左辺の絶対値をこう書く:
| x2 − y2d |
= (x + y√d) × | x − y√d | < (x + y√d) / y = x/y + √d （ク）
これは恒等式ではなく、あくまで「補題1を満たす正の整数 x, y については、これも成り立つ」という性質のもの。どうにかして一番右の分子の x を消せれば「（△）左辺は、絶対値において、定数 √d のオーダーを超えない」「だから可能な c の値は有限個」「だから無限個の解を考えると、必ず c の重複が起きる」…と論じられる。
不等式（ク）の x を消すため、（☆）を使って:
x = (x − y√d) + y√d ≤ | x − y√d | + y√d < 1/y + y√d
これを（ク）の一番右に代入すると:
| x2 − y2d |
< 1/y2 + √d + √d ≤ 1 + 2√d （クク）
y は正の整数なので（補題1参照）、1/y2 ≤ 1 である。
すなわち、補題1の条件を満たす x, y のペア（無限個存在）を（△）左辺に入れたとき、（△）の式の値は、絶対値において 1 + 2√d を超えることがない。これら無限種類の入力に対し、出力は有限の範囲の整数値なので、出力値のうち少なくとも1種類は「無限個の入力に対して、この値が出力される」という性質を持つ（これも鳩の巣原理の一種＝部屋数が有限のホテルに無限羽の鳥を入れると、少なくとも一つの部屋には無限羽の鳥が入る）。そのような出力値を c とすればいい。□
2020-06-20 野蛮な解法（その22） ペル方程式は非自明な解を持つ（後編）
d を平方数ではない正の整数とする。
定理 任意のペル方程式 x2 − dy2 = 1 は、非自明な整数解を持つ。
コメント x2 − dy2 = c （△）
が無限個の正の整数解を持つような整数 c が存在する。これについては「補題2」として前回証明した。c = 1 にできれば、これは「任意のペル方程式は無限個の解を持つ」という証明に他ならない。もし仮に（△）の両辺を c で割ることができれば、それを実現できそうだが、前回の例1…
x = 4, y = 1 ⇒ x2 − 10y2 = 6
x = 16, y = 5 ⇒ x2 − 10y2 = 6
x = 136, y = 43 ⇒ x2 − 10y2 = 6
…を眺めると、c = 6 で割り切れそうにない。でも「6 で割ると…」いう視点で見ると、3種類の x は、どれも「6 で割ると 4 余る」という共通点を持っている。y については「6 で割ると 1 または 5 余る」という特徴が見られるが、どうせ（△）では平方されるのだから y の正負はどちらでもいいわけで、第2解の y を y = −5 と書き直すと、どの y も ≡ 1 (mod 6) である。いずれにせよ、この例では mod 6 が糸口になりそうだ。この発想を一般化できるか？
定理の証明 （△）右辺の c は整数であり、負かもしれない。負の割り算は見通しが悪いので、| c | = m と置いて mod m で考えよう。もし仮に c = 0 なら √d が有理数になってしまい事実に反するので、c ≠ 0 であり m は正（0 ではない）。
√d が有理数ではないこと（容易に証明可能）は、既知の事実とする。
補題1・補題2の流れから、x, y はどちらも正の整数であるとしていい（2乗されるのだから、もちろん符号は関係ないのだが、話をシンプルにするため、正の整数に限る）。
1段に m 個の箱（引き出し）があり、それが m 段重なっている整理棚を作る。左下隅の箱を [0, 0]、右上隅の箱を [m−1, m−1] と名付け、（△）の解 (x, y) を一つ一つ分類して、対応する箱 [x′, y′] に入れる。ここで x′ は x を m で割った余り、y′ は y を m で割った余り。（△）の解は無限個ある。遅くとも、最初の m2 + 1 個の解を分類し終わった段階で、少なくともどれか一つの箱には2種類以上の解が入っている。[r, s] をそのような箱とし、そこから2組の解を取り出す。それらは次の形を持つ:
（△）の一つの解 (x1, y1) = (a1m + r, b1m + s)
上記と異なる（△）の解 (x2, y2) = (a2m + r, b2m + s)
ここで a1, a2 は、それぞれ x1, x2 を m で割ったときの商（余りを無視した整数）。b1, b2 は、それぞれ y1, y2 を m で割ったときの商（余りを無視）。
従って x2 − x1 = (a2 − a1)m。簡潔化のため A = a2 − a1 と書くと
x2 = x1 + Am
となり、同様に B = b2 − b1 と書くと
y2 = y1 + Bm
となる。m = ±c なので、この複号がマイナスの場合には A, B の符号を逆にすれば、結局:
x2 = x1 + Ac, y2 = y1 + Bc （⌘）
例2 冒頭の c = 6, d = 10 の3解のうち、(4, 1) と (136, 43) が同じ箱 [4, 1] に入る。この2種類の解について、x-座標同士の差 A = 136 − 4、y-座標同士の差 B = 43 − 1 は、いずれも c の倍数。
さて、ここからどう進めましょう？ （⌘）は「mod m で合同なものは差が m の倍数」という当たり前のことを言っているだけで、大した情報を与えてくれない。例えば2個の点 (4, 1), (136, 43) についての決定的な条件は、それらが（⌘）だけでなく（△）を満たすこと。何らかの形で条件（△）を使う必要があるはず。
補題2の証明では x ± y√d の形の無理数が重要な役割を果たした。補題1でも、この形の数が鍵となった。この形の数は（△）と密接に関係している。実際（△）を
(x + y√d)(x − y√d)
= c （△△）
と書き直すことができる。以下でもこの形を使う。
まず x2 + y2√d に（⌘）を代入する:
x2 + y2√d = (x1 + Ac) + (y1 + Bc)√d = (x1 + y1√d) + c(A + B√d)
この左端は x2, y2 の式。右端を x1, y1 の式に統一するため、（△△）の関係…
c = (x1 + y1√d)(x1 − y1√d)
…をそこに代入すると:
x2 + y2√d = (x1 + y1√d) + [(x1 + y1√d)(x1 − y1√d)](A + B√d)
= (x1 + y1√d)[1 + (x1 − y1√d)(A + B√d)]
= (x1 + y1√d)[1 + Ax1 + Bx1√d − Ay1√d − Bdy1]
= (x1 + y1√d)[(1 + Ax1 − Bdy1) + (Bx1 − Ay1)√d]
簡潔化のため
p = 1 + Ax1 − Bdy1, q = Bx1 − Ay1 （♫）
と置くと:
x2 + y2√d = (x1 + y1√d)(p + q√d)
同様に x2 − y2√d に（⌘）を代入するところから始めると、上記と同様の計算の結果:
x2 − y2√d = (x1 − y1√d)(p − q√d)
辺々掛け合わせて（△△）の関係を使うと:
(x2 + y2√d)(x2 − y2√d) = (x1 + y1√d)(x1 − y1√d)(p + q√d)(p − q√d)
c = c(p2 − dq2)
c ≠ 0 なので、両辺を c で割ることができる！
1 = p2 − dq2
これは、ペル方程式 x2 − dy2 = 1 が x = p, y = q という解を持つことを意味している。（♫）から分かるように、これは整数解。そして自明解 p = 1, q = 0 ではない（もしそうなら、
x2 + y2√d = (x1 + y1√d)(p + q√d)
なので、(x1, y1) = (x2, y2) になってしまうが、これは「箱から相異なる2解を取り出した」という前提と矛盾する）。同様に、自明解 p = −1, q = 0 でもない（もしそうなら (x1, y1) = (−x2, −y2) になってしまうが、これは「x1, y1 が正の解だけを箱に入れる」という前提と矛盾する）。従って、ペル方程式は非自明な解を持つ。□
例3 例2において x2 + 10y2 = 6 の2解 (4, 1) と (136, 43) は、同じ箱に入っている。
136 − 4 = 22 × 6 ⇒ A = 22
43 − 1 = 7 × 6 ⇒ B = 7
（♫）によると
p = 1 + 22 × 4 − 7 × 10 × 1 = 19
q = 7 × 4 − 22 × 1 = 6
は、x2 − 10y2 = 1 の解。検算: 192 − 10 × 62 = 361 − 360 = 1
つたない表現だが「2個の点の比」（のようなもの）がペル方程式の解になってる…
2020-07-12 野蛮な解法（その22-A） ペル方程式には解がある（証明の改良）
定理 d を平方数ではない正の整数とする。任意のペル方程式 x2 − dy2 = 1 は自明な解 x = ±1, y = 0 を持つが、それ以外にも解を持つ。
…既に「その22」でこれを証明したが、もっと良い方法があった。
〔定理の証明〕 まず補題2を引用する:
x2 − dy2 = c （△）
が無限個の整数解を持つような整数 c が存在する。
（△）には、次の条件を満たす2種類の整数解が存在する（この部分は易しい。その22の冒頭参照）:
c = x12 − dy12 = (x1 + y1√d)(x1 − y1√d) ←第1の解【あ】
c = x22 − dy22 = (x2 + y2√d)(x2 − y2√d) ←第2の解【い】
(x2, y2) ≠ (±x1, ±y1) ←2解は相異なる・符号だけを変えたものでもない【う】
x1 ≡ x2, y1 ≡ y2 (mod | c |) ←2解が満たす重大な付帯条件【え】
旧証明では、正の整数解の範囲で考えた。今回の証明も同じ範囲で進める方が簡単だが、あえて「正の」という限定を外すことにする。正の整数解だけでも無限個あり、上記が成り立つのだから、正に限らない範囲では当然これらが成り立つ。
うまく c で割り切れる形にしたいのだが、その22では見通しが悪かった。改善案は…【あ】と【い】の積は c2 なので、【あ】と【い】の右端の因数を組み合わせて c の倍数を作れるかも？ 例えば【あ】右端の第2因子と【い】右端の第1因子を掛け合わせると
(x1 − y1√d)(x2 + y2√d) = (x1x2 − y1y2d) + (x1y2 − x2y1)√d 【お】
【お】右辺の構成要素について【え】【あ】を使うと
x1x2 − y1y2d ≡ x1x1 − y1y1d = x12 − dy12 = c は c で割り切れ、
さらに x1y2 − x2y1 ≡ x1y1 − x1y1 = 0 も c で割り切れる！ ゆえに
P = (x1x2 − y1y2d)/c, Q = (x1y2 − x2y1)/c 【か】
と置くと、これらはいずれも整数で、【お】はこうなる:
(x1 − y1√d)(x2 + y2√d) = (cP) + (cQ)√d = c(P + Q√d) 【き】
同様に【あ】右端の第1因子と【い】右端の第2因子を掛け合わせると:
(x1 + y1√d)(x2 − y2√d) = c(P − Q√d) 【く】
【き】と【く】の積は、【あ】と【い】の積 c2 に等しい:
c2(P + Q√d)(P − Q√d) = c2 つまり P2 − dQ2 = 1
これは、ペル方程式 x2 − dy2 = 1 が解 x = P, y = Q を持つことを意味する。下記の理由から、(P, Q) は自明解 (±1, 0) ではない。□
〔例1〕 x2 − 10y2 = 6 について（c = 6, d = 10）、x1 ≡ x2, y1 ≡ y2 (mod 6) などの条件を満たす2解の例は、(x1, y1) = (4, 1), (x2, y2) = (136, 43)。【か】を適用すると P = (4 × 136 − 1 × 43 × 10)/6 = 114/6 = 19 となり、確かに c での割り算が成り立つ。同様に Q = 6。(P, Q) = (19, 6) はペル方程式 x2 − 10y2 = 1 の解。
〔例2〕 例1で (x2, y2) = (16, −5) としてもいい。結果は (P, Q) = (19, −6)。
(P, Q) が自明解 (±1, 0) ではないこと Q ≠ 0 を示す。もし仮に Q = 0 なら【か】より x1y2 − x2y1 = 0。簡単な式変形によると x1/y1 = x2/y2。これは「分数の約分ができる」ということ: 整数 k があって x2 = kx1, y2 = ky1 と書けるか、または x1 = kx2, y1 = ky2 と書ける。一般性を失うことなく、前者だと仮定できる（後者だった場合、単に変数名を入れ替えればいい）。ゆえに【い】x22 − dy22 = c は (kx1)2 − d(ky1)2 = k2(x12 − dy12) = c を含意する。この最後の丸かっこ内は【あ】より c に等しい。つまり k2(c) = c。従って k = ±1 または c = 0 だが、どちらも不合理。実際、k = ±1 は【う】に反するし、c = 0 なら、それを【あ】に代入すると d = x12/y12 となり √d = ±x1/y1 は有理数。これも事実に反する。□
「その20」の補題1（ディリクレの鳩の巣不等式）の段階から「x, y は互いに素」と制限できれば、この証明はさらに簡潔になる。けれど x2 − 23y2 = 8 の解 (10, 2) のように、一般論としては互いに素でないケースを考えたいこともあるかもしれない。
【か】の P, Q は、旧証明（♫）の p, q と等価。このことを示すには、旧証明における関係 A = ±(x1 − x2)/m, B = ±(y1 − y2)/m, m = ±c を確かめて（複号は c の正負による）、それらを（♫）の2式に代入すればいい。
2020-07-14 野蛮な解法（その22-B）
改良版の証明（22-A）では、Pete L. Clark の『Number Theory: A Contemporary Introduction』
https://web.archive.org/web/20181227114755/http://math.uga.edu/~pete/4400FULL2018.pdf
定理7.6を参考にした。この「数論: 現代的入門」は、自由にダウンロードできる形で一般公開されているのだが、その内容やスタイルはかなりユニーク。第7章ではペル方程式を扱っているが、お約束の連分数展開はバッサリ切り捨て、最初から二次体との関係を問題にする。解の存在・一般解の形は証明するものの、D = 61 のケースをどうやって解くか？といった話は一切しない。
完成度が高いとは言えない部分もある。今回の場合、「符号だけを変えたものを除外する」という【う】の条件が抜けているのが、論理的な穴になっていた。
例えば X2 − 7Y2 = 2 を考えると、相異なる2解 (X1, Y1) = (3, 1), (X2, Y2) = (−3, −1) は X1 ≡ X2, Y1 ≡ Y2 (mod 2) を満たすが、X1/Y1 = X2/Y2 となってしまい、これではペル方程式の非自明解の存在を導けない。このような穴をなくすためには「相異なる2個の整数解」では不十分であり、「絶対値が同じで符号だけ違う解については、異なると見なさない」ことが必要。
「正の整数解」に限るなら「相異なる」で足りる。定理を証明するだけなら、それが早道だろう。「正の」という限定をしたくないなら、上記のように「相異なる」の意味を調整する必要がある。
他にも α*conj(β) を計算すると言いつつ実際には conj(α)*β を計算していたり、変数名が D なのか d なのか不統一だったり、誤字脱字が散在したり、要するに「推敲された完成原稿」ではなく「講義で話す予定のことを書いてみたメモ」なのだろう。数学の教科書にありがちな「ひたすらエレガント＆クール」というものではなく、おしゃべりのようなノリで書かれている。刺激的でとがった話題も多い。ペル方程式の章で、超越数の存在に話が飛んだり…。だから悪い、という意味ではない。逆にこういうのが好き。「フェルマーのクリスマス～」も、これに刺激され、これを第1の参考文献としたものだった。
2020-06-21 野蛮な解法（その23） ペル方程式の解の数
ペル方程式が無限個の解を持つのは事実だが、前回の「解の存在証明」の方法によって「無限個の解の存在」まで証明できるだろうか。無限個の異なる「点」が入っている箱。そこから任意の組み合わせの2個の「点」を取り出して、ペル方程式の解を作れる。原料の2個の「点」が異なれば、一般には、異なる解が得られる。原料は無限種類あるのだから、得られる解も無限種類だと思えるが…。具体例で試してみる。
x2 − 10y2 = 6 の正の整数解（小さい順に8番目まで）:
P1 = (4, 1) ⇒ a = 0, b = 0
P2 = (16, 5) ⇒ a = 2, b = 0
P3 = (136, 43) ⇒ a = 22, b = 7
P4 = (604, 191) ⇒ a = 100, b = 31
P5 = (5164, 1633) ⇒ a = 860, b = 272
P6 = (22936, 7253) ⇒ a = 3822, b = 1208
P7 = (196096, 62011) ⇒ a = 32682, b = 10335
P8 = (870964, 275423) ⇒ a = 145160, b = 45903
奇数番目の解は ≡ (4, 1) mod 6、偶数番目の解は ≡ (4, 5) mod 6。a は x/6 の整数部分、b は y/6 の整数部分。
相異なる合同な2解 (x1, y1), (x2, y2) を選び、x1/6, y1/6, x2/6, y2/6 の整数部分がそれぞれ a1, b1, a2, b2 だとすると、
p = 1 + Ax1 − 10By1, q = Bx1 − Ay1
はペル方程式 p2 − 10q2 = 1 の解。ここで A = a2 − a1, B = b2 − b1。
P1, P3 から A = 22, B = 7:
p = 1 + 22 × 4 − 10 × 7 × 1 = 19, q = 7 × 4 − 22 × 1 = 6 は、確かにそのペル方程式を満たす（基本解）。
逆順に解を選択した場合、P3, P1 から A = −22, B = −7:
p = 1 + (−22) × 136 − 10 × (−7) × 43 = 19, q = (−7) × 136 − (−22) × 43 = −6。上記・基本解の符号を変えただけ。
P1, P5 から A = 860, B = 272:
p = 1 + 860 × 4 − 10 × 272 × 1 = 721, q = 272 × 4 − 860 × 1 = 228 は、確かにそのペル方程式を満たす（第2解）。
P3, P5 から A = 838, B = 265:
p = 1 + 838 × 136 − 10 × 265 × 43 = 19, q = 265 × 136 − 838 × 43 = 6 は、P1, P3 からと同じ（基本解）。
P3, P7 から A = 32660, B = 10328:
p = 1 + 32660 × 136 − 10 × 10328 × 43 = 721, q = 10328 × 136 − 32660 × 43 = 228 は、P1, P5 からと同じ（第2解）。
P1, P7 から A = 32682, B = 10335:
p = 1 + 32682 × 4 − 10 × 10335 × 1 = 27379, q = 10335 × 4 − 32682 × 1 = 8658 は、確かにそのペル方程式を満たす（第3解）。
P2, P4 から A = 98, B = 31:
p = 1 + 98 × 16 − 10 × 31 × 5 = 19, q = 31 × 16 − 98 × 5 = 6 は、P1, P3 からと同じ（基本解）。
選択の仕方によっては既出の解が反復して現れることもあるが、一つの（小さい）解を固定して、他方の解としてだんだん大きな解を選ぶなら、ペル方程式の新しい（大きい）解がいくらでも得られる。この方法だけで、事実上、無限個の解があると言えそうだ。しかし厳密に「あるところまでいったらパターンが循環して、有限個の解しか得られない」といった可能性を否定するためには、「前記の方法では p, q が単調に増加して、微妙な不等式の列を満たす」ことを示す必要がある。そんな苦労をしなくても「ペル方程式の非自明解は、1個あれば無限個ある」は簡単な方法で直接示せるので、それでいいでしょう。
追記 「p, q の増加列があること」「正の A, B をいくらでも大きくできること」を示すのは難しくないが、A, B は独立に大きくできず、これらが大きくなるとき、事実としては A/B は、より優秀な「√10 の近似分数」になっている。ペル方程式の性質からして驚くには当たらないし、どんどん誤差が小さくなるような √10 の近似分数の列があるのも事実だが、それは別の定理。「ついでにチャチャッと証明」というわけにはいかないだろう。
「野蛮な解法」第5部 バニラとチョコレート
2020-06-22 野蛮な解法（その24） d = 109 を1秒で
あえて普通の解法を封印し、ブルートフォースによる直接検索で
x2 − dy2 = 1 （●）
の正の整数解を求める。ここで「ブルートフォース」とは、上の式の x なり y なりに 1, 2, 3, … を入れ、方程式の解をしらみつぶしに探すこと。
d = 61 のケースは少し難しく、d = 109 のケースは、無理やり何とかしたものの、ブルートフォースによる攻略はこの辺が限界かと感じられた。この壁を越えるため、新しい手法を試してみたい。それは次の観察に基づく:
定理（○） 方程式 X2 − dY2 = −1 が正の整数解 X = A, Y = B を持つならば、
x = A2 + B2d, y = 2AB （＊）
は（●）の解である。
証明 仮定により: A2 − dB2 = −1
左辺を分解すると: (A + B√d)(A − B√d) = −1
その両辺を平方する: (A2 + B2d + 2AB√d)(A2 + B2d − 2AB√d) = (−1)2
（＊）の置き換えによれば: (x + y√d)(x − y√d) = 1
ゆえに x2 − dy2 = 1 が成立。□
このことから、（●）の右辺が 1 になってくれなくても、−1 になってくれれば（●）の解が得られる。
さっそく試す（PARI）。test は、与えられた数をテストする。（●）の解が見つかったら終了。done を呼ぶ（解を整形表示して1を返す）。（○）の解が見つかったら good を呼ぶ。解を表示し、定理によって（●）の解も見つかる。pell2 が検索の本体。どこかで test が成功すると、そこから1が返り、ループから復帰。成功しないままループが終わると、失敗メッセージを表示。
test(y,d) = { if( issquare(d*y^2+1,&x), done(x,y,d),
issquare(d*y^2-1,&x), good(x,y,d) ); }
done(x,y,d) = { print(x "^2-" d "*" y "^2=" x^2-d*y^2); 1; }
good(A,B,d) = { print1(A "^2-" d "*" B "^2=" A^2-d*B^2 " : ");
done(A^2+B^2*d, 2*A*B, d); }
pell2(d,yMax=10^6) = { print1("d=" d " : ");
for( y=1, yMax, test(y,d) && return );
print("\n### failed ###"); }
試しに pell2(61) とすると、一瞬で
d=61 : 29718^2-61*3805^2=-1 : 1766319049^2-61*226153980^2=1
と表示される。恐る恐る pell2(109) とタイプして [Enter] を押すと…おおっ、1秒もかからず
d=109 : 8890182^2-109*851525^2=-1 : 158070671986249^2-109*15140424455100^2=1
が出力された。d = 109 が、ブルートフォースの軍門に下った…。このケースを最初に解いたときは、めちゃくちゃ苦労したのだが…それが今や、何の最適化もしていないベタなブルートフォースで（偶数・奇数の区別さえしていない）、あっけなく解けてしまった。「負のペル方程式が解を持てば」という条件付きだが、その条件を満たすとき、大ざっぱに半分の桁数で（例えば20桁の解が10桁までの検索で）見つかる！
もちろん、この方法だって、結局は壁にぶち当たる…。偵察のため
for(d=1,100,issquare(d)||pell2(d))
for(d=101,200,issquare(d)||pell2(d))
とすると、d = 139, 151, 163, 166, 172, 181, 199 の7カ所でつまずく。このうち幾つかは検索空間の上限 10^6 を少し上げるだけで解決するし、扱う数は d = 61 のケースより小さいことが多く、やればできるだろう。
偵察を続けると、d = 331 が困難そうだ。d = 109 の場合より解が約18倍大きく、しかも今回のアルゴリズムによるショートカットが利かない。そして d = 379 のケースは、それよりさらに5倍大きい。
d = 181 は新記録の巨大解になるが、さいわい（○）が有効なケース。上記のテスト実装は遅いが（分オーダーの時間がかかる）、既に解は得られる:
pell2(181,oo)
d=181 : 1111225770^2-181*82596761^2=-1 : 2469645423824185801^2-181*183567298683461940^2=1
原理のテストは成功なので、後は高速化。最低でもあと10倍くらい速くしたい。
単純な事実（○）を再発見・利用することで、ブルートフォースで攻略可能な予想範囲が、少なくとも d = 330 まで広がった（今までは d = 181 が越えられない壁だった）。素直に連分数を使っていたら、この変な楽しさは味わえなかった（笑）。数学的な面では「d が4で割って1余る素数のとき、新記録が出やすい」という観察（その11参照）の意味も、見えてきた。大ざっぱに、4で割って1余る・3余る⇒mod d で−1が平方剰余・非剰余⇒○に解あり・なし⇒ペルの解が2乗のオーダー・1乗のオーダー⇒桁数が2倍になる・ならない…という流れだろう。
2020-06-24 野蛮な解法（その25） d = 109 を0.1秒で
【1】 x2 − 109y2 = −1 （ア）ならば x2 ≡ −1 (mod 109) （イ）。109 はバニラ素数なので（イ）は解を持つ。ここで「バニラ素数」とは、4 で割って 1 余る素数のこと。4 で割って 3 余る素数は「チョコ（レート）素数」。これらは正式な数学用語でなく、他の場所で使うとばかにされるので注意（笑）。なぜバニラだと上記が成り立ち、チョコだと駄目なのか？ これは「第一補充法則」と呼ばれる現象で、簡単に証明できることだが、話が拡散して本題がぼやけるので、今は暫定的にこれを事実と認める（次回以降に一段落したら証明する予定）。
具体的に（イ）の解は何か。怠惰に全数検索…
for(x=0,108,x^2%109==108&&print(x))
とタイプすると（数学者のひんしゅくを買うかもしれないが）、一瞬で x = 33, 76 が得られる。ゆえに（ア）を満たす x があれば、それは 109 の倍数プラス 33 or 76。従って、検索方法は:
T(x,d) = { yy=(x^2+1)/d; issquare(yy) && good(x,sqrtint(yy),d); }
forstep(n=0,oo,109, if(T(n+33,109)||T(n+76,109), break))
good と done は前回のものを再利用。実行すると、約0.1秒で
8890182^2-109*851525^2=-1 : 158070671986249^2-109*15140424455100^2=1
が出力される。
【2】 x, y が偶数か奇数か考えると、さらに2倍高速化できる。mod 4 において x ≡ 0, 1, 2, 3 のとき、それぞれ x2 ≡ 0, 1, 0, 1。d ≡ 1 つまりバニラなら、y ≡ 0, 1, 2, 3 のとき、全く同様に dy2 ≡ 0, 1, 0, 1。ゆえに x2 − dy2 ≡ −1 になるのは、その左辺第1項が 0、第2項が 1 の場合、すなわち x ≡ 0 or 2 かつ y ≡ 1 or 3 の場合に限られる。従って（イ）のような条件に加え、x が偶数の場合だけ考えれば十分。
forstep(x=0,109*2-1,2,x^2%109==108&&print(x)); （ウ）
とタイプすると、x = 76, 142 が得られ、
forstep(n=0,oo,109*2, if(T(n+76,109)||T(n+142,109), break)) （エ）
とすれば、【1】と同じ解が半分の時間で得られる。前回の「ベタな全数検索」に比べると、ざっと100整数につき1個だけテストするので、100倍オーダーの高速化となる。同じ時間で100倍の範囲を探せるので、検索上限を 10^6 から 10^8 にしてもよく、成功確率の向上が期待される。
ところで d がチョコ素数のとき（イ）従って（ア）は絶対に解を持たないので、（○）検索は不毛であり、無駄に所要時間が倍増する。一方、d が奇素数のとき（○）が解を持つとすれば d はバニラだが、その逆を証明したわけではない。原理的には「バニラ素数だが○の解がない」という可能性がある。そんな事態が起きたら、そのときはそのときでまた考えよう。
【3】 d が与えられるたびに（ウ）（エ）を手動でタイプ・計算するのは面倒なので、自動化しておく。ついでに mod 4（結果的には偶数・奇数による縛り）の代わりに mod 8（結果的には4で割った余りによる縛り）を使うことで、さらに2倍の高速化を行う。d がバニラ素数 p の場合、x2 − py2 ≡ −1 が成り立つためには、p が8で割って1余るなら x は4で割り切れ（余り r=0）、p が8で割って5余るなら x は4で割って2余る（r=2）…という性質（容易に確かめられる）に基づく。
prep2a(p,r) = { v=[]; forstep(x=r, 4*p-1, 4, x^2%p==p-1 && v=concat(v,x)); v; }
pell2a(p,nMax=11*10^8) = { v = if( p%8==1, prep2a(p,0), prep2a(p,2) );
printf("[ p=%d ]\n", p);
forstep( n=0,nMax,p*4,
for(i=1,#v, T(n+v[i],p) && return) );
print(" ### failed ###"); }
注意 pell2a は、d がバニラ素数の場合専用のルーチンであり、チョコ素数、合成数などには対応できない（それらについては後から手を打つ）。
まずは小手調べとして pell2a(61) とタイプすると、1ミリ秒もかからず答えが出る。昔はやや手ごわい相手だったが、今は瞬殺！ 因縁の pell2a(109) についても、もはや0.1秒かからず答えが表示される。しかし問題はここから。
forprime(p=2,200, p%4==1 && pell2a(p))
とすると p = 181 で失敗する。個別に無限検索 pell2a(181,oo) を実行すると、1秒ほどで
1111225770^2-181*82596761^2=-1 : 2469645423824185801^2-181*183567298683461940^2=1
が出るので、苦戦というほどではないが…。同様に200から400の範囲のバニラ素数を試すと、p = 277 と p = 397 の2カ所で失敗した。それらについて個別撃破を試みると、数秒かかって
8920484118^2-277*535979945^2=-1 : 159150073798980475849^2-277*9562401173878027020^2=1
が得られ、約10秒かかって
20478302982^2-397*1027776565^2=-1 : 838721786045180184649^2-397*42094239791738433660^2=1
が得られた（工夫すれば、もう少し速くできると思われる）。
今回の成果。d が400以下のバニラ素数については、全てのペル方程式の基本解が（ブルートフォースで！）得られた。副産物として「ペル方程式の基本解・巨大新記録」が2回更新された。まず x=277 のとき、x が21桁（1垓5915京）、y が19桁（956京）。さらに x=397 のとき、x が21桁（8垓3872京）、y が20桁（4209京）。「d ≤ 400 におけるペル・チャンピオン」のタイトルを確定させるためには、バニラ素数以外の d についても調べる必要がある。
2020-06-25 野蛮な解法（その26） 京（けい）の上は垓（がい）
x2 − 109y2 = 1 を満たす最小の正の整数解は、
(x, y) = (158070671986249, 15140424455100)
ざっと150兆と15兆（15桁と14桁）。連分数を使わない素朴な解法（ブルートフォース）では、この辺が限界かと思われた（「その9」）。しかしアルゴリズムを少し工夫して、
x2 − dy2 = 1
の d が400くらいまでの範囲は、ブルートフォースで解けるようになってきた。d について、4で割って1余る素数の場合（バニラ素数）、4で割って3余る素数の場合（チョコ素数）、それら以外の場合（合成数など）の3パターンを区別する。ここで扱っている第1のパターンでは、第2のパターンと違い「その24」の定理（○）が利用可能。この違いは大きい。
前回は 8d を法とすることにより、単純全数検索と比べ 4d 倍くらいの高速化を図った。上記 d = 109 の例ならざっと400倍。今回は 24d を法とすることにより、さらに2倍程度の高速化を図る。
d がバニラ素数 p の場合、p ≡ 1 (mod 8) または p ≡ 5 (mod 8)。ところで、○の式 x2 − dy2 = −1 が成り立つとすれば、そのときには当然、両辺が ≡ −1 (mod 8) になる。この条件を満たす組み合わせを総当たり的に調べると、p ≡ 1 (mod 8) の場合、x は4の倍数でなければならず、p ≡ 5 (mod 8) の場合、x は4で割って2余る数でなければならない。同様に、バニラ素数は p ≡ 1 (mod 3) または p ≡ 2 (mod 3) のどちらかだが、○の式が ≡ −1 (mod 3) になるためには、前者の場合、x は3の倍数でなければならず、後者の場合、x は3の倍数以外でなければならない。
これらを組み合わせると、バニラ素数 p は p ≡ 1 (mod 24) または p ≡ 17 (mod 24) であるか（これらは p ≡ 1 (mod 8) の場合）、そうでなければ p ≡ 5 (mod 24) または p ≡ 13 (mod 24) であり（これらは p ≡ 5 (mod 8) の場合）、これら4種類の可能性について、x を24で割った余りは、それぞれ {0, 12}, {4, 8, 16, 20}, {2, 10, 14, 22}, {6, 18} のいずれかでなければならない。第1のケース {0, 12} は「4の倍数」かつ「3の倍数」。第2のケース {4, 8, 16, 20} は「4の倍数」かつ「3の倍数以外」。同様に、第3のケースは「4で割って2余る」かつ「3の倍数以外」、第4は「4で割って2余る」かつ「3の倍数」。
これに加えて、○の条件は x2 ≡ −1 (mod p) …言い換えれば x2 を p で割ると p−1 余ること…なので、検索対象は次のように絞られる（PARI）。
prep3a(p) = { my( r=p%24,
u = if(r==1,[0,12], r==17,[4,8,16,20],
r==5,[2,10,14,22], r==13,[6,18]),
k=#u, v=[] );
forstep( n=0, 24*(p-1), 24,
for( i=1, k, my(x=n+u[i]); if( x^2%p==p-1, v=concat(v,x); break ) )
); v; }
例えば prep3a(109) とタイプすると、こういう結果が出る:
[294, 1014, 1602, 2322]
24p = 2616 個の自然数につき、上記4種類の x だけを試せばいいので、単純に全部の自然数を片っ端から試すより600倍以上、効率が良い。具体的には 2616 の倍数のそれぞれに上記の4個の数を足したものを x として、○を y について解こうとする（整数解が見つかれば成功）。検索を行うコード例:
pell3a(p,N=21*10^9) = {
if(isprime(p)&&p%4==1, printf("p=%d\n ",p), return);
my( v=prep3a(p), k=#v );
forstep( n=0, N, 24*p,
for( i=1, k, T(n+v[i],p) && return ));
print("### failed ###");
}
T(x,d) = { issquare((x^2+1)/d) && G(x, sqrtint((x^2+1)/d), d); }
G(A,B,d) = { printf("%d^2-%d*%d^2=1\n", A^2+B^2*d, d, 2*A*B); 1; }
pell3a(109) とタイプすると、10ミリ秒のオーダーで
158070671986249^2-109*15140424455100^2=1
と表示される。さらに d ≤ 400 の範囲の全てのバニラ素数について、
forprime(p=2,400, p%4==1 && pell3a(p))
とタイプして一括検索すると、数秒で全部解決する。このうち一番難しいのはやはり d = 397 で、これだけで所要時間の半分以上を占める。pell3a 関数の N のデフォルト値 21*10^9 は、このケースが解けるギリギリに設定されている。…全部解けるともう少しやってみたくなるもの。
forprime(p=401,500, p%4==1 && pell3a(p))
とタイプすると、結果は…
p=401 801^2-401*40^2=1 p=409 ### failed ### p=421 ### failed ### p=433 104564907854286695713^2-433*5025068784834899736^2=1 p=449 71798771299708449^2-449*3388393513402120^2=1 p=457 ### failed ### p=461 1182351890184201^2-461*55067617520620^2=1
このうち p=409 については、無限検索コマンド pell3a(409,oo) を投入すると、10秒くらいかかるかもしれないが
25052977273092427986049^2-409*1238789998647218582160^2=1
が得られる。
x は23桁、250垓5297京7273兆924億2798万6049（ざっと25セクスティリオン）、y は22桁、約12垓（1.2セクスティリオン）。ペル方程式の巨大解・新記録（？）が更新された。p=457 はもう少し小さく、それぞれ69垓台と3垓台。残った p=421 については、最小の正の整数解が30桁以上なので、この実装では実用的でない。チョコレート素数 p=331 や p=379 のケース（それぞれ16桁、17桁）を解くのが先決だろう。p=409 の最終結果は23桁だが、内部的には○の12桁で処理しているので、見掛けほど難しくない（さらなる高速化の余地はある）。
2020-06-26 野蛮な解法（その26の補足）
速度的メリットはないが、prep3a を簡潔にこう書くこともできる:
prep3A(p) = { my( r=p%24,
u = if(r==1,[0], r==17,[4,8],
r==5,[2,10], r==13,[6]),
k=#u, v=[] );
forstep( n=0, 24*(p-1), 12,
for( i=1, k, my(x=n+u[i]); if( x^2%p==p-1, v=concat(v,x); break ) )
); v; }
prep3a と比べ、内側のループ回数が半分になる代わり、外側のループ回数が2倍になる。だから本質的な計算量に違いはなく、実測でも prep3a とほぼ同じ速度だった（どちらかと言うと、わずかに遅くなるかもしれない）。
prep3A では内側のループが1～2回しか回らず、1回しかしないことをループで書くのは無駄。一方、prep3a では、配列 u の長さが2倍になって無駄（12ずらして、同じ要素を2回ずつ書いている）。もともと物量攻撃のブルートフォースであり、そんな微視的な部分の無駄は大勢に影響ない。しかし r の値で場合分けして、内側のループをやめ、そこをベタで書いた4種類の prep3A を作れば、これらの無駄を排除できるだろう。もちろん2倍・3倍といった劇的な高速化は見込めない（ビット計算量が減るわけでなく、単に制御構造の問題なので）。「コード量が4倍になって1%の高速化」みたいな、せこい話になるだろう。
2020-06-30 野蛮な解法（その27） 普通に解いてもつまらない！ 数学の縛りプレイ
x2 − 214y2 = 1 のような形の2変数方程式（「ペル方程式」と呼ばれる）の、正の整数解を求めたい。パズル的に面白いだけでなく、数学的にもいろいろと重要な意味を持つのだが、とりあえず意味なんか気にせず、遊んでみる。
x2 − dy2 = 1 の d（平方数ではない自然数）が小さい場合、大抵は簡単に解を見つけることができる。例えば
x2 − 11y2 = 1
の最小解は、容易に暗算可能。左辺に y = 1, 2, 3, … を入れると y = 3 のとき第2項は −99 になるが、そのとき x2 = 100 となり、それを満たす正の整数 x は、もちろん10。次の例題も、同様にして、暗算で簡単に解けるので、やってみよう。
x2 − 12y2 = 1
ペル方程式の解法としては、平方根の近似分数が利用されることが多い。ここでは普通の道を歩まず、上記の暗算のような「全数検索」的アルゴリズム pell4(d) をテスト実装した。pell4 は、まず小さい y について、上の例そのままの素朴な検索を試み、それでは解が見つからない場合、d が 4k+1型の素数なら pell4a(d)、d が 4k+3型の素数なら pell4b(d)、それ以外の場合、pell4c(d) をそれぞれ呼び出す。
pell4a は「その26」の pell3a の改良版。解きたい方程式の右辺を −1 に置き換えたもの（負のペル方程式）の解を見つけることで、間接的に与式を解く。まず 64d を法として解 x になり得る数を決定。それらの数について、1ステップ 64d の「ホイールシーブ」のような検索を1000万程度まで行う。これが不成功なら、今度は 64⋅81⋅25⋅49d を法として1200億まで検索。…「ホイールシーブ」の準備では、使いたい法の因子（素数べき）について、それぞれを法として「負のペル」の解 x になり得る数を決定し、得られた複数の条件に中国剰余定理を適用して、大きな法における1種類の条件にまとめておく。法 d に関しては、q2 ≡ −1 (mod d) を満たす q を求めることに相当。これも全数検索可能だが、数論的には、原始根の (d − 1)/4 乗を q として {q, q3} が求めるもの。
pell4b も pell4a とほぼ同じ実装。こちらは「負のペル」の解を探さず、直接与式を解く。まず法 256d で1億まで。それで駄目なら、法 256⋅81⋅25⋅49d で2800億まで。q2 ≡ 1 (mod d) を満たす q が必要だが、それは ±1 (mod d) に他ならない。
最後の pell4c は、合成数を扱うため、ややこしい。現在のテスト実装では、まず 128⋅81q を法として、5億まで検索。ここで q は、d を2と3でそれぞれ割れるだけ割った結果の数（言い換えれば、2とも3とも互いに素であるような q の約数うち、最大のもの）。それで駄目なら、今度は対応する「負のペル」が解を持ち得るか調べ、持ち得るならそれを探す。それも駄目な場合、法 256⋅81⋅25⋅49q で検索。この q は、d を2、3、5、7で割れるだけ割った結果の数。
以上によって、330以下の全ての d について、ブルートフォースでペル方程式を解けるようになった。所要時間は、易しいケースなら一瞬、難しいケースでもおおむね1秒以内だが、d = 214 は数秒かかり、d = 301 はさらに遅い。この実装（特に pell4c）はまだ試作段階。うまくやれば、もっと速くなるかもしれない。楽に解きたいだけなら、もっと速いやり方が複数あるのだが、あえてブルートフォースで。
pell4 がうまく対応できない最初のケース（時間がかかり過ぎる）は d = 331, 379（「その24」参照）。壁になっている d = 331 では、探したい数が16桁。今の速度（想定ターゲット最大13桁くらい）では、きつい。なにしろ全数検索なので、1桁増えるごとに探す場所が10倍、3桁増えれば1000倍…それだけ時間がかかる。「どこかで探し切れなくなる」という結末を承知の上で、あの手この手の高速化を行い、どこまで追いすがれるか…
pell4c の d = 301, 334 も時間がかかるが（13～14桁）、数十秒～5分くらい待てば解が得られるので、一応、実用になる（速度的には改善の余地あり）。pell4a が実用的でなくなる最初の数は d = 421 で、探したい数は17桁（「その27参照」）。今はまだ330だし、この全数検索はただの遊びなので、多分そのうち飽きる。本当の目標は、とりあえず一般化ペル方程式の基本解の上限について、すっきりさせること。
2020-07-01 野蛮な解法（その28） d=331の壁を突破
ペル方程式 x2 − dy2 = 1 について、普通に解かず、ゲーム感覚で全数検索をしてきた。d が4k+1型の素数なら便利なショートカットが存在するのだが、4k+3型の素数のときには同じショートカットを使えず、4k+3型の d = 331 が壁となっていた。何時間も計算すれば解は得られるだろうが、4k+3型素数に使えるうまいショートカットは無いものか…？
4k+1型のショートカットは、与式の右辺を −1 に置き換えた「負のペル」からの間接攻略だった。同様に、与式の右辺を −2 に置き換えた次の方程式を解くことで、d = 331 に対応できる:
A2 − dB2 = −2 （♡）
仮に（♡）を満たす正の整数 A, B が得られた場合、（♡）の両辺を2乗して
(A2 − dB2)2 = 4
(A2)2 − 2(A2)(dB2) + (dB2)2 = 4
(A2)2 + 2(A2)(dB2) + (dB2)2 − 4(A2)(dB2) = 4
(A2 + dB2)2 − d[4(AB)2] = 4 （♡♡）
になるが、（♡♡）の右辺と、左辺第2項・角かっこ内は、どちらも4の倍数なので、当然、左辺第1項も4の倍数。だから各項を4で割ることができて、
C − d(AB)2 = 1
の形になる。ここで C は、（♡♡）の左辺第1項を4で割った結果の整数だが、もしそれが平方数 C = x2 になってくれれば、後は y = AB と置くだけで、ペル方程式が解けたことになる。このアイデアは実を結ぶ。実際、（♡）より
dB2 = A2 + 2
であり、これを使うと、（♡♡）の左辺第1項は
(A2 + A2 + 2)2 = (2A2 + 2)2 = [2(A2 + 1)]2 = 4(A2 + 1)2
すなわち、それを4で割った結果 C とは (A2 + 1)2。結局、（♡）が解けるなら
x = A2 + 1, y = AB （＊）
が、対応するペル方程式の解。
d = 331 の場合、（♡）は A = 52778687, B = 2900979 という解を持つ（これは7～8桁なので、容易に全数検索可能）。これらを（＊）に入れると:
x = 2785589801443970,
y = 153109862634573
が、ペル方程式 x2 − 331y2 = 1 の解。検算すると、確かに合っている！
たった16桁、1000兆のオーダーだが、直接検索は難しい。d が4k+1型素数のとき、すぐに解が京や垓のオーダーになるが、あれは一種の「バブル」。内部的な検索は10桁やそこらなのだが、最終結果ではそれが大ざっぱに平方されるため、見掛け上、20桁（1000京）オーバーも珍しくない。現状、ブルートフォースで直接検索できているのは、秒単位では12～13桁くらいまで。分単位では14桁～ギリギリ15桁。あとたった1桁とはいえ、ブルートフォースでは、それは計算量的に10倍。「ギリギリできていることの10倍」は、きつい。並列処理とか、強力なCPUとか、高速なコンパイル言語とか、そういう手もあるけど…
d = 331 については
x = 8278, y = 455 ⇒ x2 − 331y2 = 9
という関係も見つかったが、これは（少なくとも簡単な方法では）役立ちそうにない。それと、一つ宿題が残っている。d = 331 のペル方程式については、上記の通り、間接的な方法で15～16桁の解を得たが、pell4 で直接検索したのは12桁まで。上記が基本解（最小の正の解）であることを言うためには「12～16桁のまだ探してない範囲に、これより小さい他の解がない」ことを示す必要がある。4k+1型関連の「巨大バブル解」にも、同様の問題が…。「12桁未満の解がなく16桁の解があるなら、それが最小解」というのは、ペル方程式の基本性質から分かることだが、このメモでは、まだそういう考察をしていない。
2020-07-02 野蛮な解法（その28の追記）
（♡）の方法は、右辺が +2 の場合にも同様に成立、（＊）が x = A2 − 1 になる。右辺 ±2 の形の（♡）を適宜使うことで、d が4k+3型素数のときの直接検索は必要なくなった。d = 331, 379 などはもはや問題でなく、4k+3型素数については、恐らく700台まで進むことができる。d が合成数でも、それを法として ±2 が平方剰余なら、同じ手を試みることが可能。その結果、速度的に遅かった d = 214, 334 もほぼ一瞬になった。新しい壁は、4k+1型素数の d = 421。
d = 421 以外の d については、525 まで全部対応できるが、合成数のケースで少し時間がかかるものが4個ある。d = 301（依然として10秒単位の時間がかかる）、d = 436（分単位の時間がかかる）など。d = 421 についても、ブルートフォースで解決できるめどが立った。
追記（2020年7月4日） 「時間がかかる4個の合成数」とは、301、436、501、508だったが、4の倍数については方法が見つかったので、436、508は解決。残った301、501についても「ケイリー法」で高速化できる見通し。現在、壁になっている合成数は649（これは恐らく対応可能）、721（合成数。恐らく対応可能）、769（4k+1型素数）。4k+3型素数については、1000まで到達可能と思われる。
2020-07-03 野蛮な解法（その29） ケイリーのノット
【1】 ペル方程式 x2 − dy2 = 1 を直接解かなくても、その右辺を −1 または ±2 にした式の解が見つかれば、間接的にペル方程式の解が得られる:
(I) A2 − dB2 = −1 ならば x2 − dy2 = 1 である。ここで x = A2 + B2d, y = 2AB。（「その24」参照）
(II) A2 − dB2 = −2 ならば x2 − dy2 = 1 である。ここで x = A2 + 1, y = AB。（その28）
(III) A2 − dB2 = 2 ならば x2 − dy2 = 1 である。ここで x = A2 − 1, y = AB。（その28の追記）
この手法により、d = 500 以上まで「ブルートフォースで」ペル方程式を解けるようになった…ただ一つの例外 d = 421 を除いて。x2 − 421y2 = 1 の最小解は34桁で、とてもブルートフォースでは解けない。対応する (I) の解も17桁で、ブルートフォースでは時間がかかり過ぎる。このケースでは (II), (III) の手法も使えない。
「ブルートフォースで」という縛りは「できるかな？」という遊びであり、普通の解法を使えば d = 421 が難しいわけではない。でも、たった一つの d のため「d = 500 までコンプ」できないというのも悔しい。このケースに使えるうまい手は無いか？
【2】 (I)～(III) を眺めていると、A2 − dB2 = ±4 からのショートカットも可能なのでは、という考えが浮かぶ。d = 421 については
4449392 − 421 × 216852 = −4
という事実もある。これを利用できないか？
19世紀の英国の数学者 Cayley（ケイリー）は、この種の状況について、こう指摘した [1]。
(IV) t2 − du2 = −4 ならば T2 − dU2 = 4 である。ここで T = t2 + 2, U = tu。
(V) T2 − dU2 = 4 ならば x2 − dy2 = 1 である。ここで x = (T3 − 3T)/2, y = (T2 − 1)U/2。
【3】 ケイリーのような大数学者は、この問題をどういう角度から眺めていたのだろうか？ それを紹介する前に、比較のため、素朴な方法で同じ結論を導いてみよう。(IV) は易しい。−4 を +4 にしたいのだから2で割って2乗するか、（同じことだが）2乗して4で割ればいい。
(t2 − du2)2 = 16 の左辺は
t4 − 2t2du2 + d2u4 = (t2 + du2)2 − 4t2du2
仮定により du2 = t2 + 4 なので上記右辺は
(t2 + t2 + 4)2 − 4t2u2d = [2(t2 + 2)]2 − 4(tu)2d = 4(t2 + 2)2 − 4(tu)2d
これが 16 に等しいというのだから 4 で割ると
(t2 + 2)2 − (tu)2d = 4
〔例1〕 d = 13, t = 3, u = 1, 32 − 13 × 12 = −4 のとき
T = 32 + 2 = 11, U = 3 × 1 = 3, 112 − 13 × 32 = 4
(V) については、仮定より dU2 = T2 − 4。その両辺に (T2 − 1)2 を掛けて
dU2(T2 − 1)2 = (T2 − 4)(T2 − 1)2 = (T2 − 4)(T4 − 2T2 + 1)
= T6 − 6T4 + 9T2 − 4 = (T3 − 3T)2 − 4
これを整理した (T3 − 3T)2 − d[(T2 − 1)U]2 = 4 の両辺を 4 で割れば、書かれている通りの結果を得る。□
〔例2〕 d = 13, T = 11, U = 3, 112 − 13 × 32 = 4 のとき
x = (113 − 3 × 11)/2 = 649, y = [(112 − 1) × 3]/2 = 180, 6492 − 13 × 1802 = 1
【4】 上記 (V) の証明は透明でない。なぜ dU2 = (T2 − 4) からスタートし、なぜ両辺に (T2 − 1)2 を掛けるのか。どの式変形も間違っていないが、必然性が感じられず、分かった気がしない。(IV) の証明にしても、
(a − b)2 = (a + b)2 − 4ab
のパターンは「よく使われるトリック」だが、なぜそれをここで使うのか。
ケイリー自身の考え方は、かみ砕いて言うと、次の通り。
t2 − du2 = (t + u√d)(t − u√d) という分解を念頭に置いて、(IV) の仮定を
(t + u√d)(t − u√d) = −4
と書く。その両辺を平方する。左辺の平方は次のようになり、これが右辺の平方 16 に等しい:
(t + u√d)2(t − u√d)2 = (t2 + du2 + 2tu√d)(t2 + du2 − 2tu√d)
= (2t2 + 4 + 2tu√d)(2t2 + 4 − 2tu√d) = (2t2 + 4)2 − (2tu)2d
2個目の等号では、du2 = t2 + 4 を使った。…最初の証明と同じようなものだが、トリッキーな面がなくなった。
(V) については、同様の手法が、証明の上で絶大な効果を発揮する。(T + U√d)(T − U√d) = 4 の両辺を立方すると:
(T + U√d)3(T − U√d)3 = 64 （ア）
左辺の前半について、機械的に計算すると:
(T + U√d)3 = 4(T3 − 3T) + 4(T2 − 1)U√d
途中経過は略したが、3乗の部分を展開して整数部分とそれ以外（√d を含む部分）をそれぞれまとめ、dU2 = T2 − 4 を使って整理しただけ。見やすくするため
A = 4(T3 − 3T), B = 4(T2 − 1)U （イ）
と書くと (T + U√d)3 = A + B√d となる。同様にして（ア）の左辺後半は
(T − U√d)3 = A − B√d となる。従って（ア）は、こうなる:
(A + B√d)(A − B√d) = A2 − B2d = 64
この真ん中の式について（イ）の省略記法を元に戻すと:
16(T3 − 3T)2 − 16[(T2 − 1)U]2d = 64
両辺を 64 = 16⋅22 で割ると:
[(T3 − 3T)/2]2 − d[(T2 − 1)U/2]2 = 1
すなわち x = (T3 − 3T)/2, y = (T2 − 1)U/2 と置けば x2 − dy2 = 1 である。
…淡々と計算するだけで、自然に結論が得られた！ 最初の証明（論理は分かるが、意味が分からない）とは大違い。
半面「整数解の話なのに、なぜ無理数が出てくるのか」「上記の計算は一体何を意味しているのか」「立方するという発想はどこから出てきたのか」といった別の疑問が生じるかもしれない。それを考えるためには少し準備が必要だが、複素数の計算で i を2乗すると整数 −1 になるように（そのため複素数の積において、実部と虚部が絡み合う）、ここでは無理数 √D を2乗すると整数 D になることが鍵を握っている。遊びのような具体例を通じて「普通の整数の背後にある少し深い世界」を垣間見るのも悪くないかと…。ケイリーは同じノットの中で、右辺 −4 のペル型方程式から、右辺 −1 のペル型方程式の解を導く方法にも言及している。
Cayley は英国人だが、英語だけでなく時々フランス語で論文を書いた。1857年のこの Note もフランス語（だから「ノート」でなく「ノット」）。19世紀にはまだ「共通語＝ラテン語」時代の余韻が強く、国際的な学術の世界でも、今より言語的多様性が大きかった。Note で言及されている「Degen の表」はラテン語の文献で、ペル方程式の解を d = 1000 まで収録。Cayley 自身、後にこれを拡張して d = 1001 から 1500 までの表を作った。
〔追記〕 Degen はドイツ生まれのデンマークの数学者。「デンマークの天才数学者デェーイェン」「八元数の掛け算表」参照。
[1] A. Cayley, Note sur l’équation x2 − Dy2 = ±4, D ≡ 5 (mod. 8)
http://www.digizeitschriften.de/dms/resolveppn/?PID=GDZPPN002149834 ［367～371ページ］
https://archive.org/details/collmathpapers04caylrich ［論文集・第4巻40～42ページ］
2020-07-07 野蛮な解法（その29の続き） d = 109 を1ミリ秒で
【5】 ケイリーの公式は、計算上は一般的に成立するが、実用上、ペル方程式の解法に役立つのは一定の場合に限られる。
〔悪い例〕 d = 11 とする。T2 − 11U2 = 4 の最小解は:
202 − 11 × 62 = 4 （ガ）
これに (V) を適用すると:
39702 − 11 × 11972 = 1 （ク）
（ク）は d = 11 の場合のペル方程式の解の一つだが、最小解
102 − 11 × 32 = 1 （カ）
とは異なる。右辺4の場合の最小解からスタートしたのに、右辺1の場合の最小解（カ）ではなく、ややこしい（ク）が出てきてしまった。何が悪かったのか。そもそも（ガ）の各項は 4 = 22 の倍数なのだから、単に（ガ）の両辺を4で割れば（カ）になっていた。(V) ではもともとの解が大ざっぱに3乗されるので、下手に使うと、最小解ではなく第3解が出てきてしまう。T = 20, U = 6 のように両方偶数の場合、この公式を使うのではなく、単にそれらの偶数を2で割ればいい（平方された結果としては、4で割ることに相当）。t, u が両方偶数の場合についても、(IV) により T, U が両方偶数になる。…ちなみに（カ）と（ク）の間にある第2解は:
1992 − 11 × 602 = 1 （キ）
〔別の種類の悪い例〕 d = 12 とする。T2 − 12U2 = 4 の最小解は:
42 − 12 × 12 = 4
これに (V) を適用すると:
262 − 12 × (15/2)2 = 1
…計算は合っているが、整数解を求めたいのだから、これでは役立たない。T, U がそれぞれ偶数・奇数の場合、(V) を適用すると y 側に2分の1の端数が生じてしまう。t, u がそれぞれ偶数・奇数の場合も (IV) を適用すると (V) において同じ状態に。では T, U ないし t, u がそれぞれ奇数・偶数 の場合はどうか。その場合、
T2 − dU2
のような式は「奇数マイナス偶数」で値が奇数。偶数である ±4 になり得ない。
結局、両方偶数も駄目だし、片方だけ偶数も駄目。ケイリーの公式が奏功するとすれば、T, U などが両方とも奇数の場合。例えば
2612 − 109 × 252 = −4
のような形が得られたなら、見込みがある。
【6】 上記のことから、公式の適用条件を限定できる。mod 8 で考えると、(V) が有用なのは
T2 − dU2 ≡ 4 (mod 8) （＊）
において T, U がどちらも 1, 3, 5, 7 (mod 8) のどれかになる場合だが、これら16通りの組み合わせに対して、d ≡ 0, 1, 2, …, 7 (mod 8) の8パターンを考えると、（＊）が成り立つ可能性があるのは d ≡ 5 (mod 8) のみ。(IV) についても同様。
〔証〕 仮定より T, U はどちらも ≡ ±1, ±3 (mod 8)、ゆえに T2, U2 はどちらも ≡ 1。従って（＊）は 1 − d ≡ 4 (mod 8) を含意する。つまり d ≡ 5。□
すなわち、ケイリーの公式が役立ち得るのは d ≡ 5 (mod 8) の場合のみであり、しかも、T, U などを奇数に限定できる（もし限定なしの検索を行い、右辺が ±4 の偶数解が得られた場合、ケイリーの公式を使わず直ちに各項を4で割れば、正または負のペル方程式の解が得られる）。
〔例〕 d = 37 ≡ 5 (mod 8) のとき、外形的には公式が役立ち得るが、実際には (IV), (V) は奇数解を持たない（ブルートフォースの立場からは、このことは「やってみないと分からない」）。検索を奇数に限定しない場合、
122 − 37 × 22 = −4
が見つかるが、この偶数解に (IV), (V) を適用すると、最初の「悪い例」と同様の結果になる。このような最小解が得られた場合、両辺を4で割り
62 − 37 × 12 = −1
としてから「負のペル」に対する通常の扱いを行えば、ペル方程式の最小解が得られる。この場合、最初から「負のペル」の検索を行っていれば、半分の検索範囲（半分の時間）で同じ結論が得られていた。
「奇数解が存在せず、偶数解なら存在するケース」があるのだから、両方探した方がいいのでは？という疑問が生じるかもしれない。マルチスレッドが許されるなら、±4 の奇数解と ±1 の整数解を平行して検索するのが最も効率的だろう。シングルスレッドの場合は微妙だが、奇数だけを探せば（奇数・偶数両方の検索に比べ）半分の時間で済む。2倍の高速化ができるポイントなので、できればそれを生かしたい。
普通のアルゴリズム（平方根の連分数展開を利用）では、上記の疑問は自然に消滅する（ケイリーのノットもそういう記述になっている）。ブルートフォースは「連分数というレーダー」を見ないで行うブラインド・アタックなので、この部分の見通しが利かず、何が最善なのか即断できない。
「8k+5型の d 専用」なので若干使い道が限られるが、うまくいくかどうかは別として、全素数の約25%に適用でき、同じ型の合成数にも適用できる。d < 1000 の範囲で実際に試すと、適用可能なケースの半数以上でうまくいくようだ。特に d = 61, 109, 188 のような「ペル方程式が記録的に大きい解」を持つ場合（言い換えれば、単純なブルートフォースでは非常に時間がかかる場合）、桁違いの高速化を期待できる。d = 109 のケースは、最初にやったとき、10分くらいかかった。それが一瞬。「負のペル」経由と比べても、実測で50倍くらい速く、ミリ秒（1000分の1秒）未満で検索完了する。
2020-07-10 野蛮な解法（その30） チャンピオンたち
x2 − dy2 = 1 の最小解 x1 が「新記録」の大きな数の場合、つまり、それより小さいどの d に対する最小解より大きい場合、便宜上その d をチャンピオンと呼ぶことにする。
〔例〕 d = 61 のとき x1 = 1766319049。これほど大きな最小解は d < 61 では見られないから、これはチャンピオン。次のチャンピオンは d = 109。
106 以下の約90個のチャンピオンたち（d の値）は、A033316 で見ることができる。コンピューターを使えば、もっと大きいチャンピオンを見つけることも難しくない。後で使うので、とりあえず2桁上の 108 までの全整数をチェックして、153個の新チャンピオンを得た。
1026061, 1027261, 1047589, 1053061, 1094461, 1117741, 1130581, 1167709, 1261549, 1331269, 1336141, 1404181, 1508389, 1536781, 1546141, 1627861, 1834309, 1910869, 1940821, 2065501, 2110189, 2302309, 2314909, 2375389, 2471869, 2656501, 2679301, 2747221, 2890189, 2927149, 3089629, 3220669, 3377221, 3382549, 3412861, 3414349, 3462229, 3941029, 4017589, 4207141, 4434901, 4574749, 4628101, 4635541, 4835821, 5106469, 5256469, 5429869, 5710189, 5945941, 5956549, 6105661, 6481861, 6760909, 6765301, 6810301, 7485949, 7934749, 8061421, 8356261, 8615149, 8658589, 9508669, 9549901, 9554869, 9670621, 9747061, 10675261, 10774261, 11486029, 12443029, 12765349, 13702621, 13778389, 14005429, 14011141, 14536141, 15167629, 15652261, 15725389, 15772381, 15890989, 16285669, 16561549, 16589101, 16671901, 17100301, 18355429, 18526069, 19020901, 19543549, 19572589, 20778781, 21030109, 22027021, 23404021, 23493661, 24619501, 25959781, 28024621, 30306901, 30384589, 30490069, 31757749, 34374421, 34423789, 34489261, 34534501, 35870389, 35904829, 36292741, 37448629, 38363341, 39387709, 39984421, 42579541, 42583381, 42765661, 44450701, 47710909, 50906101, 51588541, 51817789, 53255029, 53458981, 54632341, 55749541, 57651421, 59895949, 60532861, 60677269, 61118821, 61278541, 61856701, 61910941, 62982781, 65271229, 66375709, 68463781, 69118261, 72280261, 74563141, 76653229, 77148061, 78940261, 84482941, 86262829, 90561949, 91050829, 91908469, 92408941, 92919061, 99890389
d に対して、最小解 (x1, y1) がどのくらい大きくなり得るのか。この問題は今でも完全には解明されていないようだが、20世紀前半には、既に
x1 < d√d
くらいの評価が（Ivan Matveevich Vinogradov によって？）得られている。数値的に試すと、上記の（大ざっぱな）上限だと少し過大で、チャンピオンたち（上から押さえるという意味では最悪のケース）でも、もう少し小さめ。かといって、指数の √d = d1/2 を直接まけてもらおうとすると… d1/3 乗ではまるで押さえられず、d2/5 乗でも d4/9 乗でも駄目。d1/2 乗のオーダーは譲ってもらえないようだ。
専門的な深いことは分からないが、例えば d が6桁（100万）なら、d√d は100万の1000乗、既にざっと6000桁。指数をある程度「値引き」してもらっても、ペル方程式の最小解は本質的にでかいという結論は動かない。すぐ何万桁にもなる。最先端のアルゴリズムを駆使したところで、d が20桁、30桁…あるいは100桁の具体的計算は、できるとしても相当困難ではないか。
何万桁というような大きな数の計算を続けても、それら個々の数からは「でかい」という以外、あまり情報が得られないだろう。そろそろ内容的なこと（それらの背後にある世界）に目を向けていきたい。
- 妖精の森 ♌ ペル方程式の夏 (2020-12-27)
- x2 − 79y2 = 5 を満たす整数 (x, y) は存在しません。その証明は意外と難しく、しかも隠された深い意味を持っています。この種の問題を扱います。ハイライトは、2020年夏に発見されたばかりの「改良版コンラッドの不等式」。
