A simpler iteration for sin(theta/3) comes from "solving" the triplication equation: sin(3x) = 4 sin(x)^3 - 3 sin(x) sin(x) = (sin(3x) + 4 sin(x)^3) / 3 The iteration is to plug the current guess for sin(x) into the RHS, giving an improved guess for sin(x) on the LHS. Unfortunately, the convergence is inferior to RWGs scheme -- merely linear. If we assume out current guess is sin(x)+eps, the new guess is sin(x) + eps * (4 sin(x)^2) + ... . This gives OKish convergence for small sines. If sin(3x) = 1/2, and we take the first guess for sin(x) = sin(3x)/3 (i.e., a 0th guess of sin(x)=0), the next iteration is sin(x) = 14/81 = .1728... . The error has dropped from .0070 to .0008. Successive steps will get another factor of 8 each time. Rich -----Original Message----- From: math-fun-bounces+rschroe=sandia.gov@mailman.xmission.com [mailto:math-fun-bounces+rschroe=sandia.gov@mailman.xmission.com] On Behalf Of R. William Gosper Sent: Thursday, July 13, 2006 4:46 AM To: math-fun@mailman.xmission.com Subject: Re: [math-fun] Calculating pi Archimedes' pi method with sin(pi/6), sin(pi/12), sin(pi/24), ... requires a sqrt for each new sine. If instead, we sought sin(pi/6), sin(pi/18), sin(pi/54), ..., we'd need a cbrt, but each iteration provides an excellent starting value for a Newton iteration for the next cbrt. The triple angle formula is (c102) (sin(3*x),%% = decosify(trigexpand(%%))) 2 3 (d102) sin(3 x) = 3 sin(x) (1 - sin (x)) - sin (x) We can't solve this for sin(x) without a nasty complex cubic, which explains why there's no nice expression for sin(pi/9), etc. But letting S := sin(3x), s := sin(x), the Newton step toward s given S is (c197) newt : (8*s^3-\s)/(3*(4*s^2-1)) 3 8 s - S (d197) ------------ 2 3 (4 s - 1) Starting sith S = sin(pi/6), (c198) subst(1/2,\s,newt) 3 1 8 s - - 2 (d198) ------------ 2 3 (4 s - 1) we seek s = sin(pi/18), using sin(pi/6)/3 as our first guess: (c199) subst(1/6.0d0,s,%) (d199) 0.17361111111111d0 plugging back in (c200) subst(%,s,d198) (d200) 0.17364817658193d0 once more (c201) subst(%,s,d198) (d201) 0.17364817766693d0 This has converged, since the agreement doubles each time. So now we go for sin(pi/54) (c202) subst(%,\s,newt) 3 8 s - 0.17364817766693d0 (d202) ------------------------- 2 3 (4 s - 1) starting with sin(pi/18)/3: (c203) subst(d201/3,S,%) (d203) 0.05814481276438d0 (c204) subst(%,s,d202) (d204) 0.05814482891048d0 Already finished! Go for sin(pi/162). (c205) subst(%,\s,newt) 3 8 s - 0.05814482891048d0 (d205) ------------------------- 2 3 (4 s - 1) (c206) subst(d204/3,s,%); (d206) 0.01939133176448d0 [...] (c211) subst(%,\s,newt) 3 8 s - 0.00646413739654d0 (d211) ------------------------- 2 3 (4 s - 1) (c212) subst(d210/3,s,%) (d212) 0.00215472580425d0 Now we don't even need to iterate! (c213) subst(%,s,d211) (d213) 0.00215472580425d0 This was our fifth thirding of pi/6: (c214) (sin(%pi/6/3^5),%% = dfloat(%%)) %pi (d214) sin(----) = 0.00215472580425d0 1458 (c215) 1458*%; %pi (d215) 1458 sin(----) = 3.14159022259953d0 1458 Here's a different recursion, based on tangents: On 18 Apr 05 I sent
[...] another one of those annoying little pi iterations: (c124) tn2(x):=if equal(12,12-(x:x/2)^2) then x else ((tn2(x)-x)/2,2*x-4/(%%-1/%%))
x 2 (d124) tn2(x) := if equal(12, 12 - (x : -) ) then x 2 tn2(x) - x 4 else (----------, 2 x - -------) 2 1 %% - -- %% The symbol %% is a pronoun meaning the most recent result in the computation, in this case, (tn2(x)-x)/2. tn2(x) computes x-2*tan(x/4) recursively via the tangent double angle formula. Iterating this function converges quadratically to pi: (c125) tn2(3.0) (d125) 1.1368 (c126) tn2(2+%) (d126) 1.14158 (c127) tn2(2+dfloat(%)) (d127) 1.14159265355336d0 (c128) (bfprecision:23,tn2(2+bfloat(%))) (d128) 1.1415926535897932384623b0 Finally, in late '02 I sent a method based on an amazingly simple way to compute the vector [sin(x), sin(x/3), ..., epsilon], where sin(epsilon) ~ epsilon, again using the triple angle formula: (c6) sn3(x) := if equal(6,x+6) then [x] else (sn3(x/3),cons(%%[1]*(3-4*%%[1]^2),%%))$ Then with the help of magic rational vectors we can get pi with any desired rate of convergence. E.g., for 9th order (enneadic?): (c134) pi27(x):=x+[-1/12579840,1/5120,-243/5120,531441/465920] . part(sn3(27*x),1..4)$ (c137) pi27(22/7.b0); (d137) 3.14159265358979323846264339539b0 (In case anybody missed it, pi27(d137) would give *nine times* as many correct digits.) The magic vector is 1 (d139) ------------------------------ 1 1 n - k (-; -) 3 (9; 9) 9 9 k - 1 n - k n = 4, k=1..4, and (a;q)_n is the q-Pochhammer symbol (1-a)(1-aq)...(1-aq^(n-1)). Of course, if you have a sin function simply iterating x <- x + sin x will triple the number of digits of pi each time, as Gene Salamin pointed out many years ago in HAKMEM. --rwg _______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun