Ask Your Question
1

computing the digist of Pi

asked 2013-06-24 23:24:38 -0500

czsan gravatar image

updated 2015-08-31 08:43:02 -0500

tmonteil gravatar image

I tried to code the Salamin-Brent iteration. I calculated with exact numbers and the procedure returned the numerical value with the .n() function. The iteration was very very slow, for 20 iterations it is approximately 900 secs. (In Maple the same code is 0.06 secs). Can I modify the numerical precision in every step with some method? (Because in every step the number of correct digits is doubling )

edit retag flag offensive close merge delete

2 answers

Sort by ยป oldest newest most voted
2

answered 2013-06-25 05:01:47 -0500

tmonteil gravatar image

updated 2013-06-26 10:54:03 -0500

So, your code is:

sage: def BS(N): 
sage:     a0 = 1
sage:     b0 = 1/sqrt(2)
sage:     s0 = 1/2
sage:     B = 2 
sage:     a1 = a0
sage:     b1 = b0
sage:     s = s0 
sage:     for k in range(N+1):
sage:         a = (a1+b1)/2
sage:         b2 = a1 * b1
sage:         b = sqrt(b2)
sage:         c = a^2-b2
sage:         s -= B * c 
sage:         p = 2* a^2/s 
sage:         a1 = a
sage:         b1 = b
sage:         s1 = s
sage:         B = 2*B 
sage:     return p.n(digits=2^N)

When you write b0 = 1/sqrt(2), you define an element of the Symbolic Ring:

sage: b0 = 1/sqrt(2)
sage: b0.parent()
Symbolic Ring

Which means that all compurations are done in a symboloc way, there is no numerical approximation at all there. So, with square roots, products and sums, Sage will deal with bigger and bigger formulas, which explains your timings.

For example, if i replace return p.n(digits=2^N) by return p in your code, i got:

sage: BS(5)
-1/32*(4*sqrt(1/2)*2^(1/4) + 4*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) + 4*sqrt((4*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4))) + 4*sqrt((4*sqrt(1/2)*2^(1/4) + 4*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) + sqrt(2) + 2)*sqrt((4*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)))) + 4*sqrt((4*sqrt(1/2)*2^(1/4) + 4*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) + 4*sqrt((4*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4))) + sqrt(2) + 2)*sqrt((4*sqrt(1/2)*2^(1/4) + 4*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) + sqrt(2) + 2)*sqrt((4*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4))))) + sqrt(2) + 2)^2/(32*(sqrt(2) + 2)^2 + 16*(4*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)^2 + 8*(4*sqrt(1/2)*2^(1/4) + 4*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) + sqrt(2) + 2)^2 + 4*(4*sqrt(1/2)*2^(1/4) + 4*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) + 4*sqrt((4*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4))) + sqrt(2) + 2)^2 + 2*(4*sqrt(1/2)*2^(1/4) + 4*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) + 4*sqrt((4*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4))) + 4*sqrt((4*sqrt(1/2)*2^(1/4) + 4*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) + sqrt(2) + 2)*sqrt((4 ...
(more)
edit flag offensive delete link more

Comments

Thank's, it's very nice! Mainly I'm in trouble with variables after defining a ring.

czsan gravatar imageczsan ( 2013-06-25 07:46:08 -0500 )edit

Looking forward to the real number tutorial!

Eviatar Bach gravatar imageEviatar Bach ( 2013-06-25 09:40:08 -0500 )edit
0

answered 2013-06-25 02:26:34 -0500

tmonteil gravatar image

Could you provide your code ? My guess is that you compute your iteration in the symbolic ring, and then use .n() method to get your numbers at the end, but i am not sure.

Note that, if you want to get X digits of precision of an exact number a, you can do:

sage: a.n(digits=X)
edit flag offensive delete link more

Comments

def BS(N): a0 = 1; b0 = 1/sqrt(2); s0 = 1/2; B = 2 a1 = a0; b1 = b0; s = s0 for k in range(N+1): a = (a1+b1)/2; b2 = a1*b1; b = sqrt(b2) c = a^2-b2 s -= B*c p = 2*a^2/s a1 = a; b1 = b; s1 = s; B = 2*B return p.n(digits=2^n)

czsan gravatar imageczsan ( 2013-06-25 02:35:22 -0500 )edit

def BS(N): a0 = 1; b0 = 1/sqrt(2); s0 = 1/2; B = 2 a1 = a0; b1 = b0; s = s0 for k in range(N+1): a = (a1+b1)/2; b2 = a1*b1; b = sqrt(b2) c = a^2-b2 s -= B*c p = 2*a^2/s a1 = a; b1 = b; s1 = s; B = 2*B return p.n(digits=2^n)

czsan gravatar imageczsan ( 2013-06-25 02:35:22 -0500 )edit

it would be enough once... the last n is N, naturally so, i would like to compute in every iteration with 2^k digits

czsan gravatar imageczsan ( 2013-06-25 02:37:48 -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-06-24 23:24:38 -0500

Seen: 117 times

Last updated: Jun 26 '13