Continued fraction of pi by hand

edit

You found a bug ! Here is the explanation:

When you want the integer part of the (symbolic) number defined by:

sage: x = 1/(1/(1/(1/(1/(1/(1/(1/(1/(1/(1/(1/(1/(1/(pi - 3) - 7) - 15) - 1) - 292) - 1) - 1) - 1) - 2) - 1) - 3) - 1) - 14) - 2)

Sage calls the method sage: x.__int__(), of which you can get the source code by typing:

sage: x.__int__??

Then, since x is a symbolic number, you will see that Sage has to get a numerical extimation of x. For this, it uses RIF, the Real Interval Field with 53 bits of precision to have a certified answer.

But with your x, the endpoints of the interval are:

sage: RIF(x)
[-infinity .. +infinity]

So that we lost all precision, then we can not get the integer part, even if we simplify the expression first:

sage: y = x.full_simplify() ; y
-(25510582*pi - 80143857)/(52746197*pi - 165707065)
sage: RIF(y)
[-infinity .. +infinity]

So, this is not due to the accumulation of errors in the number of computations, but due to the fact that 53 bits is not enough.

Here is a workaround for your case:

To get a certified value of int(x), you can extend the precision of the Real Interval Field you use, and check that the upper value of the interval has the same integer part than the lower . For example:

sage: Field = RealIntervalField(100) ; Field
Real Interval Field with 100 bits of precision
sage: Field(x)
1.7069000446718?
sage: int(Field(x).upper()) == int(Field(x).lower())
True
sage: int(Field(x).upper())
1

I let you write a my_int function tha adapts the precision to the entry ;)

I changed your code slightly to use in Python with sympy and get the right continued fraction for pi up to term 1000. It can be easily extended to higher numbers. Calculation time is a little over 2 seconds on my machine. Here is the code:

import sympy as sp

p=sp.N(sp.pi, 5000)

n=2000
x=range(n+1)
a=range(n)
x[0]=p
for i in xrange(n):
    a[i] = int(x[i])
    x[i+1]=1/(x[i]-a[i])
    print(a[i]),