Finding continued fractions of pi with Stern-Brocot tree

user_16155472 0
Viewed 110

Im trying to use python to determine the continued fractions of pi by following the stern brocot tree. Its simple, if my estimation of pi is too high, take a left, if my estimation of pi is too low, take a right.

Im using mpmath to get arbitrary precision floating numbers, as python doesn't support that, but no matter what i set the decimal precision to using 'mp.dps', the continued fraction generation seems to stop once it hits 245850922/78256779.

In theory, it should only exit execution when it is equal to the current estimation for pi. So I tried increasing the decimal precision of mp.dps, however it still halts execution there.

have i reached a maximum amount of precision with mp.dps or is my approach inefficient? how can i make the continued fraction generation not cease at 245850922/78256779???

import mpmath as mp

mp.dps = 1000

def eval_stern_seq(seq):
  a,b,c,d,m,n=0,1,1,0,1,1
  for i in seq:
    if i=='L':
      c,d=m,n
    else:
      a,b=m,n
    m,n=a+c,b+d
  return m,n

seq = ''
while True:
  stern_frac = eval_stern_seq(seq)

  print(f"\n\ncurrent fraction:    {stern_frac[0]}/{stern_frac[1]}")

  print("current value:       " + mp.nstr(mp.fdiv(stern_frac[0],stern_frac[1]),n=mp.dps))
  print("pi (reference):      " + mp.nstr(mp.pi,n=mp.dps))

  if mp.fdiv(stern_frac[0],stern_frac[1]) > mp.pi:
    seq+='L'
  elif mp.fdiv(stern_frac[0],stern_frac[1]) < mp.pi:
    seq+='R'
  else:
    break

  input("\n\t(press enter to continue generation...)")
mpmath math python
2 Answers

You need to include the term prec=xxxxx in each fdiv operation, eg:mp.fdiv(stern_frac[0],stern_frac[1]) becomes mp.fdiv(stern_frac[0],stern_frac[1], prec=200)

user_146354300

You need from mpmath import mp, not import mpmath as mp. (Or you could import mpmath and then use mpmath.mp.dps = ...). See the documentation for additional information.

from mpmath import mp

mp.dps = 1000

def eval_stern_seq(seq):
    a, b, c, d, m, n = 0, 1, 1, 0, 1, 1
    for i in seq:
        if i == 'L':
            c, d = m, n
        else:
            a, b = m, n
        m, n = a + c, b + d
    return m, n

seq = ''
for _ in range(1000):
    stern_frac = eval_stern_seq(seq)
    print(f"{stern_frac[0]}/{stern_frac[1]}")
    if mp.fdiv(stern_frac[0], stern_frac[1]) > mp.pi:
        seq += 'L'
    elif mp.fdiv(stern_frac[0], stern_frac[1]) < mp.pi:
        seq += 'R'
    else:
        break

Output:

...
2810940413396914412349044750336209104844394887546/894750123057789325577450189803726071898054047077
2936686853946948334940603569095291581361715769163/934776458237287423525876798002671462625917371175
3062433294496982257532162387854374057879036650780/974802793416785521474303406201616853353780695273
user_120464090
Related