pellの行列で平方根近似
１．不定方程式
このワークシートはMath by Codeの一部です。
方程式には特別なタイプがあるね。
それは、解を整数や有理数に限定する不定方程式というものだ。
普通の代数方程式なら、１次方程式なら１個、２次方程式なら２個、３次方程式なら３個,…
というように、次数で解が最大で何個あるかが決まるし、解の公式がある。
しかし、不定方程式では１次でも解が無限にあったり、２次でも３個以上あったりする。
代数方程式向けの解の公式は使えなくても、整数条件に着目することで、
整式(多項式）の積が整数になることから、約数・倍数関係に着眼することが解法の鍵になるね。
さあ、不定方程式の世界に足を踏み入れてみよう。
＜一次式＞
質問：nが整数のとき、5x+3y=2nとなる整数x,yの組を求めるにはどうしますか。
5(xの式)=3(yの式)の形に変形してみよう。
5-3=2だから、2n=5n-3n。5x+3y=5n-3nとおける。3(y+n)=5(n-x)となる。
だから、y+n=-5k , n-x=-3k。たとえば、(x,y)=(n+3k, -n-5k)
正の整数に限定しないなら、整数解は無数にあるね。
練習問題：nが整数のとき、5x+3y=2となる整数x,yの一般項を求めよう。
5x-+3y=5-3とすると、3(y+1)=5(1-x)となるから、y+1=-5k, 1-x=-3k。
たとえば、(x,y)=(1+3k, -1-5k)
＜２次式＞
質問：nが整数のとき、xy-2x-y=4となる自然数x,yの組を求めるにはどうしますか。
因数分解しよう。
(x-1)(y-2)=6 ⇒(x-1,y-2)=(1,6),(2,3),(3,2),(6,1) ⇒(x,y)=(2,8),(3,5),(4,4),(7,3)
２次式なのに解は４組あった。
練習問題：nが整数のとき、6x2+xy-2y2-5x+6y=20となる自然数x,yの組をすべて求めよう。
先に2次式だけ整理する。6x2+xy-2y2=(2x -y)(3x +2y)
(2x -y+p)(3x +2y+q)=6x2+xy-2y2+(3p+2q)x+(2p-q)y+pq と係数比較して、p=1,q=-4を得るからpq=-4
(2x -y+1)(3x +2y-4)=20-4=16
⇒(2x -y+1,3x +2y-4)=(1,16),(2,8),(4,4),(8,2),(16,1)
⇒(4x -2y, 3x +2y)=(0,20),(2,12),(6,8),(14,6),(30,5)
⇒(7x, 3x +2y)=(20,20),(14,12),(14,8),(20,6),(35,5)
⇒(x, y)=,(2,3),(2,1),,
＜３次式＞
質問：nが整数のとき、x3+y3=91となる整数x,yの組を求めるにはどうしますか。
因数分解しよう。
x+y=p, xy=qとおくと、x3+y3=(x+y)(x2-xy+y2)=(x+y)((x+y)2-3xy)=p(p2-3q)=91=7*13。
p=1のとき、q=(1-91)/3=-30より、x+y=1,xy=-30=-5*6から、x,y=(6,-5),(-5,6)
p=7のとき、q=(49-13)/3=12より、x+y=7,xy=12=3*4から、x,y=(3,4)(4,3)
(x,y)=(6,-5),(-5,6),(3,4)(4,3)
３次方程式なのに解が４組あった。
練習問題：nが整数のとき、x3-y3=217となる整数x,yの組を求めよう。
x-y=p, -xy=qとおくと、x3-y3=(x-y)(x2+xy+y2)=p(p2-3q)=217=7*31。
p=1のとき、q=(1-217)/3=-72より、x-y=1,xy=72=9*8から、x,y=(9,8),(-8,-9)
p=7のとき、q=(49-31)/3=6より、x-y=7,xy=-6=-6*1から、x,y=(1,-6),(6,-1)
p=31のとき、q=(31^2-7)/3=318より、x+(-y)=31,x(-y)=318はt^2-31t+318=0の解で、D<0で実数解なし。
p=217のとき、q=(217^2-1)/3=15696より、x+(-y)=217,x(-y)=15696の場合も上記のように実数解なし。
(x,y)=(9,8),(-8,-9),(1,-6),(6,-1)
曲線と格子点を目で確認しよう
２．ペルの方程式
昔から有名な不定方程式があるね。
それは、ペルの方程式だ。
解が無数にある。
x2 -Dy2 = 1（Dが平方数でない。ｘ，ｙは整数)
これがペルの不定方程式で、オイラー・ペルの方程式と呼んだりもする。
ペルの方程式は最小解が見つかれば、
それから芋づる式に他の解を求めることができることがわかっている。
しかも、これの解がルートDの近似値計算にとても役立つ。
＜ペル方程式の簡単な例＞
では、簡単な例から見ていこう。
・D=2のとき、x2-２y2=1。
(x,y)=(1,0),(3,2)はすぐに見つかるでしょう。
そこで、(x1,y1)=(3,2)としよう。
(x+√2x)(x-√2y)=1と因数分解して、左辺と右辺をｎ乗しても
(x+√2x)n(x-√2y)n=1となり、(x+√2y)nを展開整理すれば、u+√2vの形にかけるね。
そこで、(x1+√2y1)n=xn+√2ynと決めることにすると、xn,ynはx2-２y2=1の解になることは、数学的帰納法で
証明できるでしょう。
次に、漸化式を作ってみよう。
xn+√2yn=(x1+√2y1)(xn-1+√2yn-1)=(3+2√2)(xn-1+√2yn-1)=(3xn-1+4yn-1)+√2(2xn-1+3yn-1)となる。
これって、変数だけ取り出すと、
xn=3xn-1+4yn-1
yn=2xn-1+3yn-1となる。
なんと、線形変換の行列計算で次々と解が出せるということだ。
質問：行列計算で求めるにはどんなコードをかきますか？
たとえばjuliaで行列の積を実行する関数と確認する関数を作ってみる。
実行する関数は、０番目の解xy[1;0]からスタートして、行列の積を返すfp2(xy)。
確認する関数は、xyの値をペル方程式の左辺に代入して、たぶん１を返すpell2(xy)。
この実行と確認を10回くり返してみよう。
#[IN]julia
#===========================================
using LinearAlgebra
A = [ 3 4; 2 3 ]
xy = [ 1; 0 ]
function fp2(xy)
A*xy
end
function pell2(xy)
(x,y)=xy
return x^2-2*y^2
end
res=[]
values=[ ]
for i in 1:10
xy=fp2(xy)
push!(res,xy)
push!(values,pell2(xy))
end
println(res)
println(values)
#===========================================
[OUT]
Any[[3, 2], [17, 12], [99, 70], [577, 408], [3363, 2378], [19601, 13860], [114243, 80782], [665857, 470832], [3880899, 2744210], [22619537, 15994428]]
Any[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
10回行列計算をしてペル方程式にいれると、毎回ちゃんと=１になってるね。
さて、次はxn÷ynの値を計算してみよう。
#[IN]julia
#===========================================
xy = [ 1; 0 ]
for i in 1:10
xy=fp2(xy)
println(xy[1]/xy[2])
end
#===========================================
1.5
1.4166666666666667
1.4142857142857144
1.4142156862745099
1.4142136248948696
1.4142135642135643
1.4142135624272734
1.4142135623746899
1.414213562373142
1.4142135623730965
さあ、もう気づきましたね。
これは√2の近似値になっています。
当然と言えば当然です。
x2-２y2=1から、x2=２y2+1。これを両辺をy2で割ると、(x/y)2=2+(1/y)2→2
nが大きくなると1/yは０に近づくので、(x/y)の値が√2に近づくからです。
ペルの方程式の解の漸化式
３．一般化と高速化
＜一般化＞
２をDに戻すとどうなるでしょうか？
D=2のとき、x2-２y2=1。最小解を(x1,y1)=(3,2)として漸化式
xn+√2yn=(x1+√2y1)(xn-1+√2yn-1)=(3+2√2)(xn-1+√2yn-1)=(3xn-1+4yn-1)+√2(2xn-1+3yn-1)を得ました。
もし、x2-Dy2=1。最小解を(x1,y1)=(p,q)として漸化式
xn+√Dyn=(x1+√Dy1)(xn-1+√Dyn-1)=(p+q√D)(xn-1+√Dyn-1)=(pxn-1+qDyn-1)+√D(qxn-1+pyn-1)を得ます。
最小解(p,q)がわかれば、
xn=pxn-1+qDyn-1
yn=qxn-1+pyn-1となる。
質問：√3、√5、√7、√11の最小解は何でしょう？
また、近似値を効率的に出すプログラムはどうやればよいでしょうか。
(x0,y0)=(1,0)は共通だね。
x2-3y2=1 最小解(x1,y1)=(2,1)から、行列は[2 3; 1 2]
x2-5y2=1 最小解(x1,y1)=(9,4)から、行列は[9 20;4 9]
x2-7y2=1 最小解(x1,y1)=(8,3)から、行列は[8 21; 3 8]
x2-11y2=1 最小解(x1,y1)=(10,3)から、行列は[10 33; 3 10]
√2,√3,√5,√7,√11の近似値を行列を14回実行して求めてみる。
Dと最小(p,q)の組、(D,p,q)をデータをして渡すと、14回かけ算したあとの解を返して表示する関数が
apx(x)だ。ついでに、標準で入っている平方根の近似値を並べて表示してみよう。
#[IN]julia
#==================================================
using LinearAlgebra
data=[(2,3,2),(3,2,1),(5,9,4),(7,8,3),(11,10,3)] #(D,p,q)
function apx(x)
(D,p,q)=x
A = [ p q*D; q p ] #漸化式の行列
xy = [ 1; 0 ] #初期値
for i in 1:14
xy = A*xy
end
println("近似値 =",BigFloat(xy[1]/xy[2]))
println("sqrt $D =",sqrt(D))
end
for i in data
apx(i)
end
#==================================================
[OUT]
近似値 =1.4142135623730951454746218587388284504413604736328125
sqrt 2 =1.4142135623730951
近似値 =1.73205080756887763726581397349946200847625732421875
sqrt 3 =1.7320508075688772
近似値 =2.236067977499789360962267892318777740001678466796875
sqrt 5 =2.23606797749979
近似値 =2.64575131106459071617109657381661236286163330078125
sqrt 7 =2.6457513110645907
近似値 =3.316624790355400254071582821779884397983551025390625
sqrt 11=3.3166247903554
14回行列計算を回しただけで、標準で入っている近似値に遜色ない値になるね。
すばらしい！
＜高速化＞
n番目のｘｙを求める定義式(x1+√Dy1)n=xn+√Dyn
を変形することで、x2n=xn2+Dyn2, y2n= 2xnyn
が得られる。（途中は相当省略）
すると、第１式x2n=xn2+Dyn2 とxn2 - Dyn2=1から、
x2n=2xn2-1が得られるね。
だから、ペルの方程式にはDの値に関係なく、
x2n=2xn2-1,
y2n= 2xnyn
というとても役立つ性質がある。これで倍々で先に進めるね。
質問：この高速計算を使って、√２の近似値を求めるためのコードをかいてみよう。
juliaでやる。
#[IN]julia
#==================================================
using LinearAlgebra
xy = [ 3; 2 ]
function fp2(xy)
(x,y)=xy
y=2x*y
x=2x^2-1
return (x,y)
end
function pell2(xy)
(x,y)=xy
return x^2-2*y^2
end
res=[]
values=[ ]
for i in 1:4
xy=fp2(xy)
push!(res,xy)
push!(values,pell2(xy))
end
println(res)
println(values)
#==================================================
[OUT]
Any[(17, 12), (577, 408), (665857, 470832), (886731088897, 627013566048)]
Any[1, 1, 1, 1]
４．平方三角数
おまけになりますが、D=2のペルの方程式の使い道がもう一つあります。
それは、「平方三角数」です。
平方数はn２が一般項、
三角数はm(m+1)/2が一般項でした。
これが両立する数はn2=m(m+1)/2ですね。それを平方三角数と呼びましょう。
D=2のペルの方程式の解の数列から平方三角数が作れるのです。
なぜか、x2-2y2=1の解ｘ，ｙに対してm=(x-1)/2, n=y/2を求めます。
すると、ｘ＝2m+1, y =2nなので、ペルの方程式に代入できます。
(2m+1)2-1=8n2
4m2+4m=8n2
m(m+1)/2=n2
質問：ペルの方程式の解から平方三角数の数列を表示するにはどうしますか。
#[IN]Julia
#======================================================
#平方三角数
using LinearAlgebra
A = [ 3 4; 2 3 ]
xy = [ 3; 2 ]
function fp2(xy)
(x,y)=xy
m = (x-1)÷2; n = y÷2;
triSqr=BigInt(n)^2
println("(x,y)=",xy,",(m,n)=(",m,",",n,"),平方三角数=",triSqr)
return A*xy
end
for i in 1:20
xy=fp2(xy)
end
#=====================================================
(x,y)=[3, 2],(m,n)=(1,1),平方三角数=1
(x,y)=[17, 12],(m,n)=(8,6),平方三角数=36
(x,y)=[99, 70],(m,n)=(49,35),平方三角数=1225
(x,y)=[577, 408],(m,n)=(288,204),平方三角数=41616
(x,y)=[3363, 2378],(m,n)=(1681,1189),平方三角数=1413721
(x,y)=[19601, 13860],(m,n)=(9800,6930),平方三角数=48024900
(x,y)=[114243, 80782],(m,n)=(57121,40391),平方三角数=1631432881
(x,y)=[665857, 470832],(m,n)=(332928,235416),平方三角数=55420693056
(x,y)=[3880899, 2744210],(m,n)=(1940449,1372105),平方三角数=1882672131025
(x,y)=[22619537, 15994428],(m,n)=(11309768,7997214),平方三角数=63955431761796
(x,y)=[131836323, 93222358],(m,n)=(65918161,46611179),平方三角数=2172602007770041
(x,y)=[768398401, 543339720],(m,n)=(384199200,271669860),平方三角数=73804512832419600
(x,y)=[4478554083, 3166815962],(m,n)=(2239277041,1583407981),平方三角数=2507180834294496361
(x,y)=[26102926097, 18457556052],(m,n)=(13051463048,9228778026),平方三角数=85170343853180456676
(x,y)=[152139002499, 107578520350],(m,n)=(76069501249,53789260175),平方三角数=2893284510173841030625
(x,y)=[886731088897, 627013566048],(m,n)=(443365544448,313506783024),平方三角数=98286503002057414584576
(x,y)=[5168247530883, 3654502875938],(m,n)=(2584123765441,1827251437969),平方三角数=3338847817559778254844961
(x,y)=[30122754096401, 21300003689580],(m,n)=(15061377048200,10650001844790),平方三角数=113422539294030403250144100
(x,y)=[175568277047523, 124145519261542],(m,n)=(87784138523761,62072759630771),平方三角数=3853027488179473932250054441
(x,y)=[1023286908188737, 723573111879672],(m,n)=(511643454094368,361786555939836),平方三角数=130889512058808083293251706896