Ask Your Question

Continued fraction of pi by hand

asked 2013-11-05 21:41:54 -0500

anonymous user


updated 2013-11-05 23:53:38 -0500

tmonteil gravatar image

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

for i in range(n):
    a[i] = int(x[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?


edit retag flag offensive close merge delete

3 answers

Sort by » oldest newest most voted

answered 2013-11-05 23:53:10 -0500

tmonteil gravatar image

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)
sage: int(Field(x).upper()) == int(Field(x).lower())
sage: int(Field(x).upper())

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

edit flag offensive delete link more



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

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

answered 2014-11-21 23:45:08 -0500

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)

for i in xrange(n):
    a[i] = int(x[i])
edit flag offensive delete link more

answered 2013-11-06 00:38:43 -0500

Marco Caliari gravatar image

Thanks for the quick answer. I extended the precision

for i in range(n):
    a[i] = int(x[i])

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

edit flag offensive delete link 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]`.

tmonteil gravatar imagetmonteil ( 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?

Marco Caliari gravatar imageMarco Caliari ( 2013-11-06 20:50:53 -0500 )edit

I opened a kind of meta-ticket about annoying things abour RIF [here]( 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.

tmonteil gravatar imagetmonteil ( 2013-11-07 07:02:13 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools


Asked: 2013-11-05 21:41:54 -0500

Seen: 251 times

Last updated: Nov 21 '14