(有理算術演算: 高精度数値計算のためのアルゴリズムとプログラミング 2023) を参考にして、連分数を計算するコードを書いた。
- 関数:
get_continued_fraction で連分数展開 [a0, a1, …] を取得する
- 関数:
get_approx_fraction で連分数展開から近似分数を計算する
import math
from fractions import Fraction
def get_continued_fraction(n):
"""入力は素数と仮定して、循環の最初の周期を返却する"""
r = 0
s = 1
a = int(math.floor(math.sqrt(n)))
a0 = a
lst_a = [a]
cnt = 100
while 0 < cnt:
_s = (n-(r-a*s)**2) / s
_r = -(r-a*s)
a = int(math.floor( (math.sqrt(n) + _r) / _s))
r = _r
s = _s
cnt -= 1
lst_a.append(a)
if a == 2*a0:
break
return lst_a
def get_approx_fraction(cfrac):
assert 0 < len(cfrac)
if 1 == len(cfrac):
return Fraction(cfrac[0], 1)
elif 2 == len(cfrac):
return Fraction(cfrac[0] * cfrac[1] + 1, cfrac[1])
else:
pi_1 = cfrac[0] * cfrac[1] + 1
pi = cfrac[2] * pi_1 + cfrac[0]
qi_1 = cfrac[1]
qi = cfrac[2] * qi_1 + 1
for ai in cfrac[3:]:
pi, pi_1 = ai * pi + pi_1, pi
qi, qi_1 = ai * qi + qi_1, qi
return Fraction(pi, qi)