Fractionで平方根近似
非正則連分数で素数の平方根を近似しよう
１．√2を分数を入れ子にして表す
分数（フラクション）を入れ子にした
連分数(continued fraction)というのは知ってますね。
最近の言語では分数パッケージがあるので、その使い方になれるという目的もかねて、
連分数で平方根を近似してみよう。
それを連分数展開(continued fraction expansion)と言うこともある。
たとえば、1+1/(2+1/(2+1/(2+1/......)))
は、割り算を分数にすると、無限につながる連分数の式になる。
この値は２の平方根の近似値になることが知られているね。
分子はいつも1なので、「1/(a」のaの部分と先頭の整数部分だけ決めればよいので、
２の平方根の連分数のことを[1,2,2,2,2,......]と書いたり、
2がサイクル１で続くだけなので、[1；2]などと書いたりできるね。
質問：√２の近似値を連分数を使ってpython,juliaで求めるにはどうしますか。
まずpythonで分数計算と小数表示をしてよう。
#[IN]Python
#===========================================================
#1+1/(2+1/(2+1/(2+1/......)))
#Fraction(numerator,denominator)は
#numerator,denominatorで分数オブジェクトを作るコンストラクター。
#整数の１、2を分数のone,twoとして作っておいて分数計算とする。
from fractions import Fraction
one=Fraction(1,1)
two=Fraction(2,1)
apx=one+one/(two+one/(two+one/two))
print(apx,":",float(apx))
apx2=one+one/(two+one/(two+one/(two+one/(two+one/(two+one/(two+one/(two+one/(two+one/two))))))))
print(apx2,":",float(apx2))
#===========================================================
[OUT]
17/12 : 1.4166666666666667
3363/2378 : 1.4142136248948696
Juliaではどうでしょうか。パッケージを宣言せずに、分数の線として//を使うだけで、
お手軽に分数計算ができたね。
#[IN]Julia
#===========================================================
one=1//1
two=2//1
apx=one+one/(two+one/(two+one/two))
println(apx,":",float(apx))
apx2=one+one/(two+one/(two+one/(two+one/(two+one/(two+one/(two+one/(two+one/(two+one/two))))))))
println(apx2,":",float(apx2))
#============================================================
[OUT]
17//12:1.4166666666666667
3363//2378:1.4142136248948696
２．回数指定で√２を近似する関数
上記のやり方では、繰り返しを
ベタで、式にしている。
もっと、自由に使えるものはできないだろうか。
質問：√２の連分数展開を関数にして、好きな回数で近似するにはどうしたよいか？
関数ｆで整数ｎから分数ｎ/1を作るようにすると、
f(1),f(2)とかくだけで、整数を分数に直せるね。
#[IN]Python
#=========================
from fractions import Fraction
f=lambda n:Fraction(n,1)
def continued_fraction(n):
part = f(2)
for i in range(1,n+1):
part= f(2)+f(1)/part
return f(1)+f(1)/part
for i in range(1,20):
p=continued_fraction(i)
print(float(p),":",p)
#=========================
[OUT]
1.4 : 7/5
1.4166666666666667 : 17/12
1.4137931034482758 : 41/29
1.4142857142857144 : 99/70
1.4142011834319526 : 239/169
1.4142156862745099 : 577/408
1.4142131979695431 : 1393/985
1.4142136248948696 : 3363/2378
1.4142135516460548 : 8119/5741
1.4142135642135643 : 19601/13860
1.4142135620573204 : 47321/33461
1.4142135624272734 : 114243/80782
1.4142135623637995 : 275807/195025
1.4142135623746899 : 665857/470832
1.4142135623728214 : 1607521/1136689
1.414213562373142 : 3880899/2744210
1.414213562373087 : 9369319/6625109
1.4142135623730965 : 22619537/15994428
1.4142135623730947 : 54608393/38613965
#[IN]Julia
#=========================
f = n -> n//1
function continued_fraction(n)
part = f(2)
for i in 1:n
part= f(2)+f(1)/part
end
return f(1)+f(1)/part
end
for i in 1:20
p=continued_fraction(i)
println(float(p),":",p)
end
#=========================
[OUT]
1.4:7//5
1.4166666666666667:17//12
1.4137931034482758:41//29
1.4142857142857144:99//70
1.4142011834319526:239//169
1.4142156862745099:577//408
1.4142131979695431:1393//985
1.4142136248948696:3363//2378
1.4142135516460548:8119//5741
1.4142135642135643:19601//13860
1.4142135620573204:47321//33461
1.4142135624272734:114243//80782
1.4142135623637995:275807//195025
1.4142135623746899:665857//470832
1.4142135623728214:1607521//1136689
1.414213562373142:3880899//2744210
1.414213562373087:9369319//6625109
1.4142135623730965:22619537//15994428
1.4142135623730947:54608393//38613965
1.4142135623730951:131836323//93222358
３．一般の平方根を近似する関数
√２だけの近似だけで使い道が限られるね。
そこで、もっとさまざまな平方根を求められるようにしてみよう。
√2の連分数展開は[1；2]で、サイクル１だったが、他の平方根ではサイクルがいろいろる。
どうしてそうなるかはさておき、連分数展開を１つの分数や小数に直す関数を作ってみよう。
調べたい分数の長さをcntとしてサイクルを回して、分数の全長をcntにしたdistを作ろう。
有限な連分数の計算順位は最後尾が最優先になるので、逆順のrevを作ろう。
あとは、逆数にしてリストの数をたす。
これをくりかえそう。これで1未満のパートのpartができる。
最後のDataの先頭にある整数部分を最後にたそう。
#[IN]Python
#=================================================================
from fractions import Fraction
from itertools import cycle
f = lambda n:Fraction(n,1)
#(n,[cycle])
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])]
D =[2,3,5,7,11,13,17,19]
def continued_fraction(Data,cnt):
n,lst = Data[0],Data[1]
dist = []
for i in cycle(lst):
dist.append(i)
if len(dist)>cnt-1 :
break
rev=list(reversed(dist))
print(dist,"-->",rev)
part=f(1)/f(rev[0])
for i in range(1,cnt):
part=f(rev[i])+part
part=f(1)/part
return f(n)+part
cnt=20
for i in range(len(D)):
p=continued_fraction(datum[i],cnt)
print(float(p),":",p,":√",D[i])
#=================================================================
[OUT]
[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2] --> [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
1.4142135642135643 : 19601/13860 :√ 2
[1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1] --> [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1]
1.7320512820512821 : 1351/780 :√ 3
[4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4] --> [4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]
2.236067977499794 : 16692641/7465176 :√ 5
[1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1] --> [1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1]
2.645751633986928 : 2024/765 :√ 7
[3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3] --> [3, 6, 3, 6, 3, 6, 3, 6, 3, 6, 3]
3.3166247903554016 : 31521799/9504180 :√ 11
[1, 1, 1, 1, 6, 1, 1, 1, 1, 6, 1] --> [1, 6, 1, 1, 1, 1, 6, 1, 1, 1, 1]
3.6055514974433893 : 4936/1369 :√ 13
[8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8] --> [8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8]
4.123105625617661 : 41270070401/10009462320 :√ 17
[2, 1, 3, 1, 2, 8, 2, 1, 3, 1, 2] --> [2, 1, 3, 1, 2, 8, 2, 1, 3, 1, 2]
4.3588989441930615 : 57799/13260 :√ 19
juliaでもほとんど変わらない。
#[IN]Julia
#=================================================================
f = n -> n//1
#(n,[cycle])
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])]
D =[2,3,5,7,11,13,17,19]
function continued_fraction(Data,cnt)
(n,lst) = Data
dist=[ ]
for (i, v) in enumerate(Iterators.cycle(lst))
push!(dist,v)
i > cnt-1 && break
end
rev=reverse(dist)
println(dist,"-->",rev)
part=f(1)/f(rev[1])
for i in 2:cnt
part=f(rev[i])+part
part=f(1)/part
end
return f(n)+part
end
cnt=20
for i in 1:length(D)
p=continued_fraction(datum[i],cnt)
println(float(p),":",p,":√",D[i])
end
#=================================================================
[OUT]
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]
1.4142135623730947:54608393//38613965:√2
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]
1.732050807565499:716035//413403:√3
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]
2.23606797749979:7331474697802//3278735159921:√5
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]
2.645751311063895:2388325//902702:√7
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]
3.3166247903554:31539640338297//9509559365899:√11
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]
3.6055512754637564:5564523//1543321:√13
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]
4.123105625617661:6355271576489378132//1541379764079492545:√17
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]
4.358898943540674:512445947//117563163:√19
4。Dから√Dの連分数情報を作り出す。
これまでは、連分数情報をリストで与えると、それから平方根の近似値を分数、小数で表示する関数を考えてきた。
しかし、そもそも、その情報を作り出すことができないと、情報を入手して手入力したり、文字情報を数値情報に変換するなどの手間がかかってしまう。
Dだけから、√Dの連分数情報を作り出すにはどうしたらよいだろうか。
作り方に着目してみよう。
たとえば√2=1+1/(2+1/(2+1/(2+........)))だったから、[1；2]を得た。
有限分数99/70では、99/70=1+29/70=1+1/70/29=1+1/(2+12/29)
=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))))
これを無理数√2でやると、
√2/1=1+(√2-1)/1=1+1/(√2+1)=1+1/(2+(√2-1))=1+1/(2+1/(√2+1))=.....
もとの数A/Bを整数部分＋逆数にして、分数部分をまた整数部分＋逆数にしてという操作は
AとBについてユークリッドの互除法をやっていると同じだとわかるね。
手計算で逆数を有理化していけば、単位分数を作っていけば、連分数情報が得られることがわかるね。
しかし、いちいち、逆数にして有理化するのはコード作成が大変そうだ。
そこで、サイクルの規則性についてしらべよう。
＜D=n2+1のとき＞
Dが平方数n2に1をたした形で、[ n；(2n)]
例）D=2 ならn=1から [1；(2)]
D=5 ならn=2から [2；(4)]
D=10 ならn=3から [3；(6)]
D=17 ならn=4から [4；(8)]
D=26 ならn=5から [5；(10)]
D=37 ならn=6から [6；(12)]
理由）
x=√(n2+1)=n+(√(n2+1) -n) =n + 1/x1とする。
x1=1/(√(n2+1) -n)=(√(n2+1)+n)=2n+(√(n2+k)-n)
= d + 1/x2とおける。だから、サイクルの先頭は2nだね。
x2=1/(√(n2+1)-n)=x1となるから、サイクル長は１になるね。
＜D=n2+kでd=2n/kが整数のとき＞
Dが平方数n2にkをたした形で、d=2n/kが整数なら、[n；(d, 2n)]
例）D=3 ならn=1, k=1から2n/k=1 サイクル=(1,2)
D=11ならn=3, k=2から2n/k=3 サイクル=(3,6)
これは、Dが素数じゃなくても使える。
D=6ならn=2,k=2から2n/k=2 サイクル(2,4)
D= 8ならn=2,k=4から2n/k=1 サイクル(1,4)
D=12ならn=3,k=3から2n/k=2 サイクル(2,6)
D=15ならn=3,k=6から2n/k=1 サイクル(1,6)
D=18ならn=4,k=2から2n/k=4 サイクル(4,8)
D=20ならn=4,k=4から2n/k=2 サイクル(2,8）
D=24ならn=4,k=8から2n/k=1 サイクル(1,8)
理由）
x=√(n2+k)=n+(√(n2+k) -n) =n + 1/x1とする。
x1=1/(√(n2+k) -n)=(√(n2+k)+n)/(√(n2+k) +n)(√(n2+k) -n)=(√(n2+k)+n)/k=2n/k+(√(n2+k)-n)/k
= d + 1/x2とおける。だから、サイクルの先頭はdだね。
x2=k/(√(n2+k)-n)=k(√(n2+k)+n)/(√(n2+k)-n)(√(n2+k)+n)=k(√(n2+k)+n)/k=√(n2+k)+n
=2n+√(n2+k)-n=2n+1/x3とおける。
x3=1/(√(n2+k) -n)=x1となるから、2nがサイクルの末尾になる。これで周期は２となる。
＜D=n2-kでd=2n/kが整数のとき＞
Dが平方数n2にkをたした形で、d=2n/kが整数なら、[ n-1；(1,d-2,1, 2(n-1))]
例）D=7 ならn=3, k=2からd=2n/k=3 。これからn-1=2, d-2=1より、 [2；(1,1,1,4)]
D=14ならn=4, k=2からd=2n/k=4 。これからn-1=3, d-2=2より、 [3；(1,2,1,6)]
D=23ならn=5, k=2からd=2n/k=5 。これからn-1=4, d-2=3より、 [4；(1,3,1,8)]
D=32ならn=6, k=4からd=2n/k=3 。これからn-1=5, d-2=1より、 [5；(1,1,1,10)]
D=33ならn=6, k=3からd=2n/k=4 。これからn-1=5, d-2=2より、 [5；(1,2,1,10)]
D=34ならn=6, k=2からd=2n/k=6 。これからn-1=5, d-2=4より、 [5；(1,4,1,10)]
D=47ならn=7, k=2からd=2n/k=7 。これからn-1=6, d-2=5より、 [6；(1,5,1,12)]
D=32ならn=4, k=2からd=2n/k=4 。これからn-1=3, d-2=2より、 [3；(1,2,1,6)
理由）
上記と同様の計算をすると、x5=x1となり、サイクル長が４になる。
＜D=n2+1のとき＞
Dが平方数n2から1を引いた形で、[ n-1；(1,2(n-1))]
例）D=3 ならn=2から [1；(1,2)]
D=8 ならn=3から [2；(1,4)]
D=15 ならn=4から [3；(1,6)]
D=24 ならn=5から [4；(1,8)]
D=35 ならn=6から [5；(1,10)] D=48 ならn=7から [6；(1,12)]
理由）
上記と同様の計算で、k=1にすると、d=2n/k=2nより、d-2=2n-2=2(n-1)となるので。
次に１にもどる。
このように、Dに近い平方数ｎに着目することで、規則性が使えるが、2n/kが整数にならない場合がある。
D=13,19,21,22,28,29,.....などは上記のタイプからはずれるので、規則性につながらない。
５．非正則な連分数で平方根を近似する
これまでの連分数は逆数利用で、分子を１とする形に限定してきた。
しかし、そうすることで見えてくるおもしろさはあるが、平方根の近似値を求めるための規則性としては
あまり、単純化は簡単そうではない。
そこで、上記のような分子が１ではない形でのくり返しをしてみよう。
√n=√n(m+√n)/(m+√n)=(n+m√n)/(m+√n)=(m2+m√n+n-m2)/(m+√n)=(m(m+√n)+(n-m2))/(m+√n)
= m + (n-m2)/(m + √n)
と変形できる。
たとえば、
n=11なら、3^2<11<4^2から、m=3。t=n-m^2=2。2m=6から、
√11=3+(√11-3)=(3+2/(3+√11))となる。 √11=3+2/(3+(3+2/(3+√11)))==> 3+2/(6+2/(6+2/....))
すると、
part=2/6, part=2/(6+part)をくり返し、最後に3+partする。
一般には、
mはsqr(n)の整数部分。t=n-m^2とすると、
part=t/2m, part=t/(2m+part)をくり返し、最後にm+partをする。
質問：20以下の素数の平方根の近似値を求めるにはどんなコードをかきますか？
#[IN]python
#============================================
from fractions import Fraction
f = lambda n:Fraction(n,1)
Ds = [2,3,5,7,11,13,17,19] #([n])
def irregular_fraction(n ,cnt):
m=0
while n-m**2>0:
m +=1
break
t = n-m**2
part=f(t)/f(2*m)
for i in range(1,cnt):
part=f(t)/(f(2*m)+part)
return f(m)+part
cnt2=20
cnt=11
for i in range(len(Ds)):
q=irregular_fraction(Ds[i],cnt2)
print(float(q),":",q,":非正則",Ds[i])
p=continued_fraction(datum[i],cnt)
print(float(p),":",p,":正則",D[i])
#============================================
[OUT]
1.4142135623730947 : 54608393/38613965 :非正則 2
1.4142135642135643 : 19601/13860 :正則 2
1.732050807565499 : 716035/413403 :非正則 3
1.7320512820512821 : 1351/780 :正則 3
2.236067970034716 : 12238/5473 :非正則 5
2.236067977499794 : 16692641/7465176 :正則 5
2.6457510161477704 : 306335887/115784095 :非正則 7
2.645751633986928 : 2024/765 :正則 7
3.3166108052804133 : 10634068259/3206305739 :非正則 11
3.3166247903554016 : 31521799/9504180 :正則 11
3.605505252249399 : 20241905/5614166 :非正則 13
3.6055514974433893 : 4936/1369 :正則 13
4.122853067245746 : 378927653/91909085 :非正則 17
4.123105625617661 : 41270070401/10009462320 :正則 17
4.3584203886767705 : 998235750499/229036132699 :非正則 19
4.3588989441930615 : 57799/13260 :正則 19
最初は正則連分数と同じくcnt=11回でまわしたのだけど、
非正則連分数の精度が悪いため、cnt2=20と回す数を倍ちかくにしたら、どうにか同じくらいの精度に
なった。
非正則の方がコードにしやすい分、精度が下がってしまったということだ。
分かりやすいさは、精度や速度に比例しないかも
という教訓？