def compute_convergents(sum_term, quotient_term):
"""
compute next term in convergents 1/(s + q)
returns a tuple (num, denum)
"""
num, denum = quotient_term
num += sum_term * denum
num, denum = denum, num
return (num, denum)
compute_convergents(2, (1, 2))
(2, 5)
def compute_square_root(order):
"""
1 should return 3/2,
2 should return 7/5, ...
"""
num, denum = 1, 2
for i in range(1, order):
num, denum = compute_convergents(2, (num, denum))
return (denum + num, denum)
from IPython.html.widgets import interact
interact(lambda order: print(compute_square_root(order)),
order=(1, 10))
(41, 29)
We can plot the convergence of the rational approximations compared to the "true value":
%matplotlib inline
from pylab import plot, sqrt, hlines
n = range(1, 101)
div = lambda t: t[0] / float(t[1])
div_sqrt = lambda order: div(compute_square_root(order))
plot(n, list(map(div_sqrt, n)))
hlines(sqrt(2), 1, 100, linestyles='dotted')
<matplotlib.collections.LineCollection at 0x106dc5128>
Indeed, this thing converges nicely to the square root of 2.
Let's try to get to our next step: evaluating continuous fractions with sequences other than just always the same number.
def compute_continued_fraction(order, starting_number, repeating_sequence):
"""
evaluate continued fraction defined √2 = [1;(2)] or √23 = [4;(1,3,1,8)]
"""
num, denum = 1, repeating_sequence[(order - 1) % (len(repeating_sequence))]
for i in range(order - 1, 0, -1):
lhs_factor = repeating_sequence[(i-1) % len(repeating_sequence)]
num, denum = compute_convergents(lhs_factor, (num, denum))
return (starting_number * denum + num, denum)
pdb.runcall(compute_continued_fraction, 3, 2, (1, 4, 3))
> <ipython-input-122-968fcd2ef05f>(5)compute_continued_fraction() -> num, denum = 1, repeating_sequence[(order - 1) % (len(repeating_sequence))] (Pdb) n > <ipython-input-122-968fcd2ef05f>(6)compute_continued_fraction() -> for i in range(order - 1, 0, -1): (Pdb) num, denum (1, 3) (Pdb) n > <ipython-input-122-968fcd2ef05f>(7)compute_continued_fraction() -> lhs_factor = repeating_sequence[(i-1) % len(repeating_sequence)] (Pdb) n > <ipython-input-122-968fcd2ef05f>(8)compute_continued_fraction() -> num, denum = compute_convergents(lhs_factor, (num, denum)) (Pdb) lhs_factor 4 (Pdb) n > <ipython-input-122-968fcd2ef05f>(6)compute_continued_fraction() -> for i in range(order - 1, 0, -1): (Pdb) num, denum (3, 13) (Pdb) n > <ipython-input-122-968fcd2ef05f>(7)compute_continued_fraction() -> lhs_factor = repeating_sequence[(i-1) % len(repeating_sequence)] (Pdb) n > <ipython-input-122-968fcd2ef05f>(8)compute_continued_fraction() -> num, denum = compute_convergents(lhs_factor, (num, denum)) (Pdb) lhs_factor 1 (Pdb) n > <ipython-input-122-968fcd2ef05f>(6)compute_continued_fraction() -> for i in range(order - 1, 0, -1): (Pdb) n > <ipython-input-122-968fcd2ef05f>(9)compute_continued_fraction() -> return (starting_number * denum + num, denum) (Pdb) n --Return-- > <ipython-input-122-968fcd2ef05f>(9)compute_continued_fraction()->(45, 16) -> return (starting_number * denum + num, denum) (Pdb) c
(45, 16)
compute_continued_fraction(3, 2, (1, 4, 3))
(45, 16)
Let's see if our code works for the square root of 2.
interact(lambda order: print(compute_continued_fraction(order, 1, (2,))),
order=(1, 10))
(577, 408)
Yes, seems to work. Let's see how this goes with the square root of 23.
interact(lambda order: print(compute_continued_fraction(order, 4, (1,3,1,8))),
order=(1, 10))
(211, 44)
import pdb
pdb.runcall(compute_continued_fraction, 1, 4, (1,3,1,8))
> <ipython-input-55-25f396e5f3b0>(5)compute_continued_fraction() -> repeating_sequence = repeating_sequence[::-1] (Pdb) q
sqrt(23)
4.7958315233127191
19/4
4.75
24/ 5
4.8
211 /44
4.795454545454546
Now on to the exponential:
e = [2; 1,2,1, 1,4,1, 1,6,1 , ... , 1,2k,1, ...]
The convergents are:
2, 3, 8/3, 11/4, 19/7, 87/32, 106/39, 193/71, 1264/465, 1457/536
interact(lambda order: print(compute_continued_fraction(order, 2, (1,2,1,1,4,1))),
order=(1, 10))
(492, 181)
Interlude: we write a function that generates the right coefficients until rank $n$.
def exponential_coefs(n):
coefs = []
for i in range(n):
if (i % 3) == 0 or (i % 3) == 2:
coefs.append(1)
else:
coefs.append(int((i + 2) / 3 * 2))
return coefs
exponential_coefs(9)
[1, 2, 1, 1, 4, 1, 1, 6, 1]
Now, let's compute the 10th partial of the exponential.
compute_continued_fraction(9, 2, exponential_coefs(9))
(1457, 536)
Now the 100th partial:
compute_continued_fraction(99, 2, exponential_coefs(99))
(6963524437876961749120273824619538346438023188214475670667, 2561737478789858711161539537921323010415623148113041714756)
And the solution is:
sum([int(s) for s in "6963524437876961749120273824619538346438023188214475670667"])
272