So if we don't have log2 handy, we need to do this twice with different starting parameters, which will ultimately yield both log2 and pi? Rich ----- Quoting Bill Gosper <billgosper@gmail.com>:
I'm pretty sure I sent this before, but since it took me days to dig it back up, here it is again, that quadratically convergent π iteration using "Newton's method for √1". (= Coth double angle formula.) There's also a √ in the loop, but you need only four of them for 34 digits. You also need ln 2 and ten assorted powers of 2:
In[323]:= sum[k_, m_] := Sum[HoldForm @@ {1/2^(k*n*(n - 1)/2)}, {n, m}]^2*k/2^(k/4)
Note these sums represent trivial effort--just sprinkling m 1s into a string of 0s.
In[326]:= sum[4, 7]/sum[8, 5]
Out[326]= 1 1 1 1 Power[-------------------------- + ------------------- + ------------- + -------- + 19342813113834066795298816 1152921504606846976 1099511627776 16777216
1 1 ---- + -- + 1, 2] / 4096 16
1 1 1 1 2 (------------------------- + --------------- + -------- + --- + 1) 1208925819614629174706176 281474976710656 16777216 256
IIn[363]:= N[ReleaseHold[%326], 34]
Out[363]= 1.120652900418845687263888977941327
In[364]:= Sqrt[(% + 1/%)/2]
Out[364]= 1.0032422086165441482848802264099015
In[365]:= Sqrt[(% + 1/%)/2]
Out[365]= 1.0000026194828265385324048741581395
In[366]:= Sqrt[(% + 1/%)/2]
Out[366]= 1.0000000000017154180761229126907979
In[367]:= Sqrt[(% + 1/%)/2]
Out[367]= 1.00000000000000000000000073566479397
In[368]:= %*%%*%%%*%%%%*ReleaseHold[Numerator[%326]]*Log[16]
Out[368]= 3.141592653589793238462643383279502
With three more powers of 2, but no more √ s, we get ~49 digits: In[381]:= sum[4, 9]/sum[8, 6]
In[382]:= N[ReleaseHold[%], 50]
In[383;;386]:= Sqrt[(% + 1/%)/2]
In[387]:= π-Times @@ (Out /@ Range[383, 386])*ReleaseHold[Numerator[%381]]*Log[16]
Out[387]= 4.*10^-49
But preoccupation with minimizing √ s is anachronistic. They should be no costlier than division. So keeping the sum lengths 7 and 5, but doubling k, we get ~67 digits with six √ s: In[370]:= sum[8, 7]/sum[16, 5]
In[371]:= N[ReleaseHold[%], 69]
Out[371]= \2.01559424540388920698155253866341100838034807969448027057273589434627
In[372;;377]:= Sqrt[(% + 1/%)/2]
Out[377]= 1.0000000000000000000000000000000000000000000000001353006722721156829769
In[379]:= Times @@ (Out /@ Range[372, 377])*ReleaseHold[Numerator[%370]]*Log[4]
Out[379]= 3.14159265358979323846264338327950288419716939937510582097494459230751
In[380]:= % - π
Out[380]= -3.1*10^-67 --rwg _______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun