分数（フラクション）を入れ子にした[br][br][color=#0000ff][b][size=150]連分数(continued fraction)[/size][/b][/color]というのは知ってますね。[br][br]最近の言語では分数パッケージがあるので、その使い方になれるという目的もかねて、[br]連分数で平方根を近似してみよう。[br]それを連分数展開(continued fraction expansion)と言うこともある。[br][br]たとえば、[b][color=#0000ff][size=200]1+1/(2+1/(2+1/(2+1/......)))[/size][/color][/b][br]は、割り算を分数にすると、無限につながる連分数の式になる。 [br] [br][math]1+\frac{1}{2+\frac{1}{2+\frac{1}{2+\frac{1}{.....}}}}[/math][br][br]この値は２の平方根の近似値になることが知られているね。[br][br]分子はいつも1なので、「1/(a」のaの部分と先頭の整数部分だけ決めればよいので、[br]２の平方根の連分数のことを[b][color=#0000ff][size=150][1,2,2,2,2,......][/size][/color][/b]と書いたり、[br]2がサイクル１で続くだけなので、[b][size=150][color=#0000ff][1；[u]2[/u]][/color][/size][/b]などと書いたりできるね。[br][br][color=#9900ff][u][b][size=150]質問：√２の近似値を連分数を使ってpython,juliaで求めるにはどうしますか。[br][/size][/b][/u][/color][br]まずpythonで分数計算と小数表示をしてよう。[br]#[IN]Python[br]#===========================================================[br]#1+1/(2+1/(2+1/(2+1/......)))[br]#Fraction(numerator,denominator)は[br]#numerator,denominatorで分数オブジェクトを作るコンストラクター。[br]#整数の１、2を分数のone,twoとして作っておいて分数計算とする。[br][color=#0000ff][b]from fractions import Fraction[br]one=Fraction(1,1)[br]two=Fraction(2,1)[br]apx=one+one/(two+one/(two+one/two))[br]print(apx,":",float(apx))[br]apx2=one+one/(two+one/(two+one/(two+one/(two+one/(two+one/(two+one/(two+one/(two+one/two))))))))[br]print(apx2,":",float(apx2))[br]#===========================================================[br][/b][/color][OUT][br]17/12 : 1.4166666666666667[br]3363/2378 : 1.4142136248948696[br][br]Juliaではどうでしょうか。パッケージを宣言せずに、分数の線として//を使うだけで、[br]お手軽に分数計算ができたね。[br][color=#0000ff][b]#[IN]Julia[br]#===========================================================[br]one=1//1[br]two=2//1[br]apx=one+one/(two+one/(two+one/two))[br]println(apx,":",float(apx))[br]apx2=one+one/(two+one/(two+one/(two+one/(two+one/(two+one/(two+one/(two+one/(two+one/two))))))))[br]println(apx2,":",float(apx2))[br]#============================================================[br][/b][/color][OUT][br]17//12:1.4166666666666667[br]3363//2378:1.4142136248948696[br]

上記のやり方では、繰り返しを[br]ベタで、式にしている。[br]もっと、自由に使えるものはできないだろうか。[br][br][size=150][color=#9900ff][u]質問：√２の連分数展開を関数にして、好きな回数で近似するにはどうしたよいか？[br][/u][br][/color][color=#1e84cc][b]関数ｆで整数ｎから分数ｎ/1を作るようにすると、[br]f(1),f(2)とかくだけで、整数を分数に直せるね。[br][br][/b][/color][color=#0000ff]#[IN]Python[br][/color][/size][color=#0000ff]#=========================[br]from fractions import Fraction[br]f=lambda n:Fraction(n,1)[br]def continued_fraction(n):[br] part = f(2)[br] for i in range(1,n+1):[br] part= f(2)+f(1)/part[br] return f(1)+f(1)/part[br]for i in range(1,20):[br] p=continued_fraction(i)[br] print(float(p),":",p)[br]#=========================[br][/color][OUT][br]1.4 : 7/5[br]1.4166666666666667 : 17/12[br]1.4137931034482758 : 41/29[br]1.4142857142857144 : 99/70[br]1.4142011834319526 : 239/169[br]1.4142156862745099 : 577/408[br]1.4142131979695431 : 1393/985[br]1.4142136248948696 : 3363/2378[br]1.4142135516460548 : 8119/5741[br]1.4142135642135643 : 19601/13860[br]1.4142135620573204 : 47321/33461[br]1.4142135624272734 : 114243/80782[br]1.4142135623637995 : 275807/195025[br]1.4142135623746899 : 665857/470832[br]1.4142135623728214 : 1607521/1136689[br]1.414213562373142 : 3880899/2744210[br]1.414213562373087 : 9369319/6625109[br]1.4142135623730965 : 22619537/15994428[br]1.4142135623730947 : 54608393/38613965[br][br][color=#0000ff][b]#[IN]Julia[br]#=========================[br]f = n -> n//1[br]function continued_fraction(n)[br] part = f(2)[br] for i in 1:n[br] part= f(2)+f(1)/part[br] end[br] return f(1)+f(1)/part[br]end[br]for i in 1:20[br] p=continued_fraction(i)[br] println(float(p),":",p)[br]end[br]#=========================[br][/b][/color][OUT][br]1.4:7//5[br]1.4166666666666667:17//12[br]1.4137931034482758:41//29[br]1.4142857142857144:99//70[br]1.4142011834319526:239//169[br]1.4142156862745099:577//408[br]1.4142131979695431:1393//985[br]1.4142136248948696:3363//2378[br]1.4142135516460548:8119//5741[br]1.4142135642135643:19601//13860[br]1.4142135620573204:47321//33461[br]1.4142135624272734:114243//80782[br]1.4142135623637995:275807//195025[br]1.4142135623746899:665857//470832[br]1.4142135623728214:1607521//1136689[br]1.414213562373142:3880899//2744210[br]1.414213562373087:9369319//6625109[br]1.4142135623730965:22619537//15994428[br]1.4142135623730947:54608393//38613965[br]1.4142135623730951:131836323//93222358

√２だけの近似だけで使い道が限られるね。[br]そこで、もっとさまざまな平方根を求められるようにしてみよう。[br][br]√2の連分数展開は[1；[u]2[/u]]で、サイクル１だったが、他の平方根ではサイクルがいろいろる。[br][u]どうしてそうなるかはさておき[/u]、連分数展開を１つの分数や小数に直す関数を作ってみよう。[br]調べたい分数の長さをcntとしてサイクルを回して、分数の全長をcntにしたdistを作ろう。[br]有限な連分数の計算順位は最後尾が最優先になるので、逆順のrevを作ろう。[br]あとは、逆数にしてリストの数をたす。[br]これをくりかえそう。これで1未満のパートのpartができる。[br]最後のDataの先頭にある整数部分を最後にたそう。[br][br][color=#0000ff]#[IN]Python[br]#=================================================================[br]from fractions import Fraction[br]from itertools import cycle[br]f = lambda n:Fraction(n,1)[br]#(n,[cycle])[br]datum = [(1,[2]),(1,[1,2]),(2,[4]),(2,[1,1,1,4]),(3,[3,6]),(3,[1,1,1,1,6]),(4,[8]),(4,[2,1,3,1,2,8])][br]D =[2,3,5,7,11,13,17,19][br]def continued_fraction(Data,cnt):[br]  n,lst = Data[0],Data[1][br] dist = [][br] for i in cycle(lst):[br] dist.append(i)[br] if len(dist)>cnt-1 :[br] break[br] rev=list(reversed(dist))[br] print(dist,"-->",rev)[br] part=f(1)/f(rev[0])[br] for i in range(1,cnt):[br] part=f(rev[i])+part[br] part=f(1)/part[br] return f(n)+part[br]cnt=20[br]for i in range(len(D)):[br] p=continued_fraction(datum[i],cnt)[br] print(float(p),":",p,":√",D[i])[br][br]#=================================================================[br][/color][OUT][br][2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2] --> [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2][br]1.4142135642135643 : 19601/13860 :√ 2[br][1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1] --> [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1][br]1.7320512820512821 : 1351/780 :√ 3[br][4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4] --> [4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4][br]2.236067977499794 : 16692641/7465176 :√ 5[br][1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1] --> [1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1][br]2.645751633986928 : 2024/765 :√ 7[br][3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3] --> [3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3][br]3.3166247903554016 : 31521799/9504180 :√ 11[br][1, 1, 1, 1, 6, 1, 1, 1, 1, 6, 1] --> [1, 6, 1, 1, 1, 1, 6, 1, 1, 1, 1][br]3.6055514974433893 : 4936/1369 :√ 13[br][8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8] --> [8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8][br]4.123105625617661 : 41270070401/10009462320 :√ 17[br][2, 1, 3, 1, 2, 8, 2, 1, 3, 1, 2] --> [2, 1, 3, 1, 2, 8, 2, 1, 3, 1, 2][br]4.3588989441930615 : 57799/13260 :√ 19[br][br][color=#0000ff]juliaでもほとんど変わらない。[br]#[IN]Julia[br]#=================================================================[br]f = n -> n//1[br]#(n,[cycle])[br]datum = [(1,[2]),(1,[1,2]),(2,[4]),(2,[1,1,1,4]),(3,[3,6]),(3,[1,1,1,1,6]),(4,[8]),(4,[2,1,3,1,2,8])][br]D =[2,3,5,7,11,13,17,19][br]function continued_fraction(Data,cnt)[br]  (n,lst) = Data[br] dist=[ ][br] for (i, v) in enumerate(Iterators.cycle(lst))[br] push!(dist,v)[br] i > cnt-1 && break[br] end[br] rev=reverse(dist)[br] println(dist,"-->",rev)[br] part=f(1)/f(rev[1])[br] for i in 2:cnt[br] part=f(rev[i])+part[br] part=f(1)/part[br] end[br] return f(n)+part[br]end[br]cnt=20[br]for i in 1:length(D)[br] p=continued_fraction(datum[i],cnt)[br] println(float(p),":",p,":√",D[i])[br]end[br][/color]#=================================================================[color=#0000ff][br][OUT][/color][br]Any[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]-->Any[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2][br]1.4142135623730947:54608393//38613965:√2[br]Any[1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2]-->Any[2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1][br]1.732050807565499:716035//413403:√3[br]Any[4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]-->Any[4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4][br]2.23606797749979:7331474697802//3278735159921:√5[br]Any[1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4]-->Any[4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1][br]2.645751311063895:2388325//902702:√7[br]Any[3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3, 6]-->Any[6, 3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3][br]3.3166247903554:31539640338297//9509559365899:√11[br]Any[1, 1, 1, 1, 6, 1, 1, 1, 1, 6, 1, 1, 1, 1, 6, 1, 1, 1, 1, 6]-->Any[6, 1, 1, 1, 1, 6, 1, 1, 1, 1, 6, 1, 1, 1, 1, 6, 1, 1, 1, 1][br]3.6055512754637564:5564523//1543321:√13[br]Any[8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8]-->Any[8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8][br]4.123105625617661:6355271576489378132//1541379764079492545:√17[br]Any[2, 1, 3, 1, 2, 8, 2, 1, 3, 1, 2, 8, 2, 1, 3, 1, 2, 8, 2, 1]-->Any[1, 2, 8, 2, 1, 3, 1, 2, 8, 2, 1, 3, 1, 2, 8, 2, 1, 3, 1, 2][br]4.358898943540674:512445947//117563163:√19

これまでは、連分数情報をリストで与えると、それから平方根の近似値を分数、小数で表示する関数を考えてきた。[br]しかし、そもそも、その情報を作り出すことができないと、情報を入手して手入力したり、文字情報を数値情報に変換するなどの手間がかかってしまう。[br][br]Dだけから、√Dの連分数情報を作り出すにはどうしたらよいだろうか。[br][br]作り方に着目してみよう。[br]たとえば√2=1+1/(2+1/(2+1/(2+........)))だったから、[1；2]を得た。[br]有限分数99/70では、99/70=1+29/70=1+1/70/29=1+1/(2+12/29)[br] =1+1/(2+1/29/12)=1+1/(2+1/(2+5/12))=1+1/(2+1/(2+1/(2+2/5))=1+1/(2+1/(2+1/(2+1/(2+1/2))))[br]これを無理数√2でやると、[br]√2/1=1+(√2-1)/1=1+1/(√2+1)=1+1/(2+(√2-1))=1+1/(2+1/(√2+1))=.....[br]もとの数A/Bを整数部分＋逆数にして、分数部分をまた整数部分＋逆数にしてという操作は[br]AとBについてユークリッドの互除法をやっていると同じだとわかるね。[br]手計算で逆数を有理化していけば、単位分数を作っていけば、連分数情報が得られることがわかるね。[br]しかし、いちいち、逆数にして有理化するのはコード作成が大変そうだ。[br]そこで、サイクルの規則性についてしらべよう。[br][br][b][size=150]＜D=n[sup]2[/sup]+1のとき＞[br][/size][/b]Dが平方数n[sup]2[/sup]に1をたした形で、[b][size=150][ n；(2n)][br][/size][/b]例）D=2 ならn=1から [1；(2)][br] D=5 ならn=2から [2；(4)][br] D=10 ならn=3から [3；(6)][br] D=17 ならn=4から [4；(8)][br] D=26 ならn=5から [5；(10)][br] D=37 ならn=6から [6；(12)][br]理由）[br]　　x=√(n[sup]2[/sup]+1)=n+(√(n[sup]2[/sup]+1) -n) =n + 1/x[sub]1[/sub]とする。[br] x[sub]1[/sub]=1/(√(n[sup]2[/sup]+1) -n)=(√(n[sup]2[/sup]+1)+n)=2n+(√(n[sup]2[/sup]+k)-n)[br] = d + 1/x[sub]2[/sub]とおける。だから、サイクルの先頭は2nだね。[br]　　x[sub]2[/sub]=1/(√(n[sup]2[/sup]+1)-n)=x[sub]1[/sub]となるから、サイクル長は１になるね。[br][br][b][size=150]＜D=n[sup]2[/sup]+kでd=2n/kが整数のとき＞[br][/size][/b]Dが平方数n[sup]2[/sup]にkをたした形で、d=2n/kが整数なら、[b][size=150][n；(d, 2n)][br][/size][/b]例）D=3 ならn=1, k=1から2n/k=1　サイクル=(1,2)[br] D=11ならn=3, k=2から2n/k=3 　サイクル=(3,6)[br]これは、Dが素数じゃなくても使える。[br] D=6ならn=2,k=2から2n/k=2 サイクル(2,4)[br] D= 8ならn=2,k=4から2n/k=1 サイクル(1,4)[br] D=12ならn=3,k=3から2n/k=2 サイクル(2,6)[br] D=15ならn=3,k=6から2n/k=1 サイクル(1,6)[br] D=18ならn=4,k=2から2n/k=4 サイクル(4,8)[br] D=20ならn=4,k=4から2n/k=2 サイクル(2,8）[br] D=24ならn=4,k=8から2n/k=1 サイクル(1,8)[br]理由）[br]　　x=√(n[sup]2[/sup]+k)=n+(√(n[sup]2[/sup]+k) -n) =n + 1/x[sub]1[/sub]とする。[br] x[sub]1[/sub]=1/(√(n[sup]2[/sup]+k) -n)=(√(n[sup]2[/sup]+k)+n)/(√(n[sup]2[/sup]+k) +n)(√(n[sup]2[/sup]+k) -n)=(√(n[sup]2[/sup]+k)+n)/k=2n/k+(√(n[sup]2[/sup]+k)-n)/k[br] = d + 1/x[sub]2[/sub]とおける。だから、サイクルの先頭はdだね。[br]　　x[sub]2[/sub]=k/(√(n[sup]2[/sup]+k)-n)=k(√(n[sup]2[/sup]+k)+n)/(√(n[sup]2[/sup]+k)-n)(√(n[sup]2[/sup]+k)+n)=k(√(n[sup]2[/sup]+k)+n)/k=√(n[sup]2[/sup]+k)+n[br] =2n+√(n[sup]2[/sup]+k)-n=2n+1/x[sub]3[/sub]とおける。[br]　　x3=1/(√(n[sup]2[/sup]+k) -n)=x1となるから、2nがサイクルの末尾になる。これで周期は２となる。[br][br][b][size=150]＜D=n[sup]2[/sup]-kでd=2n/kが整数のとき＞[br][/size][/b]Dが平方数n[sup]2[/sup]にkをたした形で、d=2n/kが整数なら、[b][size=150][ n-1；(1,d-2,1, 2(n-1))][br][/size][/b]例）D=7 ならn=3, k=2からd=2n/k=3 。これからn-1=2, d-2=1より、 [2；(1,1,1,4)][br] D=14ならn=4, k=2からd=2n/k=4 。これからn-1=3, d-2=2より、 [3；(1,2,1,6)][br] D=23ならn=5, k=2からd=2n/k=5 。これからn-1=4, d-2=3より、 [4；(1,3,1,8)][br] D=32ならn=6, k=4からd=2n/k=3 。これからn-1=5, d-2=1より、 [5；(1,1,1,10)][br] D=33ならn=6, k=3からd=2n/k=4 。これからn-1=5, d-2=2より、 [5；(1,2,1,10)][br] D=34ならn=6, k=2からd=2n/k=6 。これからn-1=5, d-2=4より、 [5；(1,4,1,10)][br] D=47ならn=7, k=2からd=2n/k=7 。これからn-1=6, d-2=5より、 [6；(1,5,1,12)][br] D=32ならn=4, k=2からd=2n/k=4 。これからn-1=3, d-2=2より、 [3；(1,2,1,6)[br]理由）[br]　上記と同様の計算をすると、x5=x1となり、サイクル長が４になる。[br][b][size=150][br]＜D=n[sup]2[/sup]+1のとき＞[br][/size][/b]Dが平方数n[sup]2[/sup]から1を引いた形で、[b][size=150][ n-1；(1,2(n-1))][br][/size][/b]例）D=3 ならn=2から [1；(1,2)][br] D=8 ならn=3から [2；(1,4)][br] D=15 ならn=4から [3；(1,6)][br] D=24 ならn=5から [4；(1,8)][br] D=35 ならn=6から [5；(1,10)] D=48 ならn=7から [6；(1,12)][br]理由）[br]　上記と同様の計算で、k=1にすると、d=2n/k=2nより、d-2=2n-2=2(n-1)となるので。[br]　次に１にもどる。[br][br]このように、Dに近い平方数ｎに着目することで、規則性が使えるが、2n/kが整数にならない場合がある。[br]D=13,19,21,22,28,29,.....などは上記のタイプからはずれるので、規則性につながらない。[br]

これまでの連分数は逆数利用で、分子を１とする形に限定してきた。[br]しかし、そうすることで見えてくるおもしろさはあるが、平方根の近似値を求めるための規則性としては[br]あまり、単純化は簡単そうではない。[br][br][math]\sqrt{n}=m+\frac{n-m^2}{m+\sqrt{n}}[/math] [br][br]そこで、上記のような分子が１ではない形でのくり返しをしてみよう。[br]√n=√n(m+√n)/(m+√n)=(n+m√n)/(m+√n)=(m[sup]2[/sup]+m√n+n-m[sup]2[/sup])/(m+√n)=(m(m+√n)+(n-m[sup]2[/sup]))/(m+√n)[br] = m + (n-m[sup]2[/sup])/(m + √n)[br]と変形できる。[br][br]たとえば、[br]n=11なら、3^2<11<4^2から、m=3。t=n-m^2=2。2m=6から、[br]√11=3+(√11-3)=(3+2/(3+√11))となる。　√11=3+2/(3+(3+2/(3+√11)))==> 3+2/(6+2/(6+2/....))[br]すると、[br]part=2/6, part=2/(6+part)をくり返し、最後に3+partする。[br]一般には、[br]mはsqr(n)の整数部分。t=n-m^2とすると、[br]part=t/2m, part=t/(2m+part)をくり返し、最後にm+partをする。[br][br][color=#9900ff][b][u][size=150]質問：20以下の素数の平方根の近似値を求めるにはどんなコードをかきますか？[br][/size][/u][/b][/color][br]#[IN]python[br]#============================================[br]from fractions import Fraction[br]f = lambda n:Fraction(n,1)[br]Ds = [2,3,5,7,11,13,17,19] #([n])[br]def irregular_fraction(n ,cnt):[br] m=0[br] while n-m**2>0:[br] m +=1[br] break[br] t = n-m**2 [br] part=f(t)/f(2*m)[br] for i in range(1,cnt):[br] part=f(t)/(f(2*m)+part)[br] return f(m)+part[br]cnt2=20[br]cnt=11[br]for i in range(len(Ds)):[br] q=irregular_fraction(Ds[i],cnt2)[br] print(float(q),":",q,":非正則",Ds[i])[br] p=continued_fraction(datum[i],cnt)[br] print(float(p),":",p,":正則",D[i])[br]#============================================[br][OUT][br]1.4142135623730947 : 54608393/38613965 :非正則 2[br]1.4142135642135643 : 19601/13860 :正則 2[br]1.732050807565499 : 716035/413403 :非正則 3[br]1.7320512820512821 : 1351/780 :正則 3[br]2.236067970034716 : 12238/5473 :非正則 5[br]2.236067977499794 : 16692641/7465176 :正則 5[br]2.6457510161477704 : 306335887/115784095 :非正則 7[br]2.645751633986928 : 2024/765 :正則 7[br]3.3166108052804133 : 10634068259/3206305739 :非正則 11[br]3.3166247903554016 : 31521799/9504180 :正則 11[br]3.605505252249399 : 20241905/5614166 :非正則 13[br]3.6055514974433893 : 4936/1369 :正則 13[br]4.122853067245746 : 378927653/91909085 :非正則 17[br]4.123105625617661 : 41270070401/10009462320 :正則 17[br]4.3584203886767705 : 998235750499/229036132699 :非正則 19[br]4.3588989441930615 : 57799/13260 :正則 19[br][br]最初は正則連分数と同じくcnt=11回でまわしたのだけど、[br]非正則連分数の精度が悪いため、cnt2=20と回す数を倍ちかくにしたら、どうにか同じくらいの精度に[br]なった。[br][br]非正則の方がコードにしやすい分、精度が下がってしまったということだ。[br][color=#0000ff][b][size=150]分かりやすいさは、精度や速度に比例しない[/size][/b][/color]かも[br]という教訓？[br]