前報の無理数の連分数近似というのは，それ自体はあまりありがたいものでもない。
しかし，ペル方程式 Pell's equation を解くときに力を発揮する。
ペル方程式とは，
n を平方数ではない自然数として、未知整数 x, y について
x^2 − n*y^2 = 1
の形のディオファントス方程式(整係数多変数高次不定方程式)である。
1 つの解 (x, y) が得られたとき，
x_k + y_k √n = (x + y √n)^k
は，全てペル方程式の解になる。
最小解は、√n の連分数展開の結果を利用する方法がある。
周期が m のとき，a_0, a_1, ..., a_m とすると，a_m までの近似分数を P/Q としたとき，P, Q がペル方程式の解である。
ただし，m が奇数の場合は得られるのは x^2 − n*y^2 = −1 の解なので，それを 2 乗する。または，二倍の周期，すなわち a_2m までの近似分数 P/Q の P, Q が解である。
function regularcontinuedfractions2(d, n=10) # n は a0 も含めた項数
x = sqrt(d)
v = 1
u = a = floor(Int, x)
an = [a]
println("v = $v, a = $a, u = $u")
for i = 1:n-1
v2 = (d - u^2) / v
a2 = floor(Int, (x + u) / v2)
u2 = a2 * v2 - u
v, a, u = v2, a2, u2
println("v = $v, a = $a, u = $u")
append!(an, a)
end
an
end
function approximation2(a)
p0 = BigInt(a[1])
q0 = BigInt(1)
println("0: $p0 / $q0 = $(p0 / q0)")
p1 = BigInt(a[1]*a[2] + 1)
q1 = BigInt(a[2])
println("1: $p1 / $q1 = $(p1 / q1)")
result = [(p0, q0), (p1, q1)]
for i = 3:length(a)
p2 = BigInt(a[i]) * p1 + p0
q2 = BigInt(a[i]) * q1 + q0
println("$(i-1): $p2 / $q2 = $(p2/q2)")
append!(result, [(p2, q2)])
p0, q0, p1, q1 = p1, q1, p2, q2
end
result
end
a = regularcontinuedfractions2(61, 25);
for i = 1:length(a)
println("a[$(i-1)] = $(a[(i)])")
end
#=
a[0] = 7
a[1] = 1
a[2] = 4
a[3] = 3
a[4] = 1
a[5] = 2
a[6] = 2
a[7] = 1
a[8] = 3
a[9] = 4
a[10] = 1
a[11] = 14
a[12] = 1
a[13] = 4
a[14] = 3
a[15] = 1
a[16] = 2
a[17] = 2
a[18] = 1
a[19] = 3
a[20] = 4
a[21] = 1
a[22] = 14
a[23] = 1
a[24] = 4
=#
周期は 11（1,4,3,1,2,2,1,3,4,1,14） で，奇数なので，a[22] までの解を用いる。
b = approximation2(a);
for i = 1:25
println("$i, $(b[i][1]), $(b[i][2]), $((b[i][1])^2 - (b[i][2])^2 * 61)")
end
#=
1, 7, 1, -12
2, 8, 1, 3
3, 39, 5, -4
4, 125, 16, 9
5, 164, 21, -5
6, 453, 58, 5
7, 1070, 137, -9
8, 1523, 195, 4
9, 5639, 722, -3
10, 24079, 3083, 12
11, 29718, 3805, -1
12, 440131, 56353, 12
13, 469849, 60158, -3
14, 2319527, 296985, 4
15, 7428430, 951113, -9
16, 9747957, 1248098, 5
17, 26924344, 3447309, -5
18, 63596645, 8142716, 9
19, 90520989, 11590025, -4
20, 335159612, 42912791, 3
21, 1431159437, 183241189, -12
22, 1766319049, 226153980, 1
23, 26159626123, 3349396909, -12
24, 27925945172, 3575550889, 3
25, 137863406811, 17651600465, -4
=#
a[22] までの連分数近似では，(P, Q) = (x, y) = (1766319049, 226153980)
1766319049^2 - 61 * 226153980^2 # 1
「x_k + y_k √n = (x + y √n)^k」というのは，単純な数値計算をするのではない。
(x + √n * y)^2 # 1.2479531931441058e19 という1つの数になる
(x + √n * y)^2 を数式として展開すると
x^2 + n*y^2 + 2xy√n
で，新たな解 x は x^2 + n*y^2，y は 2xy である。
つまり，√n を含まない項（の和）が x，√n を含む項の係数が y ということである。
x = 1766319049^2 + 61 * 226153980^2 # 6239765965720528801
y = 2 * 1766319049 * 226153980 # 798920165762330040
つまり，(x, y) = (6239765965720528801, 798920165762330040) である。
x^2 - 61 * y^2 # 1
