# Continued fraction of pi by hand

Hello, I'm a newbie in Sage. I tried to compute (sagenb.org) the continued fraction of pi by hand with the following commands

n=15
x=range(n+1);a=range(n)
x[0]=pi
for i in range(n):
a[i] = int(x[i])
x[i+1]=1/(x[i]-a[i])
print(a[i])


the answer (14 terms) is correct, but then I get

Traceback (click to the left of this block for traceback)
...
ValueError: Calling floor() on infinity or NaN

Am I doing something wrong?

Cheers

edit retag close merge delete

Sort by » oldest newest most voted

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 ;)

more

1

Hopefully, when ticket #14567 will be reviewed continued_fraction(pi) will work like a charm ;-)

( 2014-11-22 23:40:24 -0500 )edit

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]),

more

Thanks for the quick answer. I extended the precision


R=RealIntervalField(100)
n=28
x=range(n+1);a=range(n)
x[0]=R(pi).upper()
for i in range(n):
a[i] = int(x[i])
x[i+1]=1/(x[i]-a[i])
print(a[i])


but now a[28] is wrong (compared to continued_fraction). So 100 bits is not enough. How I discover it?

more

When you write x[0]=R(pi).upper() , x[0] stop being the exact pi but a rational approximation. Hence at some point you will get a wrong answer, since the rational and pi do not have the same expression. What i was suggesting is : keep x[0] = pi as a symbolic (exact) number, an use RealIntervalField(precision) to compute the integer part of x[i].

( 2013-11-06 01:57:03 -0500 )edit

Thank you for the explanation. I see that this question was tagged "confirmed_bug". Is there a place where I can follow what happens to this bug? Any maintainer working on it?

( 2013-11-06 20:50:53 -0500 )edit

I opened a kind of meta-ticket about annoying things abour RIF [here](http://trac.sagemath.org/ticket/15360). I am supposed to write a patch since i already know where are some problems in the source code, but it is likely that i will wait the new git framework to work on it, since such a patch will touch quite a lot of files.

( 2013-11-07 07:02:13 -0500 )edit