[math-fun] clean 16-th order AGM iteration
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: ? t = AGM16(1.0, 0.5) [1.00000000000000, 0.292893218813452] \\ == [a, a-b] [0.853460904507906, 2.94235252084151 E-17] [0.853460904507906, 3.16563558120117 E-273] [0.853460904507906, 1.02028230345421 E-4368] [0.853460904507906, 1.38315649004407 E-69896] [0.853460904507906, 0.E-100009] \\ using a precision of 100,000 digits 0.728395515523453 \\ returned value ? t-agm(a,b) -4.08056715461103 E-100010 \\ check
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
[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.
Order 16 product for hypergeom([+1/2,+1/2],[1],k,N), this one may actually be new: N=300; L=ceil(log(N)/log(16)) \\ need this many steps k='k+O('k^N); K=hypergeom([+1/2,+1/2],[1],k,N); \\ needs hypergeom.inc.gp, else use AGM below \\K=1/agm(1,sqrt(1-k)); k='k+O('k^N); tt=1; { for (i=1, L, skp = (1-k)^(1/4); P = (1 + skp)/2; M = (1 - skp)/2; R = (P^4-M^4)^(1/4); tt *= 4/(P+R)^2; k = (( P - R ) / ( P + R ))^(4) + O('k^N); ); } tt-K \\ == zero Some savings are possible, e.g. as follows: tt=1; { for (i=1, L, kp = (1-k)^(1/2); \\ SQRT skp = kp^(1/2); \\ SQRT P = (1 + skp)/2; R = (skp*(1+kp)/2)^(1/4); \\ MULT & 2 x SQRT m = 1/(P+R); \\ INVERSE tt *= 2*m; \\ MULT k = ((P - R)*m)^4 + O('k^N); \\ 2 x SQUARE ); \\ total: 4 x SQRT, 1 x INVERSE, 2 x MULT, 2 x SQUARE tt *= tt; \\ final square } @RWG: congrats on that 163 thing!
participants (2)
-
Fred lunnon -
Joerg Arndt