[sorry for delay b/c I was cut off my mail] * Fred lunnon <fred.lunnon@gmail.com> [Jul 16. 2009 16:07]:
On 7/15/09, Joerg Arndt <arndt@jjj.de> wrote:
AGM16(a,b)= \\ return AGM(a,b) { local(P,M,R); a = sqrt(a); b = sqrt(b); while ( abs(a-b) > eps, \\ eps == 10^(-precision) P=(a+b)/2; M=(a-b)/2; R=(P^4-M^4)^(1/4); a=(P+R)/2; b=(P*R*(P^2+R^2)/2)^(1/4); ); return( a^2 ); }
The rate of convergence is _very_ impressive:
Looking at the code though, I can't help wondering whether it might possibly be massaged into a pair of fourth-order iterations laid end-to-end [no-one would be in the least surprised, as Dorothy Parker might acerbically have observed.] WFL
This is true indeed as I found in the hour after sending: AGM16(a,b)= \\ return AGM(a,b) { local(P,M,R, c); a = sqrt(a); b = sqrt(b); while ( abs(a-b) > eps, \\ eps == 10^(-precision) print( [a, a-b] ); P=(a+b)/2; M=(a-b)/2; R=(P^4-M^4)^(1/4); \\ R=(a*b*(a^2+b^2)/2)^(1/4); a=(P+R)/2; c=(P-R)/2; \\ b=(P*R*(P^2+R^2)/2)^(1/4); b=(a^4-c^4)^(1/4); ); print( [a, a-b] ); return( a^2 ); } /* ----- */ Tried shortly (pencil & paper) to combine the steps but things get messy (and Schoenhage's AGM is unbeatable anyway). Maybe the corresponding rule k(q) --> k(q^16) is more interesting: sqrt(k(q^16)) = ( P - (P^4-M^4)^(1/4) ) / ( P + (P^4-M^4)^(1/4) ) where P:=1-sqrt(k'(q)) and M = 1 + sqrt(k'(q)) also sqrt(k'(q^16)) = (8*P*R*(P^2+R^2))^(1/4) / (P+R) where R:=(P^4-M^4)^(1/4) Same waste-of-time session, write E_k for E(q^k): +-+ 4 2 6 4 \|q E E E 8 2 4 k = ---------------, 12 4 8 E + 4 q E E 4 2 8 12 4 8 E -4 q E E 4 2 8 k' = ------------- 12 4 8 E +4 q E E 4 2 8 The last one I like. From that one follows: [ 12 [ 12 4 8 ] ] | E | E -4 q E E | | | 2 | 4 2 8 | |1/8 E_1 = | --- | ------------- | | | 4 | 12 4 8 | | | E | E +4 q E E | | [ 4 [ 4 2 8 ] ] Also, Tk for Theta_K: T3(q) = 1 + T2(q^4) + T2(q^16) + T2(q^64) + T2(q^256) + ... T4(q) = 1 - T2(q^4) + T2(q^16) + T2(q^64) + T2(q^256) + ... have not seen these, surely known. These ones are easy, but I missed them so far: T3(q^16) + T2(q^16) = 1/2 * ( T3(q) + T4(q) ) T2(q)^2 = 2*T3(q^2)*T2(q^2) EOT from the elliptic nonsense department.