Ask Your Question
1

Continued fraction of pi by hand

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

anonymous user

Anonymous

updated 2017-04-19 17:32:02 -0500

tmonteil gravatar image

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 flag offensive close merge delete

3 answers

Sort by ยป oldest newest most voted
1

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

edit flag offensive delete link more

Comments

1

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
0

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)

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]),
edit flag offensive delete link more
0

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

Marco Caliari gravatar image

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?

edit flag offensive delete link more

Comments

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](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.

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

Stats

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

Seen: 281 times

Last updated: Nov 21 '14