[b]このワークシートは[url=https://www.geogebra.org/m/twxxx3yq]Math by Code[/url]の一部です。[br][/b][br]方程式には特別なタイプがあるね。[br]それは、解を整数や有理数に限定する不定方程式というものだ。[br]普通の代数方程式なら、１次方程式なら１個、２次方程式なら２個、３次方程式なら３個,…[br]というように、次数で解が最大で何個あるかが決まるし、解の公式がある。[br]しかし、不定方程式では１次でも解が無限にあったり、２次でも３個以上あったりする。[br][br]代数方程式向けの解の公式は使えなくても、整数条件に着目することで、[br]整式(多項式）の積が整数になることから、約数・倍数関係に着眼することが解法の鍵になるね。[br][br]さあ、不定方程式の世界に足を踏み入れてみよう。[br][br][b][size=150]＜一次式＞[/size][/b][br][color=#9900ff][u][b][size=150]質問：nが整数のとき、5x+3y=2nとなる整数x,yの組を求めるにはどうしますか。[br][/size][/b][/u][/color][br]5(xの式)=3(yの式)の形に変形してみよう。[br]5-3=2だから、2n=5n-3n。5x+3y=5n-3nとおける。3(y+n)=5(n-x)となる。[br]だから、y+n=-5k , n-x=-3k。たとえば、(x,y)=(n+3k, -n-5k)[br]正の整数に限定しないなら、整数解は無数にあるね。[br][br]練習問題：nが整数のとき、5x+3y=2となる整数x,yの一般項を求めよう。[br]5x-+3y=5-3とすると、3(y+1)=5(1-x)となるから、y+1=-5k, 1-x=-3k。[br]たとえば、(x,y)=(1+3k, -1-5k)[br][br][b][size=150]＜２次式＞[br][/size][/b][u][b][size=150][color=#9900ff]質問：nが整数のとき、xy-2x-y=4となる自然数x,yの組を求めるにはどうしますか。[br][/color][/size][/b][/u][br]因数分解しよう。[br](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) [br]２次式なのに解は４組あった。[br][br]練習問題：nが整数のとき、6x[sup]2[/sup]+xy-2y[sup]2[/sup]-5x+6y=20となる自然数x,yの組をすべて求めよう。[br]先に2次式だけ整理する。6x[sup]2[/sup]+xy-2y[sup]2[/sup]=(2x -y)(3x +2y)[br](2x -y+p)(3x +2y+q)=6x[sup]2[/sup]+xy-2y[sup]2[/sup]+(3p+2q)x+(2p-q)y+pq と係数比較して、p=1,q=-4を得るからpq=-4[br](2x -y+1)(3x +2y-4)=20-4=16[br]⇒(2x -y+1,3x +2y-4)=(1,16),(2,8),(4,4),(8,2),(16,1) [br]⇒(4x -2y, 3x +2y)=(0,20),(2,12),(6,8),(14,6),(30,5) [br][size=150]⇒(7x, 3x +2y)=(20,20),(14,12),(14,8),(20,6),(35,5) [br][b]⇒(x, y)=,(2,3),(2,1),, [br][br]＜３次式＞[br][/b][/size][color=#9900ff][u][b][size=150]質問：nが整数のとき、x[sup]3[/sup]+y[sup]3[/sup]=91となる整数x,yの組を求めるにはどうしますか。[br][/size][/b][/u][/color][br]因数分解しよう。[br]x+y=p, xy=qとおくと、x[sup]3[/sup]+y[sup]3[/sup]=(x+y)(x[sup]2[/sup]-xy+y[sup]2[/sup])=(x+y)((x+y)[sup]2[/sup]-3xy)=p(p[sup]2[/sup]-3q)=91=7*13。[br]p=1のとき、q=(1-91)/3=-30より、x+y=1,xy=-30=-5*6から、x,y=(6,-5),(-5,6)[br]p=7のとき、q=(49-13)/3=12より、x+y=7,xy=12=3*4から、x,y=(3,4)(4,3)[br](x,y)=(6,-5),(-5,6),(3,4)(4,3)[br]３次方程式なのに解が４組あった。[br][br][size=150][size=100]練習問題：nが整数のとき、x[sup]3[/sup]-y[sup]3[/sup]=217となる整数x,yの組を求めよう。[/size][/size][br]x-y=p, -xy=qとおくと、x[sup]3[/sup]-y[sup]3[/sup]=(x-y)(x[sup]2[/sup]+xy+y[sup]2[/sup])=p(p[sup]2[/sup]-3q)=217=7*31。[br]p=1のとき、q=(1-217)/3=-72より、x-y=1,xy=72=9*8から、x,y=(9,8),(-8,-9)[br]p=7のとき、q=(49-31)/3=6より、x-y=7,xy=-6=-6*1から、x,y=(1,-6),(6,-1)[br]p=31のとき、q=(31^2-7)/3=318より、x+(-y)=31,x(-y)=318はt^2-31t+318=0の解で、D<0で実数解なし。[br]p=217のとき、q=(217^2-1)/3=15696より、x+(-y)=217,x(-y)=15696の場合も上記のように実数解なし。[br](x,y)=(9,8),(-8,-9),(1,-6),(6,-1)[br]

昔から有名な不定方程式があるね。[br][br]それは、ペルの方程式だ。[br]解が無数にある。[br][br][size=200]x[sup]2 [/sup]-Dy[sup]2 [/sup]= 1（Dが平方数でない。ｘ，ｙは整数)[br][/size][br]これが[b][color=#0000ff]ペルの不定方程式[/color][/b]で、オイラー・ペルの方程式と呼んだりもする。[br][br]ペルの方程式は最小解が見つかれば、[br]それから芋づる式に他の解を求めることができることがわかっている。[br]しかも、これの解がルートDの近似値計算にとても役立つ。[br][br][b][size=150]＜ペル方程式の簡単な例＞[br][/size][/b]では、簡単な例から見ていこう。[br]・D=2のとき、x[sup]2[/sup]-２y[sup]2[/sup]=1。[br](x,y)=(1,0),(3,2)はすぐに見つかるでしょう。[br]そこで、(x1,y1)=(3,2)としよう。[br](x+√2x)(x-√2y)=1と因数分解して、左辺と右辺をｎ乗しても[br](x+√2x)[sup]n[/sup](x-√2y)[sup]n[/sup]=1となり、(x+√2y)nを展開整理すれば、u+√2vの形にかけるね。[br]そこで、(x[sub]1[/sub]+√2y[sub]1[/sub])[sup]n[/sup]=x[sub]n[/sub]+√2y[sub]n[/sub]と決めることにすると、x[sub]n[/sub],y[sub]n[/sub]はx[sup]2[/sup]-２y[sup]2[/sup]=1の解になることは、数学的帰納法で[br]証明できるでしょう。[br]次に、漸化式を作ってみよう。[br]x[sub]n[/sub]+√2y[sub]n[/sub]=(x[sub]1[/sub]+√2y[sub]1[/sub])(x[sub]n-1[/sub]+√2y[sub]n-1[/sub])=(3+2√2)(x[sub]n-1[/sub]+√2y[sub]n-1[/sub])=([b]3[/b]x[sub]n-1[/sub]+[b]4[/b]y[sub]n-1[/sub])+√2([b]2[/b]x[sub]n-1[/sub]+[b]3[/b]y[sub]n-1[/sub])となる。[br]これって、変数だけ取り出すと、[br][color=#9900ff][br][size=150]x[sub]n[/sub]=[b]3[/b]x[sub]n-1[/sub]+[b]4[/b]y[sub]n-1[/sub][br]y[sub]n[/sub]=[b]2[/b]x[sub]n-1[/sub]+[b]3[/b]y[sub]n-1[/sub]となる。[br][/size][/color][br]なんと、線形変換の行列計算で次々と解が出せるということだ。[br][br][color=#9900ff][u][b][size=150]質問：行列計算で求めるにはどんなコードをかきますか？[br][/size][/b][/u][/color][br]たとえばjuliaで行列の積を実行する関数と確認する関数を作ってみる。[br]実行する関数は、０番目の解xy[1;0]からスタートして、行列の積を返すfp2(xy)。[br]確認する関数は、xyの値をペル方程式の左辺に代入して、たぶん１を返すpell2(xy)。[br]この実行と確認を10回くり返してみよう。[br][br][color=#0000ff][b]#[IN]julia[br]#===========================================[br]using LinearAlgebra[br]A = [ 3 4; 2 3 ][br]xy = [ 1; 0 ][br]function fp2(xy)[br] A*xy[br]end[br]function pell2(xy)[br] (x,y)=xy[br] return x^2-2*y^2[br]end[br]res=[][br]values=[ ][br]for i in 1:10[br] xy=fp2(xy)[br] push!(res,xy)[br] push!(values,pell2(xy))[br]end[br]println(res)[br]println(values)[br][/b][/color][color=#0000ff][b]#===========================================[br][/b][/color][OUT][br]Any[[3, 2], [17, 12], [99, 70], [577, 408], [3363, 2378], [19601, 13860], [114243, 80782], [665857, 470832], [3880899, 2744210], [22619537, 15994428]][br]Any[1, 1, 1, 1, 1, 1, 1, 1, 1, 1][br][br]10回行列計算をしてペル方程式にいれると、毎回ちゃんと=１になってるね。[br][br]さて、次はx[sub]n[/sub]÷y[sub]n[/sub]の値を計算してみよう。[br][br][color=#0000ff][b]#[IN]julia[br]#===========================================[br][/b]xy = [ 1; 0 ][br]for i in 1:10[br] xy=fp2(xy)[br] println(xy[1]/xy[2])[br]end[br][b]#===========================================[br][/b][/color]1.5[br]1.4166666666666667[br]1.4142857142857144[br]1.4142156862745099[br]1.4142136248948696[br]1.4142135642135643[br]1.4142135624272734[br]1.4142135623746899[br]1.414213562373142[br]1.4142135623730965[br][br]さあ、もう気づきましたね。[br]これは[b][color=#9900ff]√2の近似値[/color][/b]になっています。[br]当然と言えば当然です。[br][br]x[sup]2[/sup]-２y[sup]2[/sup]=1から、x[sup]2[/sup]=２y[sup]2[/sup]+1。これを両辺をy2で割ると、(x/y)[sup]2[/sup]=2+(1/y)[sup]2[/sup]→2[br][br]nが大きくなると1/yは０に近づくので、(x/y)の値が√2に近づくからです。[br]

[b][size=150]＜一般化＞[br][/size][/b][br]２をDに戻すとどうなるでしょうか？[br][br]D=2のとき、x[sup]2[/sup]-２y[sup]2[/sup]=1。最小解を(x1,y1)=(3,2)として漸化式[br]x[sub]n[/sub]+√2y[sub]n[/sub]=(x[sub]1[/sub]+√2y[sub]1[/sub])(x[sub]n-1[/sub]+√2y[sub]n-1[/sub])=(3+2√2)(x[sub]n-1[/sub]+√2y[sub]n-1[/sub])=(3x[sub]n-1[/sub]+4y[sub]n-1[/sub])+√2(2x[sub]n-1[/sub]+3y[sub]n-1[/sub])を得ました。[br][br]もし、x[sup]2[/sup]-Dy[sup]2[/sup]=1。最小解を(x1,y1)=(p,q)として漸化式[br]x[sub]n[/sub]+√Dy[sub]n[/sub]=(x[sub]1[/sub]+√Dy[sub]1[/sub])(x[sub]n-1[/sub]+√Dy[sub]n-1[/sub])=(p+q√D)(x[sub]n-1[/sub]+√Dy[sub]n-1[/sub])=(px[sub]n-1[/sub]+qDy[sub]n-1[/sub])+√D(qx[sub]n-1[/sub]+py[sub]n-1[/sub])を得ます。[br][color=#9900ff][b][size=200]最小解(p,q)がわかれば、[br]x[sub]n[/sub]=px[sub]n-1[/sub]+qDy[sub]n-1[/sub][br]y[sub]n[/sub]=qx[sub]n-1[/sub]+py[sub]n-1[/sub]となる。[br][/size][/b][/color][br][br][color=#9900ff][u][b][size=150]質問：√3、√5、√7、√11の最小解は何でしょう？[br]また、近似値を効率的に出すプログラムはどうやればよいでしょうか。[br][/size][/b][/u][/color][br](x0,y0)=(1,0)は共通だね。[br]x[sup]2[/sup]-3y[sup]2[/sup]=1 最小解(x1,y1)=(2,1)から、行列は[2 3; 1 2][br]x[sup]2[/sup]-5y[sup]2[/sup]=1 最小解(x1,y1)=(9,4)から、行列は[9 20;4 9][br]x[sup]2[/sup]-7y[sup]2[/sup]=1 最小解(x1,y1)=(8,3)から、行列は[8 21; 3 8][br]x[sup]2[/sup]-11y[sup]2[/sup]=1 最小解(x1,y1)=(10,3)から、行列は[10 33; 3 10][br]√2,√3,√5,√7,√11の近似値を行列を14回実行して求めてみる。[br]Dと最小(p,q)の組、(D,p,q)をデータをして渡すと、14回かけ算したあとの解を返して表示する関数が[br]apx(x)だ。ついでに、標準で入っている平方根の近似値を並べて表示してみよう。[br][br][color=#0000ff]#[IN]julia[br]#==================================================[br][/color][color=#0000ff]using LinearAlgebra[br]data=[(2,3,2),(3,2,1),(5,9,4),(7,8,3),(11,10,3)] #(D,p,q) [br]function apx(x)[br] (D,p,q)=x[br] A = [ p q*D; q p ] #漸化式の行列[br] xy = [ 1; 0 ] #初期値[br] for i in 1:14[br] xy = A*xy[br] end[br] println("近似値 =",BigFloat(xy[1]/xy[2]))[br] println("sqrt $D =",sqrt(D))[br]end[br]for i in data[br] apx(i)[br]end[br]#==================================================[br][/color][OUT][br]近似値 =1.4142135623730951454746218587388284504413604736328125[br]sqrt 2 =1.4142135623730951[br]近似値 =1.73205080756887763726581397349946200847625732421875[br]sqrt 3 =1.7320508075688772[br]近似値 =2.236067977499789360962267892318777740001678466796875[br]sqrt 5 =2.23606797749979[br]近似値 =2.64575131106459071617109657381661236286163330078125[br]sqrt 7 =2.6457513110645907[br]近似値 =3.316624790355400254071582821779884397983551025390625[br]sqrt 11=3.3166247903554[br][br]14回行列計算を回しただけで、標準で入っている近似値に遜色ない値になるね。[br]すばらしい！

[b][size=150]＜高速化＞[/size][/b][br][br]n番目のｘｙを求める定義式(x[sub]1[/sub]+√Dy[sub]1[/sub])[sup]n[/sup]=x[sub]n[/sub]+√Dy[sub]n[/sub][br]を変形することで、x[sub]2n[/sub]=x[sub]n[/sub][sup]2[/sup]+Dy[sub]n[/sub][sup]2[/sup], y[sub]2n[/sub]= 2x[sub]n[/sub]y[sub]n[br][/sub]が得られる。（途中は相当省略）[br][br]すると、第１式x[sub]2n[/sub]=x[sub]n[/sub][sup]2[/sup]+Dy[sub]n[/sub][sup]2 [/sup]とx[sub]n[/sub][sup]2 [/sup]- Dy[sub]n[/sub][sup]2[/sup]=1から、[br][br]x[sub]2n[/sub]=2x[sub]n[/sub][sup]2[/sup]-1が得られるね。[br]だから、ペルの方程式にはDの値に関係なく、[br][br][color=#0000ff][b][size=150]x[sub]2n[/sub]=2x[sub]n[/sub][sup]2[/sup]-1,[br]y[sub]2n[/sub]= 2x[sub]n[/sub]y[sub]n[/sub][br][/size][/b][/color][br]というとても役立つ性質がある。[color=#0000ff][b]これで倍々で先に進めるね。[br][/b][/color][br][color=#9900ff][u][b][size=150]質問：この高速計算を使って、√２の近似値を求めるためのコードをかいてみよう。[br][/size][/b][/u][/color][br]juliaでやる。[br][color=#0000ff]#[IN]julia[br]#==================================================[br]using LinearAlgebra[br]xy = [ 3; 2 ][br]function fp2(xy)[br] (x,y)=xy[br] y=2x*y[br] x=2x^2-1[br] return (x,y)[br]end[br]function pell2(xy)[br] (x,y)=xy[br] return x^2-2*y^2[br]end[br]res=[][br]values=[ ][br]for i in 1:4[br] xy=fp2(xy)[br] push!(res,xy)[br] push!(values,pell2(xy))[br]end[br]println(res)[br]println(values)[br]#==================================================[br][/color][OUT][br]Any[(17, 12), (577, 408), (665857, 470832), (886731088897, 627013566048)][br]Any[1, 1, 1, 1]

おまけになりますが、D=2のペルの方程式の使い道がもう一つあります。[br][br]それは、「[b]平方三角数[/b]」です。[br][br]平方数はn[sup]２[/sup]が一般項、[br]三角数はm(m+1)/2が一般項でした。[br][br]これが両立する数は[color=#0000ff][b][size=150]n[sup]2[/sup]=m(m+1)/2[/size][/b][/color]ですね。それを平方三角数と呼びましょう。[br][br]D=2のペルの方程式の解の数列から平方三角数が作れるのです。[br]なぜか、x[sup]2[/sup]-2y[sup]2[/sup]=1の解ｘ，ｙに対して[b][color=#0000ff]m=(x-1)/2, n=y/2[/color][/b]を求めます。[br]すると、ｘ＝2m+1, y =2nなので、ペルの方程式に代入できます。[br](2m+1)[sup]2[/sup]-1=8n[sup]2[/sup][br]4m[sup]2[/sup]+4m=8n[sup]2[/sup][br]m(m+1)/2=n[sup]2[/sup][br][br][color=#9900ff][u][b][size=150]質問：ペルの方程式の解から平方三角数の数列を表示するにはどうしますか。[br][/size][/b][/u][/color]#[IN]Julia[br]#======================================================[br]#平方三角数[br]using LinearAlgebra[br]A = [ 3 4; 2 3 ][br]xy = [ 3; 2 ][br][color=#0000ff][b]function fp2(xy)[br] (x,y)=xy[br] m = (x-1)÷2; n = y÷2;[br] triSqr=BigInt(n)^2[br] println("(x,y)=",xy,",(m,n)=(",m,",",n,"),平方三角数=",triSqr)[br] return A*xy[br]end[br]for i in 1:20[br] xy=fp2(xy)[br]end[br][/b][/color]#=====================================================[br](x,y)=[3, 2],(m,n)=(1,1),平方三角数=1[br](x,y)=[17, 12],(m,n)=(8,6),平方三角数=36[br](x,y)=[99, 70],(m,n)=(49,35),平方三角数=1225[br](x,y)=[577, 408],(m,n)=(288,204),平方三角数=41616[br](x,y)=[3363, 2378],(m,n)=(1681,1189),平方三角数=1413721[br](x,y)=[19601, 13860],(m,n)=(9800,6930),平方三角数=48024900[br](x,y)=[114243, 80782],(m,n)=(57121,40391),平方三角数=1631432881[br](x,y)=[665857, 470832],(m,n)=(332928,235416),平方三角数=55420693056[br](x,y)=[3880899, 2744210],(m,n)=(1940449,1372105),平方三角数=1882672131025[br](x,y)=[22619537, 15994428],(m,n)=(11309768,7997214),平方三角数=63955431761796[br](x,y)=[131836323, 93222358],(m,n)=(65918161,46611179),平方三角数=2172602007770041[br](x,y)=[768398401, 543339720],(m,n)=(384199200,271669860),平方三角数=73804512832419600[br](x,y)=[4478554083, 3166815962],(m,n)=(2239277041,1583407981),平方三角数=2507180834294496361[br](x,y)=[26102926097, 18457556052],(m,n)=(13051463048,9228778026),平方三角数=85170343853180456676[br](x,y)=[152139002499, 107578520350],(m,n)=(76069501249,53789260175),平方三角数=2893284510173841030625[br](x,y)=[886731088897, 627013566048],(m,n)=(443365544448,313506783024),平方三角数=98286503002057414584576[br](x,y)=[5168247530883, 3654502875938],(m,n)=(2584123765441,1827251437969),平方三角数=3338847817559778254844961[br](x,y)=[30122754096401, 21300003689580],(m,n)=(15061377048200,10650001844790),平方三角数=113422539294030403250144100[br](x,y)=[175568277047523, 124145519261542],(m,n)=(87784138523761,62072759630771),平方三角数=3853027488179473932250054441[br](x,y)=[1023286908188737, 723573111879672],(m,n)=(511643454094368,361786555939836),平方三角数=130889512058808083293251706896[br]