まずは感触をつかむために、1000以下の x, y の組について
n = (x2-1)/y2, (x2+1)/y2 を計算し、
n が 100以下の自然数、かつ 2乗因子を含まないケースを探す。
任意の n について解を持つことが保証されている(参考文献)。
x2 - ny2 = -1 の解がわかれば、+1 の場合の解を構成することができる。
左辺を実数の範囲で因数分解すると、
(x - sqrt(n)y)(x + sqrt(n)y) = -1
両辺を 2乗すると、
(x - sqrt(n)y)2(x + sqrt(n)y)2 = 1
ここで、
(x - sqrt(n)y)2 = x2 + ny2 - 2xy * sqrt(n)
(x + sqrt(n)y)2 = x2 + ny2 + 2xy * sqrt(n)
なので、
X = x2 + ny2
Y = 2xy
とおくと、
1 = (x - sqrt(n)y)2(x + sqrt(n)y)2
= {(x2 + ny2) - 2xy * sqrt(n)}{(x2 + ny2) + 2xy * sqrt(n)}
= (X - Ysqrt(n))(X + Ysqrt(n))
= X2 - (Ysqrt(n))2
= X2 - nY2
となり、確かに +1 の場合の解となっている。
よって探査は、1か - 1のどちらかの解が見つかった時点で終了してよい。
nについての解が見つかったかどうかを判定するために、配列 M%(100) をとり、
0で初期化、解が見つかったら 1をセットする。
x, y を for 〜 next で 1 から 1000 ぐらいまでまわして、
n = (x2 +1) / y2 が整数であり、かつ 100 以下の場合
n = (x2 -1) / y2 が整数であり、かつ 100 以下の場合
に、n, x, y, x2 - ny2 を print し、M%(n) = 1とする。
プログラムは以下のとおり。
10 ' pell1 20 dim M%(100) 30 for I%=2 to 100 40 if moeb(I%)=0 then M%(I%)=1 else M%(I%)=0 50 next I% 60 ' 70 for X%=1 to 1000:X2=X%^2 80 for Y%=1 to 1000:Y2=Y%^2 90 if (X2+1)@Y2=0 then W=(X2+1)\Y2 else 120 100 if or{W<1,W>=100} then 120 110 if M%(W)=0 then print W;X%;Y%;X2-W*Y2:M%(W)=1 120 if (X2-1)@Y2=0 then W=(X2-1)\Y2 else 150 130 if or{W<1,W>=100} then 150 140 if M%(W)=0 then print W;X%;Y%;X2-W*Y2:M%(W)=1 150 next Y% 160 next X% 170 ' 180 for I%=2 to 100 190 if M%(I%)=0 then print I%; 200 next I%
40行の moeb(I%) は、メビウス関数といい、I% が2乗因子を持つとき0を返す。
1000×1000 回のループが終わったら、180〜200行で、解が得られなかった nの値を返す。
実行結果は以下のとおり。
n | x | y | x2 - ny2 |
---|---|---|---|
2 | 1 | 1 | -1 |
3 | 2 | 1 | 1 |
5 | 2 | 1 | -1 |
6 | 5 | 2 | 1 |
7 | 8 | 3 | 1 |
10 | 3 | 1 | -1 |
11 | 10 | 3 | 1 |
13 | 18 | 5 | -1 |
14 | 15 | 4 | 1 |
15 | 4 | 1 | 1 |
17 | 4 | 1 | -1 |
19 | 170 | 39 | 1 |
21 | 55 | 12 | 1 |
22 | 197 | 42 | 1 |
23 | 24 | 5 | 1 |
26 | 5 | 1 | -1 |
29 | 70 | 13 | -1 |
30 | 11 | 2 | 1 |
33 | 23 | 4 | 1 |
34 | 35 | 6 | 1 |
35 | 6 | 1 | 1 |
37 | 6 | 1 | -1 |
38 | 37 | 6 | 1 |
39 | 25 | 4 | 1 |
41 | 32 | 5 | -1 |
42 | 13 | 2 | 1 |
47 | 48 | 7 | 1 |
51 | 50 | 7 | 1 |
53 | 182 | 25 | -1 |
55 | 89 | 12 | 1 |
57 | 151 | 20 | 1 |
58 | 99 | 13 | -1 |
59 | 530 | 69 | 1 |
62 | 63 | 8 | 1 |
65 | 8 | 1 | -1 |
66 | 65 | 8 | 1 |
70 | 251 | 30 | 1 |
74 | 43 | 5 | -1 |
77 | 351 | 40 | 1 |
78 | 53 | 6 | 1 |
79 | 80 | 9 | 1 |
82 | 9 | 1 | -1 |
83 | 82 | 9 | 1 |
85 | 378 | 41 | -1 |
87 | 28 | 3 | 1 |
89 | 500 | 53 | -1 |
95 | 39 | 4 | 1 |
見つからなかったのは、
31, 43, 46, 61, 67, 69, 71, 73, 86, 91, 93, 94, 97
の13個。
次に、見つからなかった n について x, y の値を大きくして探してみる。
とりあえず 10000 ぐらいまでとする。n は毎回入力するようにする。
上の解からわかるように x と y では y の方が小さいので、y でループを回した方が、回数が少なくてすむ。
x2 = ny2 ±1 のルートを取り、整数値となる場合を探す。
プログラムは以下のとおり。
10 ' pell2 20 input N:if moeb(N)=0 then 20 30 ' 40 for Y=2 to 10000:Y2=Y^2 50 W=N*Y2-1:X=isqrt(W):if res=0 then print N;X;Y;W-N*Y2:end 60 W=N*Y2+1:X=isqrt(W):if res=0 then print N;X;Y;W-N*Y2:end 70 next Y
実行結果は以下のとおり。
n | x | y | x2 - ny2 |
---|---|---|---|
31 | 1520 | 273 | 1 |
43 | 3482 | 531 | 1 |
46 | 24335 | 3588 | 1 |
61 | 29718 | 3805 | -1 |
67 | 48842 | 5967 | 1 |
69 | 7775 | 936 | 1 |
71 | 3480 | 413 | 1 |
73 | 1068 | 125 | -1 |
86 | 10405 | 1122 | 1 |
91 | 1574 | 165 | 1 |
93 | 12151 | 1260 | 1 |
97 | 5604 | 569 | -1 |
残るは n=94 の場合だけである。
上記のプログラムで、yの上限値を 10^9 ぐらいにしておく。
実行すると、割と簡単にケリがつく。
94 : 2143295 , 221064 : 1
221064回のループの後に、解が見つかった。
これで 100以下の解についてはすべて求められた。
前 | この章の目次 | 次 |
---|