1．brute force method による解

brute force method による探査

まずは感触をつかむために、1000以下の x, y の組について n = (x²-1)/y², (x²+1)/y² を計算し、
n が 100以下の自然数、かつ 2乗因子を含まないケースを探す。
任意の n について解を持つことが保証されている（参考文献）。

x² - ny² = -1 の解がわかれば、+1 の場合の解を構成することができる。
左辺を実数の範囲で因数分解すると、

(x - sqrt(n)y)(x + sqrt(n)y) = -1

両辺を 2乗すると、

(x - sqrt(n)y)²(x + sqrt(n)y)² = 1

ここで、

(x - sqrt(n)y)² = x² + ny² - 2xy * sqrt(n)
(x + sqrt(n)y)² = x² + ny² + 2xy * sqrt(n)

なので、

X = x² + ny²
Y = 2xy

とおくと、

1 = (x - sqrt(n)y)²(x + sqrt(n)y)²
  = {(x² + ny²) - 2xy * sqrt(n)}{(x² + ny²) + 2xy * sqrt(n)}
  = (X - Ysqrt(n))(X + Ysqrt(n))
  = X² - (Ysqrt(n))²
  = X² - nY²

となり、確かに +1 の場合の解となっている。
よって探査は、1か - 1のどちらかの解が見つかった時点で終了してよい。

プログラム

nについての解が見つかったかどうかを判定するために、配列 M%(100) をとり、
0で初期化、解が見つかったら 1をセットする。
x, y を for 〜 next で 1 から 1000 ぐらいまでまわして、

n = (x² +1) / y² が整数であり、かつ 100 以下の場合
n = (x² -1) / y² が整数であり、かつ 100 以下の場合

に、n, x, y, x² - ny² を 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% が２乗因子を持つとき０を返す。
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 でループを回した方が、回数が少なくてすむ。
x² = ny² ±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以下の解についてはすべて求められた。