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!