Julia で ペル方程式，Pell's equation

前報の無理数の連分数近似というのは，それ自体はあまりありがたいものでもない。

しかし，ペル方程式 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

Julia で連分数展開，連分数近似

連分数展開とは，たとえば x の平方根を以下のような連続する分数の形式で表現することである。

√x = a0 + 1/(a1 + 1/(a2 + 1/(a3 + 1/(a4 + 1 / (a5 + ε)))))
εはさらに 1/(a6 + 1/(a7 + 1/...)) と同じように無限に繰り返されるが，だんだんと近似に寄与する割合が小さくなるので，いつかの時点で ε = 0 とすることで近似できる。

例えば x = 31 で a0 から a8 までを使うと，以下のような連分数による近似が得られる。
√31 ≒ 5 + 1/(1 + 1/(1 + 1/(3 + 1/(5 + 1 / (3 + 1/(1 + 1/(1 + 1/(10 + 0))))))))
    = 16063 / 2885 = 5.5677642980935875

なお，全ての分数の分子が 1 であるような場合は，正則連分数とよぶ。

√31 の連分数展開は以下のように行われる。

d = 31
z = sqrt(d) # 5.5677643628300215

1. z の整数部を取る。これが a0 である
    a0 = 5
2. z の小数部の逆数（0.5677643628300215）をとり，新たな z とする
    z = 1 / 0.5677643628300215 # 1.7612940604716716
3. z の整数部を取る。これが a1 である
    a1 = 1
4. z の小数部（0.7612940604716716）の逆数をとり，新たな z とする
    z = 1 / 0.7612940604716716 # 1.3135528725660022
この後も 3.，4. の手順を繰り返す（結局は 1. 2. でもあるが）をくりかえす。

なお，計算の途中で z が 2. のときの z の小数部と等しくなることがある。
そうすると，それ以降の z の整数部は以前のz の整数部の繰り返しになることになる。
√31 の展開の場合は，5, 1, 1, 3, 5, 3, 1, 1, 10, 1, 1, 3, 5, 3, 1, 1, 10, 1, 1, ...
のようになり，1, 1, 3, 5, 3, 1, 1, 10 が循環する。

手順の中に 「z の小数部」が必要になるが，手順を繰り返していくと精度に問題が生じる可能性がある。
そこで，対処法の一つは長精度実数（Julia では BigFloat) を使うことである（問題がないこともあるし，不十分なこともある）。

この計算プログラムでは，平方根の連分数近似だけではなく，引数に与えられる数値（有理数，無理数とも）の連分数近似を行うことができる。

function regularcontinuedfractions(z, n=10) # n は a0 も含めた項数
    an = []
    for i =0:n-1
        ai = Int(floor(z))
        append!(an, ai)
        println("$i: a[$i] = $ai,  z - ai = $(z - ai)")
        z = 1 / (z - ai)
    end
    an
end

a = regularcontinuedfractions(sqrt(big(31)), 11)
#=
0: a[0] = 5,  z - ai = 0.5677643628300219221194712989185495204763933775704143039684325856035898392542376
1: a[1] = 1,  z - ai = 0.7612940604716703203532452164864249200793988962617357173280720976005983065423672
2: a[2] = 1,  z - ai = 0.3135528725660043844238942597837099040952786755140828607936865171207179678508544
3: a[3] = 3,  z - ai = 0.1892547876100073073731570996395165068254644591901381013228108618678632797513434
4: a[4] = 5,  z - ai = 0.2838821814150109610597356494592747602381966887852071519842162928017949196290188
5: a[5] = 3,  z - ai = 0.5225881209433406407064904329728498401587977925234714346561441952011966130611749
6: a[6] = 1,  z - ai = 0.9135528725660043844238942597837099040952786755140828607936865171207179679371505
7: a[7] = 1,  z - ai = 0.09462739380500365368657854981975825341273222959506905066140543093393163977229678
8: a[8] = 10,  z - ai = 0.5677643628300219221194712989185495204763933775704143039684325856035898508027295
9: a[9] = 1,  z - ai = 0.7612940604716703203532452164864249200793988962617357173280720976005982707171392
10: a[10] = 1,  z - ai = 0.3135528725660043844238942597837099040952786755140828607936865171207180296644554

11-element Vector{Any}:
  5
  1
  1
  3
  5
  3
  1
  1
 10
  1
  1
=#

もう一つの対処法は，z の小数部を使わない方法である。
for 文の中の 2 行目，floor(Int, (x + u) / v2) で整数部をとるが，小数部はこのプログラムのどこにも出てこない。
ただし，以下のプログラムは，引数 d の平方根 √d の連分数近似に限定される。

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

rcf = regularcontinuedfractions2(31)
#=
v = 1, a = 5, u = 5
v = 6.0, a = 1, u = 1.0
v = 5.0, a = 1, u = 4.0
v = 3.0, a = 3, u = 5.0
v = 2.0, a = 5, u = 5.0
v = 3.0, a = 3, u = 4.0
v = 5.0, a = 1, u = 1.0
v = 6.0, a = 1, u = 5.0
v = 1.0, a = 10, u = 5.0
v = 6.0, a = 1, u = 1.0

10-element Vector{Int64}:
  5
  1
  1
  3
  5
  3
  1
  1
 10
  1
=#

連分数近似の結果

次に，連分数展開の結果による近似値を計算する関数を定義しよう。
regularcontinuedfractions() が返す結果ベクトル a_i, i = 0, 1, ..., n-1 を与える。

function approximation(a)
    x = a[end]
    for i in reverse(a[2:end-1])
        x = i + 1/x
    end
    a[1] + 1 / x
end

rcf = regularcontinuedfractions2(31, 9) # n は a0 も含めた項数
#=
v = 1, a = 5, u = 5
v = 6.0, a = 1, u = 1.0
v = 5.0, a = 1, u = 4.0
v = 3.0, a = 3, u = 5.0
v = 2.0, a = 5, u = 5.0
v = 3.0, a = 3, u = 4.0
v = 5.0, a = 1, u = 1.0
v = 6.0, a = 1, u = 5.0
v = 1.0, a = 10, u = 5.0

9-element Vector{Int64}:
  5
  1
  1
  3
  5
  3
  1
  1
 10
=#

approximation(rcf) # 5.5677642980935875

途中の近似値を計算するための分数を全て返す関数を定義しよう。
連分数近似 regularcontinuedfractions() が返す結果ベクトル a_i, i = 0, 1, ..., n-1 を与える。
(p(i), q(i)), i = 0, 1, ... のタプルのベクトル。
i までの結果を使う近似値は p(i) / q(i) で得られる。

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

rcf = regularcontinuedfractions2(31, 9) # n は a0 も含めた項数
b = approximation2(rcf)
#=
0: 5 / 1 = 5.0
1: 6 / 1 = 6.0
2: 11 / 2 = 5.5
3: 39 / 7 = 5.571428571428571428571428571428571428571428571428571428571428571428571428571419
4: 206 / 37 = 5.567567567567567567567567567567567567567567567567567567567567567567567567567558
5: 657 / 118 = 5.567796610169491525423728813559322033898305084745762711864406779661016949152556
6: 863 / 155 = 5.567741935483870967741935483870967741935483870967741935483870967741935483870972
7: 1520 / 273 = 5.567765567765567765567765567765567765567765567765567765567765567765567765567756
8: 16063 / 2885 = 5.567764298093587521663778162911611785095320623916811091854419410745233968804128

9-element Vector{Tuple{BigInt, BigInt}}:
 (5, 1)
 (6, 1)
 (11, 2)
 (39, 7)
 (206, 37)
 (657, 118)
 (863, 155)
 (1520, 273)
 (16063, 2885)
 =#

b[7][1] / b[7][2] # 5.567741935483870967741935483870967741935483870967741935483870967741935483870972

regularcontinuedfractions() は無理数も連分数近似できる。

π の近似

a = regularcontinuedfractions(π)
approximation2(a)
big(π)
#=                    3.141592653589793238462643383279502884197169399375105820974944592307816406286198
9: 1146408 / 364913 = 
 3.14159265359140397848254241421927966391989323482583519907484797746312134673195
=#

ネイピア数

a = regularcontinuedfractions(ℯ, 50);
approximation2(a)
big(ℯ)
#=
 2.718281828459045235360287471352662497757247093699959574966967627724076630353555
9: 1457 / 536 =
 2.718283582089552238805970149253731343283582089552238805970149253731343283582103
49: 3471379623627236787997 / 1277049196033919722542 =
 2.718281828459045067661362633906433282499985585104088805075501397944665141326474
=#

有理数の連分数近似

a = regularcontinuedfractions(big"1.2345", 10);
approximation(a) # 1.2345
approximation2(a)
# 9: 2469 / 2000 = 1.234499999999999999999999999999999999999999999999999999999999999999999999999991