A friend of mine recently asked me what pi is all about. This was rather remarkable, as he's failing out of high school, and people treat him like a failure. Turns out he's fascinated by physics and math. Anyway. So I described pi in terms of the experimental method of drawing circles and noticing that their circumference always seems a little more than 3 times their diameter. That led to a great discussion about experimental vs. theoretical science, and eventually he asked me how we went from experimental methods of calculating pi to computational approaches. I checked the wikipedia, but its approaches all seem to assume advanced math to prove that the approximations describe pi, and he hasn't gotten past algebra and geometry. At this point, the simplest approach I can think of is to write a program to choose random points in a circumscribed circle, then calculate the proportion of the points which fall within it. But that still has a heuristic element to it. Are there any plain arguments for the deterministic approximations? (For that matter, how about doing the same thing for the trig functions? I took the same approach for sin/cos/tan, talking about building a table of angle/ratio measurements, but I had no idea how to describe computational approaches.) -J
Surely best to follow Archimedes and calculate the perimeters of inscribed and circumscribed regular 6-, 12-, 24-, ...-gons? R. On Wed, 12 Jul 2006, Jason Holt wrote:
A friend of mine recently asked me what pi is all about. This was rather remarkable, as he's failing out of high school, and people treat him like a failure. Turns out he's fascinated by physics and math. Anyway. So I described pi in terms of the experimental method of drawing circles and noticing that their circumference always seems a little more than 3 times their diameter. That led to a great discussion about experimental vs. theoretical science, and eventually he asked me how we went from experimental methods of calculating pi to computational approaches.
I checked the wikipedia, but its approaches all seem to assume advanced math to prove that the approximations describe pi, and he hasn't gotten past algebra and geometry. At this point, the simplest approach I can think of is to write a program to choose random points in a circumscribed circle, then calculate the proportion of the points which fall within it. But that still has a heuristic element to it. Are there any plain arguments for the deterministic approximations?
(For that matter, how about doing the same thing for the trig functions? I took the same approach for sin/cos/tan, talking about building a table of angle/ratio measurements, but I had no idea how to describe computational approaches.)
-J
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
From: "Richard Guy" <rkg@cpsc.ucalgary.ca>
Surely best to follow Archimedes and calculate the perimeters of inscribed and circumscribed regular 6-, 12-, 24-, ...-gons? R.
Another approach is using a polynomial approximation of sin. Linear approximation sin(x) = x is not very good. For x=pi/6 it gives 1/2=pi/6 and pi=3. However, cubic approximation sin(x) = x + Ax^3 is already pretty good (quadratic term is not present because sin is odd). From sin(pi/6)=1/2 and sin(pi/4)=sqrt(2)/2=.707 we get a system of linear equations pi/6 + A pi^3/6^3 = 0.5 pi/4 + A pi^3/4^3 = 0.707 Multyplying first equation by 6^3 and second - by 4^3, we get 36 pi + A pi^3 = 108 16 pi + A pi^3 = 45.2 Subtracting, 20 pi = 62.8. and pi = 3.14 Alec Mihailovs
On Wed, 12 Jul 2006, Richard Guy wrote:
Surely best to follow Archimedes and calculate the perimeters of inscribed and circumscribed regular 6-, 12-, 24-, ...-gons? R.
Working that out computationally seems to require sines, so that the problem of calculating pi reduces to the problem of calculating sines, yes? -J
On 7/12/06, Jason Holt <jason@lunkwill.org> wrote:
Working that out computationally seems to require sines, so that the problem of calculating pi reduces to the problem of calculating sines, yes?
-J
No, just the pythagorean theorem. Start with an equilateral triangle; bisect it to get triangles of length 1,1/2, sqrt(3)/2. Then you attach two more triangles below it with lengths 1/2, 1-sqrt(3)/2, sqrt(5-2sqrt(3)/4) to get the sides of the dodecagon. Something vaguely like this: /|\ / | \ --+-- \|/ -- Mike Stay metaweta@gmail.com http://math.ucr.edu/~mike
On 7/12/06, Mike Stay <mike@math.ucr.edu> wrote:
sqrt(5-2sqrt(3)/4) to get the sides of the dodecagon. Oops. I meant sqrt((5-2sqrt(3))/4) -- Mike Stay metaweta@gmail.com http://math.ucr.edu/~mike
Yes, but no need to know what a sine is. E.g. for the 6-gon, we only need to know that the base of an equilateral triangle is equal to the other two sides. So half of it gives `sin pi/6 = 1/2', but this statement adds nothing to the understanding of a beginner. The essentials are Pythagoras's theorem and the angle bisector theorem. Suppose OAB is a triangle, right-angled at B and angle BOA = 30 degrees. If OB = 1, then AB = 1/2 and the inscribed 6-gon has perimeter 6. If OA = 1, then AB = 2 / root 3, and the circumscribed 6-gon has perimeter 4 root 3 To get the 12-gon, we need to bisect the angle. For the circumscribed 12-gon it's `tan pi/12' we want, and for the inscribed one, `sin pi/12', but the names and notation don't help. Bisect the angle BOA and let the bisector meet AB in C. Then AC : CB = OA : OB = root 3 : 2 and the (`tangent') length we want is 2 - root 3, twenty-four of which give approximately 6.43, while the inscribed (`sine') length is (root(2 - root 3)) / 2, twenty-four of which give approximately 6.21. I don't recall how many stages Archimedes went through, but he used (better and better) rational approximations to root 3, and gave 3 1/7 > pi > 3 10/71 whose mean is 3.14185. No need to know any trigonometry or analysis. R. On Wed, 12 Jul 2006, Jason Holt wrote:
On Wed, 12 Jul 2006, Richard Guy wrote:
Surely best to follow Archimedes and calculate the perimeters of inscribed and circumscribed regular 6-, 12-, 24-, ...-gons? R.
Working that out computationally seems to require sines, so that the problem of calculating pi reduces to the problem of calculating sines, yes?
-J
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
Thanks for all the suggestions! I also just realized that my nondeterministic solution of choosing random points and finding out whether they're in the circle works just as well with deterministic uniform point selection. In other words, counting cells inside a circle drawn on graph paper instead of counting darts thrown at a dart board. Let's see if I can write it as a formula instead of an algorithm: n n 4 ___ ___ --- \ \ pi ~ n*n * / / 1 if sqrt(x*x+y*y) < n --- --- x=0 y=0 Does that look about right? I'm surprised I didn't run across anything like that. Maybe I should add it to the wikipedia article for the benefit of readers with more elementary math backgrounds. I suppose it's not quite as nice as a series, since you have to decide when to stop before you begin. Perhaps we could get a run-until-you-get-board series if we did something tricky like progressively subdividing the coordinates rather than drawing bigger and bigger circles. I'm also curious whether a similar algorithm could be constructed for approximating sines. Maybe something doing recursive angle bisections as was suggested for the Archimedes approach? -J
On Thu, 13 Jul 2006, Jason Holt wrote:
n n 4 ___ ___ --- \ \ pi ~ n*n * / / 1 if sqrt(x*x+y*y) < n --- --- x=0 y=0
I just noticed that this approach is sort of an intuitive way of looking at the discrete form of the integral from 0..1 of sqrt(1-x^2) dx. The discrete version could also be computed as: n 4 ___ --- * \ pi ~ n*n / sqrt(n*n-x*x) --- x=0 Regarding the subdivision approach to avoid having to choose n at the outset, perhaps we could use a function (call it subdiv()) that, given n, outputs the nth element of the breadth-first traversal of this tree: 0.5 0.25 0.75 0.125 0.375 0.625 0.875 ... Huh. That's the sizes in my SAE socket wrench set: 1/2, 1/4, 3/4, 1/8, 3/8, 5/8, 7/8 [1/16, 3/16, 5/16, 7/16, 9/16, 11/16, 13/16, 15/16] ... So, it's something like... subdiv(n): (1+2*n) mod 2^int(1+log2(n)) -------------------------- 2^int(1+log2(n)) And we could call order(n)=2^int(1+log2(n)). This sort of function must have been invented before; anybody know what it's called? Our series for pi then becomes: n ___ \ order(n)*pi ~ 4* / sqrt(1-subdiv(i)^2) --- i=0 So that's closer, anyway: run as long as you want, then divide the result by order(n) at the end to get your approximation of pi. Of course, it works much better if n=2^k-1 for some k. -J
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
Beginners are often startled by (c64) (-1/2)! (d64) sqrt(%pi) but if we are to have (n+1/2)! = (n+1/2)(n-1/2)! for all n, then (c86) (n+1/2)! = h*product(k-1/2,k,1,n+1) n + 1 /===\ 1 | | 1 (d86) (n + -)! = h | | (k - -) 2 | | 2 k = 1 and sqrt(pi) is the only value of h with 1 (n + -)! n! 2 (n + 1)! (d68) -------- < -------- < -------- 1 n! 1 (n - -)! (n + -)! 2 2 as n -> oo. E.g., using the products for (n-1/2)!,...,(n+1)!, (c77) [h*prod(k-1/2,k,1,n),prod(k,k,1,n),h*prod(k-1/2,k,1,n+1),prod(k,k,1,n+1)] n n n + 1 n + 1 /===\ /===\ /===\ /===\ | | 1 | | | | 1 | | (d77) [h | | (k - -), | | k, h | | (k - -), | | k] | | 2 | | | | 2 | | k = 1 k = 1 k = 1 k = 1 the ratios in (d68) become (c78) rest(%)/rest(%,-1) n n + 1 n + 1 /===\ /===\ /===\ | | | | 1 | | | | k h | | (k - -) | | k | | | | 2 | | k = 1 k = 1 k = 1 (d78) [---------------, ---------------, ---------------] n n n + 1 /===\ /===\ /===\ | | 1 | | | | 1 h | | (k - -) | | k h | | (k - -) | | 2 | | | | 2 k = 1 k = 1 k = 1 (c81) sfloat(sqrt(%pi)) (d81) 1.77245 Choosing h too small, n = 22 is large enough to fail: (c79) apply_nouns(subst([n = 22,h = 1.5d0],d78)) (d79) [5.5739412602892d0, 4.03664103177732d0, 5.69780662162896d0] Increasing h not quite enough, (c80) apply_nouns(subst([n = 22,h = 1.75d0],d78)) (d80) [4.77766393739075d0, 4.70941453707354d0, 4.88383424711054d0] With h too big, (c82) apply_nouns(subst([n = 22,h = 1.95d0],d78)) (d82) [4.28764712329939d0, 5.24763334131052d0, 4.38292817048382d0] With h very close, (c83) apply_nouns(subst([n = 22,h = 1.77d0],d78)) (d83) [4.72367903414339d0, 4.76323641749724d0, 4.82864967934658d0] and we'd need to increase n>22 to see which way h is off. This process produces the constant 1.7724..., but doesn't tell us that it's sqrt(pi). But the above arguments lead to the generalization (c4) (product(n*(1+1/n)^z/(n+z),n,1,inf),%% = closedform(%%)) inf 1 z /===\ (- + 1) n | | n (d4) | | ---------- = z! | | z + n n = 1 and for z = -1/2, (c91) subst(-1/2,z,%^2) inf /===\ | | n 2 2 (d91) ( | | -------------------) = (-1/2)! | | 1 1 n = 1 sqrt(- + 1) (n - -) n 2 (after squaring) we can get Wallis's product (c92) factor(intoprod(%)) inf /===\ 3 | | 4 n (d92) | | ------------------ = %pi | | 2 n = 1 (n + 1) (2 n - 1) (c93) factor(prodfudge(%,n/(n-1/2)))/2 inf /===\ 2 | | 4 n %pi (d93) | | ------------------- = --- | | (2 n - 1) (2 n + 1) 2 n = 1 which he presumably got from Euler's infinite product for sin(x)/x. But, as Rich and I mentioned a couple of weeks ago, factorial can be q-extended, inf /===\ i | | 1 - q | | ---------- | | z + i i = 1 1 - q (d97) QF(z, q) = ---------------- z (1 - q) (which becomes (d4) when q -> 1) and when z = -1/2 leads to a pi-like function of q: pi(q) := q^(1/4) QF(-1/2,q^2)^2 . E.g., the q-Wallis product becomes (c1) pi(q)/(1+q) =q^(1/4)*product((1-q^(2*n))/(1-q^(2*n-1)).(1-q^(2*n))/(1-q^(2*n+1)),n,1,inf) inf /===\ 2 n 2 n pi(q) 1/4 | | 1 - q 1 - q (d1) ----- = q | | (------------ . ------------). q + 1 | | 2 n - 1 2 n + 1 n = 1 1 - q 1 - q pi(q) can also be written as a rapidly convergent series: (c3) q^(1/4)*(1-q^2)*sum(q^binomial(n,2),n,1,inf)^2 inf (n - 1) n ==== --------- 1/4 2 \ 2 2 (d3) q (1 - q ) ( > q ) / ==== n = 1 Comparing the Taylor expansions, (c4) taylor(% = (1+q)*rhs(d2),q,0,11) 1/4 5/4 17/4 21/4 25/4 29/4 33/4 (d4)/T/ q + 2 q + q - 2 q + q + 2 q - 3 q 41/4 1/4 5/4 17/4 21/4 25/4 29/4 + 2 q + . . . = q + 2 q + q - 2 q + q + 2 q 33/4 41/4 - 3 q + 2 q + . . . pi(q) can also be written as etas: (c19) pi(q) = (1-q^2)*eta(q^2)^4/eta(q)^2 2 4 2 (1 - q ) eta (q ) (d19) pi(q) = -----------------. 2 eta (q) Recall a week or two back our "least abstruse" eta relation: (c99) 16*eta(q)^8*eta(q^4)^16+eta(q)^16*eta(q^4)^8=eta(q^2)^24 8 16 4 16 8 4 24 2 (d99) 16 eta (q) eta (q ) + eta (q) eta (q ) = eta (q ) This leads to a pi(q) relation: (c5) pi(q) = sqrt(pi(q^2)*(4*pi(q^4)/(q^4+1)+(q^4+1)*pi(q^2)^2/pi(q^4)))/(q^2+1) 4 4 2 2 2 4 pi(q ) (q + 1) pi (q ) sqrt(pi(q ) (-------- + ----------------)) 4 4 q + 1 pi(q ) (d5) pi(q) = ------------------------------------------ 2 q + 1 As Gene Salamin pointed out to me, we can now initialize q to a small number, quickly evaluate the series (d3) for q^2 and q^4, and then use (d5) for q, sqrt q, q^(1/4), ..., 1-epsilon, and pi(1-epsilon) ~~ pi. Let's see how well: To use the recurrence (d5), define a[n] := pi(2^-(2^(4-n))), b[n] := 1+2^-(2^(4-n)): (c23) (bfprecision:25,bfloat(subst(2^-8,q,d16))) (d23) pi(3.90625b-3) = 2.519531250577538086392374b-1 (c24) bfloat(subst(2^-16,q,d16)) (d24) pi(1.52587890625b-5) = 6.250190734863281250338803b-2 (c25) (a[1] : rhs(d23),a[0] : rhs(%), b[0] : 1.0b0+2^-16) (d25) 1.0000152587890625b0 (c26) b[n] := 1+sqrt(b[n-1]-1) (d26) b := 1 + sqrt(b - 1) n n - 1 Then (d5) becomes (c27) a[n] := sqrt(a[n-1]*(b[n-2]*a[n-1]^2/a[n-2]+4*a[n-2]/b[n-2]))/b[n-1] 2 b a 4 a n - 2 n - 1 n - 2 sqrt(a (------------- + --------)) n - 1 a b n - 2 n - 2 (d27) a := --------------------------------------- n b n - 1 and now let's try for pi: (c29) a[2] (d29) 5.625067088994892759372818b-1 (c30) a[3] (d30) 1.06226912564319452424773b0 (c31) a[4] (d31) 1.699635053182282093917339b0 (c32) a[5] (d32) 2.26618007091359690481384b0 . . . (c48) a[21] (d48) 3.141576039984956345244343b0 (c49) a[22] (d49) 3.141584346772731862424423b0 Pretty disappointing! But now look what happens if we divide a[n] by 2^(n-4) (1-2^(5-n)): (c105) makelist(a[n+4]/2^n/(1-2.0b0^-2^(1-n)),n,-4,2) (d105) [1.000030517810962749189647b0, 2.015655756928048047972053b0, 2.258850470247360857097163b0, 2.26617413470548165172849b0, 2.266180070909709458556452b0, 2.26618007091359690481384b0, 2.266180070913596904813842b0] Quadratic convergence! But to what? (c106) %*2*log(2.0b0) (d106) [1.386336667789142000816372b0, 2.794292209788197907696069b0, 3.1314316695169296120584b0, 3.141584424257956628976669b0, 3.141592653584404093636867b0, 3.141592653589793238462641b0, 3.141592653589793238462643b0] I.e., pi/(2 ln 2). --rwg PS, speaking of pi/(2 ln 2) the first identity in http://www.tweedledum.com/rwg/idents.htm is equivalent to %pi limit 1/(n log(---) - log(atan(n)) - log(atan(n + 1))- . . . n -> inf 2 %pi - log(atan(2 n - 1))) = -------- 2 log(2)
participants (5)
-
Alec Mihailovs -
Jason Holt -
Mike Stay -
R. William Gosper -
Richard Guy