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 )
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)
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 ...Asked: 2013-06-25 06:24:38 +0200
Seen: 382 times
Last updated: Jun 26 '13
Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
