Re: [math-fun] Calculating pi
Jason writes: << . . . 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? . . .
Yes. In fact one of the old classic math books (Geom. & the Imagination? What is Mathematics? Math. & the Imagination?) uses this observation to obtain the series pi/4 = 1 - 1/3 + 1/5 - 1/7 + . . . . --Dan
On Wed, 12 Jul 2006, dasimov@earthlink.net wrote:
Jason writes:
<< . . . 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? . . .
Yes. In fact one of the old classic math books (Geom. & the Imagination? What is Mathematics? Math. & the Imagination?) uses this observation to obtain the series pi/4 = 1 - 1/3 + 1/5 - 1/7 + . . . .
Fascinating! Do you remember how they got from one to the other? -J
On Thu, 13 Jul 2006, Jason Holt wrote:
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
We coded up the algorithm and ran it for varying values of n. Here's the result for n=3,000,000, which took 23 hours on our Athlon64 3000+: ... area = 7068586468447 pi=3.141593985976 (pi/realpi=1.000000424112) Command exited with non-zero status 45 79435.25user 490.80system 23:26:48elapsed Not the most quickly converging series, but it did get to within 0.00004% of the true value! Next his younger brother is going to code up a guess-and-check square root program, so that we can run the whole program using only elementary operations (+-*/). Then we'll try some optimizations on the pi algorithm, like the discrete integration I mentioned, and finish up with a sin/cos/tan calculator using only the pythagorean theorem. All of this makes a great introduction to programming, too. -J
Am I missing something? Why do you need square root at all? sqrt(x*x+y*y) < n is the same as x*x+y*y < n*n. Also, you don't need to check every point (an O(n^2) process). If you just follow around the edge of the circle (an O(n) process), you can easily infer the interior points. --ms Jason Holt wrote:
On Thu, 13 Jul 2006, Jason Holt wrote:
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
We coded up the algorithm and ran it for varying values of n. Here's the result for n=3,000,000, which took 23 hours on our Athlon64 3000+:
... area = 7068586468447 pi=3.141593985976 (pi/realpi=1.000000424112) Command exited with non-zero status 45
79435.25user 490.80system 23:26:48elapsed
Not the most quickly converging series, but it did get to within 0.00004% of the true value!
Next his younger brother is going to code up a guess-and-check square root program, so that we can run the whole program using only elementary operations (+-*/). Then we'll try some optimizations on the pi algorithm, like the discrete integration I mentioned, and finish up with a sin/cos/tan calculator using only the pythagorean theorem. All of this makes a great introduction to programming, too.
-J
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
On 7/18/06, Mike Speciner <speciner@ll.mit.edu> wrote:
Am I missing something? Why do you need square root at all? sqrt(x*x+y*y) < n is the same as x*x+y*y < n*n.
Also, you don't need to check every point (an O(n^2) process). If you just follow around the edge of the circle (an O(n) process), you can easily infer the interior points.
--ms
In fact, there's a clever algorithm for poltting the mandelbrot set that uses that fact: since points in the set take the longest to compute, we want to compute as few as possible. So divide up the space into quadrants and compute the edge values: if all the edge values are the same, then the interior is the same (assuming you don't have the entire set inside one quadrant). -- Mike Stay metaweta@gmail.com http://math.ucr.edu/~mike
On 7/18/06, Mike Stay <mike@math.ucr.edu> wrote:
In fact, there's a clever algorithm for poltting the mandelbrot set that uses that fact: since points in the set take the longest to compute, we want to compute as few as possible. So divide up the space into quadrants and compute the edge values: if all the edge values are the same, then the interior is the same (assuming you don't have the entire set inside one quadrant).
And if the edge values aren't all the same, recurse. -- Mike Stay metaweta@gmail.com http://math.ucr.edu/~mike
I just read a few of recent responces. An interesting thing is that nobody mentioned things that are really used for calculating pi, such as arctan or other hypergeometric series, or Simon Plouffe's formula etc. In addition to Richardson's method that I briefly described at the beginning of the thread, another not very often used approach is to use better than sin(x)<x<tan(x) inequalities. For example, sin(x)*(14+cos(x))/(9+6*cos(x)) < x < 8*sin(x)/9/(1+cos(x))+25*sin(x)/(36+9*cos(x)) Both inequalities can be certainly improved, and Richardson's method applied to them would give even better results, but even the cited inequalities give for x = pi/6 3.14156< pi < 3.14161 and for x=pi/12, 3.141592< pi < 3.141593 Compare that with Archimedes calculations; 96*sin(pi/96) < pi < 96*tan (pi/96) requires a lot of work, but gives only 3.141 < pi < 3.1427 Alec Mihailovs
On Wednesday 19 July 2006 03:43, Alec Mihailovs wrote:
In addition to Richardson's method that I briefly described at the beginning of the thread, another not very often used approach is to use better than sin(x)<x<tan(x) inequalities.
For example,
sin(x)*(14+cos(x))/(9+6*cos(x)) < x < 8*sin(x)/9/(1+cos(x))+25*sin(x)/(36+9*cos(x))
Are these less painful to find and prove than they look? :-) (I suppose you'd prove them by clearing denominators, expressing everything in terms of sin(kx) and cos(kx), and comparing terms in the Taylor series, or something like that. How do you find them?) -- g
From: "Gareth McCaughan" <gareth.mccaughan@pobox.com>
On Wednesday 19 July 2006 03:43, Alec Mihailovs wrote:
sin(x)*(14+cos(x))/(9+6*cos(x)) < x < 8*sin(x)/9/(1+cos(x))+25*sin(x)/(36+9*cos(x))
(I suppose you'd prove them by clearing denominators, expressing everything in terms of sin(kx) and cos(kx), and comparing terms in the Taylor series, or something like that. How do you find them?)
In Maple, I found (the beginning of) Taylor series for sin(x)*(a+b*cos(x))/(c+d*cos(x)), then solved the system of equations, the coefficient at x equals 1 and the coefficients at x^3 and x^5 equal 0. That gave series x - x^7/2100 + O(x^9) That proofs (or can be used to prove) the first inequality. Similarly for the second one, with the resulting series x + x^7/2800 + O(x^9) Alec
I said inf /===\ 2 | | 4 n %pi
(d93) | | ------------------- = --- | | (2 n - 1) (2 n + 1) 2 n = 1
which (Wallis) presumably got from Euler's infinite product for sin(x)/x.
Dick Askey chided that this formula preceded Euler by several decades, and pointed me at https://ceres.math.wisc.edu/horde/imp/attachment.php?u=askey&t=1153038348&f=... Sergey V (not N!) Khrushchev's reconstruction of Lord Brouncker's 1655 derivation of the first continued fraction for pi, which I would notate 4 2 1 -- + 1 = K(2, (2 k + 1) ) = 2 + -----------, pi 9 2 + ------- 25 2 + --- . . . whose reciprocal is termwise equivalent to Gregory's series, pi/4 = 1 - 1/3 + 1/5 - ... Using the matrix methods of gosper.org/stanfordn2.dvi, the corresponding CF for 1/atan(x) 2 2 2 2 K(x - 2 k (x - 1) + 1, (2 k + 1) x ) = --------------------------------------- - x. x The CF transformation in the next section, gosper.org/stanfordn3.dvi, makes this 2 2 K((2 k + 1) (1 - x) (x + 1), 4 (k + 1) x ) -------------------------------------------, x instead of the usual 2 2 K(2 k + 1, (k + 1) x ) -----------------------. x But note that x -> 1 in the previous gives Wallis's product, modulo regrouping. The transformation applied to the usual atan CF gives asec(x) 1 ------------ = ---------------------------------------. 2 2 2 sqrt(x - 1) (2 k + 1) (x - 1) 1 - x K(x + 2 k, -------------------) + ----- 4 2 A&S seems to give a CF for asin(x)/sqrt(1-x^2) but not just asin. The term-for-term formula, after a tweak, is 1 1 ------- = - asin(x) x x - ---------------------------------------------------------------, 1 2 2 3 3 2 2 K(2 (k + -) x + (k + 1) (2 k + 3), - 4 (k + 1) (k + -) x ) 2 2 which may have been too ugly for them. The transformation formula only applies to K(linear,quadratic), but I can probably find one that works. Also, Alec M. remarked of no discussion of the hypergeometric approach, which I believe most practical for enormous precision. My favorite is the Chudnovsky formula inf /===\ [ - (2 k - 1) (6 k - 5) (6 k - 1) 545140134 k - 531548725 ] | | [ ] | | [ 3 3 ] | | [ 0 72 53360 k ] k = 1 sqrt(10005) (homogeneously) = -------------- 25625606400 pi with which they extracted biwyuns of Saganic verses (MLB's term). Other series with higher "digits/term" have a quadratic surd in their term ratio, thus requiring > twice the computation/term. The AGM methods (incuding that pi_q method I mentioned) are asymptotically better, but the crossover point is astronomical. But I claim anyone computing digits instead of partial quotients is wasting computer time. --rwg
On Tue, 18 Jul 2006, Mike Speciner wrote:
Am I missing something? Why do you need square root at all? sqrt(x*x+y*y) < n is the same as x*x+y*y < n*n.
Yes, but don't tell the younger brother. :) Square roots were the limit of his skills, plus sticking with ordinary pythagoras made things easier on the older brother.
Also, you don't need to check every point (an O(n^2) process). If you just follow around the edge of the circle (an O(n) process), you can easily infer the interior points.
Indeed. I was ready to guide him toward my integration solution, where you just calculate the height of the quarter circle at each x, but he came up with a very simple alternative on his own. a. Start at (0,y+1) b. move down until (0,y) is interior (y--) c. move right (x++) d. goto b unless x>=r It was great to watch his expression when he saw how much faster it was than the n^2 solution. -J
participants (7)
-
Alec Mihailovs -
dasimov＠earthlink.net -
Gareth McCaughan -
Jason Holt -
Mike Speciner -
Mike Stay -
R. William Gosper