Revision history [back]
 | 1 | initial version |
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*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)))) + sqrt(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*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 - 256*(sqrt(2) + 2)*sqrt(1/2)*2^(1/4) - 128*(4*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) - 64*(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))) - 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))) + 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)))) - 16*(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)))) + sqrt(2) + 2)*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))))) - 256*sqrt(2) - 128)
I let you bet the size of BS(20) !!!
So, you have to work numerically from the beginnig. I am currently writing a tutorial about different representations of real numbers in Sage. In the meanwhile, you can have a look at this ask question.
You suggest to "modify the numerical precision in every step". WARNING !!! You can not add precision along the algorithm in such an algorithm, as it can be done for example with algorithms looking for attractove fixed points.
See http://wwwmaths.anu.edu.au/~brent/pub/pub028.html page 11, comment 1 : "Unlike Newton's iteration, the AG mean iteration is not self-correcting. Thus, we cannot start with low precision and increase it, as was possible in Section 2."
So you have to fix your precision at the beginning of your algorithm
(depending on N). Typically of the order of 2^N.
You can use RealField, but if you are not able to guess the loss of
precision due to the roundings, you should use RealIntervalField, that
will provide a guaranty on the accuracy of the result.
sage: def BS_RIF(N):
sage: prec = 2^N
sage: R = RealIntervalField(prec)
sage: a0 = R(1)
sage: b0 = 1/sqrt(R(2))
sage: s0 = 1/R(2)
sage: s = s0
sage: B = 2
sage: a1 = a0
sage: b1 = b0
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
sage: p = BS_RIF(10) ; p
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127373?
sage: p.absolute_diameter()
5.689536106941499008575555960546130516249177454201987114043163738254750611863536965981360580426580804692525828365655879098284825675433813032945912229534680802819610637241265279059667943112160783272264161040682451473408157600065113693301079013467783152091488508497612611025084280206747911787727795216158247186e-303
Which means that i got 303 digits of precision with N=10 iterations. With N = 11, i got 611, then 1227, 2459, and so on, so i can guarantee that the number of accurate
digits of precision is indeed doubled at each iteration.
| | 2 | No.2 Revision |
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*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)))) + sqrt(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*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 - 256*(sqrt(2) + 2)*sqrt(1/2)*2^(1/4) - 128*(4*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) - 64*(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))) - 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))) + 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)))) - 16*(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)))) + sqrt(2) + 2)*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))))) - 256*sqrt(2) - 128)
which has length 3294. I let you bet the size of BS(20) !!!
So, you have to work numerically from the beginnig. I am currently writing a tutorial about different representations of real numbers in Sage. In the meanwhile, you can have a look at this ask question.
You suggest to "modify the numerical precision in every step". WARNING !!! You can not add precision along the algorithm in such an algorithm, as it can be done for example with algorithms looking for attractove fixed points.
See http://wwwmaths.anu.edu.au/~brent/pub/pub028.html page 11, comment 1 : "Unlike Newton's iteration, the AG mean iteration is not self-correcting. Thus, we cannot start with low precision and increase it, as was possible in Section 2."
So you have to fix your precision at the beginning of your algorithm
(depending on N). Typically of the order of 2^N.
You can use RealField, but if you are not able to guess the loss of
precision due to the roundings, you should use RealIntervalField, that
will provide a guaranty on the accuracy of the result.
sage: def BS_RIF(N):
sage: prec = 2^N
sage: R = RealIntervalField(prec)
sage: a0 = R(1)
sage: b0 = 1/sqrt(R(2))
sage: s0 = 1/R(2)
sage: s = s0
sage: B = 2
sage: a1 = a0
sage: b1 = b0
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
sage: p = BS_RIF(10) ; p
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127373?
sage: p.absolute_diameter()
5.689536106941499008575555960546130516249177454201987114043163738254750611863536965981360580426580804692525828365655879098284825675433813032945912229534680802819610637241265279059667943112160783272264161040682451473408157600065113693301079013467783152091488508497612611025084280206747911787727795216158247186e-303
Which means that i got 303 digits of precision with N=10 iterations. With N = 11, i got 611, then 1227, 2459, and so on, so i can guarantee that the number of accurate
digits of precision is indeed doubled at each iteration.
| | 3 | No.3 Revision |
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*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)))) + sqrt(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*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 - 256*(sqrt(2) + 2)*sqrt(1/2)*2^(1/4) - 128*(4*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) - 64*(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))) - 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))) + 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)))) - 16*(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)))) + sqrt(2) + 2)*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))))) - 256*sqrt(2) - 128)
which has length 3294. I let you bet the size of BS(20) !!!
So, you have to work numerically from the beginnig. I am currently writing a tutorial about different representations of real numbers in Sage. In the meanwhile, you can have a look at this ask question.
You suggest to "modify the numerical precision in every step". WARNING !!! You can not add precision along the algorithm in such an algorithm, as it can be done for example with algorithms looking for attractove fixed points.
See http://wwwmaths.anu.edu.au/~brent/pub/pub028.html page 11, comment 1 : "Unlike Newton's iteration, the AG mean iteration is not self-correcting. Thus, we cannot start with low precision and increase it, as was possible in Section 2."
So you have to fix your precision at the beginning of your algorithm
(depending on N). Typically of the order of 2^N.
You can use RealField, but if you are not able to guess the loss of
precision due to the roundings, you should use RealIntervalField, that
will provide a guaranty on the accuracy of the result.
sage: def BS_RIF(N):
sage: prec = 2^N
sage: R = RealIntervalField(prec)
sage: a0 = R(1)
sage: b0 = 1/sqrt(R(2))
sage: s0 = 1/R(2)
sage: s = s0
sage: B = 2
sage: a1 = a0
sage: b1 = b0
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
sage: p = BS_RIF(10) ; p
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127373?
sage: p.absolute_diameter()
5.689536106941499008575555960546130516249177454201987114043163738254750611863536965981360580426580804692525828365655879098284825675433813032945912229534680802819610637241265279059667943112160783272264161040682451473408157600065113693301079013467783152091488508497612611025084280206747911787727795216158247186e-303
Which means that i got 303 an error of at most 5.7e-303, that is i have correct 302 digits of precision with N=10 iterations. With N = 11, i got 611, 610, then 1227, 2459, 1226, 2458, and so on, so i can guarantee that the number of accurate
digits of precision is indeed doubled at each iteration.
| | 4 | No.4 Revision |
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*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)))) + sqrt(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*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 - 256*(sqrt(2) + 2)*sqrt(1/2)*2^(1/4) - 128*(4*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) - 64*(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))) - 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))) + 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)))) - 16*(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)))) + sqrt(2) + 2)*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))))) - 256*sqrt(2) - 128)
which has length 3294. I let you bet the size of BS(20) !!!
So, you have to work numerically from the beginnig. I am currently writing
a tutorial about different representations of real numbers in Sage. In the
meanwhile, Meanwhile, you can have a look at this ask question. to see what are the possibilities.
You suggest to "modify the numerical precision in every step". WARNING !!! You can not add precision along the algorithm in such an algorithm, as it can be done for example with algorithms looking for attractove fixed points.
See http://wwwmaths.anu.edu.au/~brent/pub/pub028.html page 11, comment 1 : "Unlike Newton's iteration, the AG mean iteration is not self-correcting. Thus, we cannot start with low precision and increase it, as was possible in Section 2."
So you have to fix your precision at the beginning of your algorithm
(depending on N). Typically of the order of 2^N.
You can use RealField, but if you are not able to guess the loss of
precision due to the roundings, you should use RealIntervalField, that
will provide a guaranty on the accuracy of the result.
sage: def BS_RIF(N):
sage: prec = 2^N
sage: R = RealIntervalField(prec)
sage: a0 = R(1)
sage: b0 = 1/sqrt(R(2))
sage: s0 = 1/R(2)
sage: s = s0
sage: B = 2
sage: a1 = a0
sage: b1 = b0
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
sage: p = BS_RIF(10) ; p
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127373?
sage: p.absolute_diameter()
5.689536106941499008575555960546130516249177454201987114043163738254750611863536965981360580426580804692525828365655879098284825675433813032945912229534680802819610637241265279059667943112160783272264161040682451473408157600065113693301079013467783152091488508497612611025084280206747911787727795216158247186e-303
Which means that i got an error of at most 5.7e-303, that is i have correct 302 digits of precision with N=10 iterations. With N = 11, i got 610, then 1226, 2458, and so on, so i can guarantee that the number of accurate
digits of precision is indeed doubled at each iteration.
| | 5 | No.5 Revision |
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*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)))) + sqrt(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*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 - 256*(sqrt(2) + 2)*sqrt(1/2)*2^(1/4) - 128*(4*sqrt(1/2)*2^(1/4) + sqrt(2) + 2)*sqrt((sqrt(2) + 2)*sqrt(1/2)*2^(1/4)) - 64*(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))) - 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))) + 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)))) - 16*(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)))) + sqrt(2) + 2)*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))))) - 256*sqrt(2) - 128)
which has length 3294. I let you bet the size of BS(20) !!!
So, you have to work numerically from the beginnig. I am currently writing a tutorial about different representations of real numbers in Sage. Meanwhile, you can have a look at this ask question to see what are the possibilities.
WARNING : You suggest to "modify "modify the numerical precision in every step". WARNING !!!
step".
You can not add precision along the algorithm in such an algorithm, as it
can be done for example with algorithms looking for attractove attractive fixed
points.
points.
See http://wwwmaths.anu.edu.au/~brent/pub/pub028.html page 11, comment 1 :
"Unlike "Unlike Newton's iteration, the AG mean iteration is not self-correcting.
Thus, we cannot start with low precision and increase it, as was possible in
Section 2."2."
So Therefore, you have to fix your precision at the beginning of your algorithm
(depending on N). Typically of the order of 2^N.
You .
For that, you can use RealField, but if you are not able to guess the loss of
precision due to the roundings, you should use RealIntervalField, that
will provide a guaranty on the accuracy of the result.
sage: def BS_RIF(N):
sage: prec = 2^N
sage: R = RealIntervalField(prec)
sage: a0 = R(1)
sage: b0 = 1/sqrt(R(2))
sage: s0 = 1/R(2)
sage: s = s0
sage: B = 2
sage: a1 = a0
sage: b1 = b0
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
sage: p = BS_RIF(10) ; p
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127373?
sage: p.absolute_diameter()
5.689536106941499008575555960546130516249177454201987114043163738254750611863536965981360580426580804692525828365655879098284825675433813032945912229534680802819610637241265279059667943112160783272264161040682451473408157600065113693301079013467783152091488508497612611025084280206747911787727795216158247186e-303
Which means that i got an error of at most 5.7e-303, that is i have correct 302 digits of precision with N=10 iterations. With N = 11, i got 610, then 1226, 2458, and so on, so i can guarantee that the number of accurate
digits of precision is indeed doubled at each iteration.
