Finding continued fractions of pi with Stern-Brocot tree-CodePudding

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...)")

CodePudding user response：

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)

CodePudding user response：

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