91/03/03
To: Dr. Michael W. Ecker
909 Violet Terrace, Clarks Summit, PA 18411
From: Harry J. Smith, 19628 Via Monte Dr., Saratoga, CA 95070
Enclosure: MS-DOS disk (m2n2, revision 1.10).
Dear Mike,
Your "Of Pi" piece in EASY PIECES, January 1991 Algorithm, was
interesting. Especially to learn that the value of pi could be
estimated by finding, on average, over a very large collection of
integers from 1 to n, the number of ways an integer can be written
as the sum of two squares.
My program m2n2.Pas more or less verifies this, though it is not
a proof. An outline of the program is:
s_Max : Max s used for current estimate
n : n used in s = n*n + m*m
m : m used in s = n*n + m*m
n2 : n squared
m2 : m squared
s : Sum of two squares = n*n + m*m
t : t = total number of ways integers <= s_Max can be
written as the sum of two squares
a : Average number of ways, t / s_Max;
Error : Computed error = Pi - a
Done : Boolean done flag
s_Max <- 1;
repeat
n <- 0; t <- 0;
repeat
m <- -1; n <- n + 1; n2 <- n * n; Done <- True;
repeat
m <- m + 1; m2 <- m * m; s <- n2 + m2;
if (s <= s_Max) {
Done <- False;
if (m = 0) or (m = n)
then t <- t + 4
else t <- t + 8;
}
else
m <- n;
until (m = n);
until Done or (interrupted by operator);
a <- t / s_Max; Error <- Pi - a;
output a, Error, s_Max;
s_Max <- 2 * s_Max;
until (interrupted by operator);
To see how this program works, consider the triangular array of
s = m*m + n*n, n > 0, 0 <= m <= n, s <= s_Max = 128:
m \ n | 1 2 3 4 5 6 7 8 9 10 11 12
----- | ----------------------------------------------
0 | 1 4 9 16 25 36 49 64 81 100 121 ---
1 | 2 5 10 17 26 37 50 65 82 101 122
2 | 8 13 20 29 40 53 68 85 104 125
3 | 18 25 34 45 58 73 90 109 ---
4 | 32 41 52 65 80 97 116
5 | 50 61 74 89 106 125
6 | 72 85 100 117 ---
7 | 98 113 ---
8 | 128
Only the numbers in this array can be represented as the sum of
the squares of two integers. When m = n there are 4 ways for each
entry: (n, n), (n, -n), (-n, n), (-n, -n). When m = 0 there are 4
ways for each entry: (0, n), (0, -n), (n, 0), (-n, 0). When m > 0
and m <> n there are 8 ways for each entry: (m, n), (m, -n), (-m,
n), (-m, -n), (n, m), (n, -m), (-n, m), (-n, -m). In the case of
s_Max = 128 there are 8 entries with m = n, 11 entries with m = 0,
and 41 entries with m <> n. The estimate for pi then is (4*(8+11)
+ 8*41) / 128 = 3.15625.
While working on this program I noticed that for s_Max > 10000
or so, the error in the estimate of pi was less than an error
estimate of:
ErrEst = 0.5 / �(s_Max).
The output of the program looks like this:
Compute Pi by determining ways an integer = sum of 2 squares (m2n2)
Error estimate = ErrEst = 0.5 / SqRt(s_Max)
Pi~=4.00000000 Error=-0.85840735 ErrEst=0.50000000 s_Max=1
Pi~=4.00000000 Error=-0.85840735 ErrEst=0.35355339 s_Max=2
Pi~=3.00000000 Error= 0.14159265 ErrEst=0.25000000 s_Max=4
...
Pi~=3.14257813 Error=-0.00098547 ErrEst=0.00552427 s_Max=8192
...
Pi~=3.14158630 Error= 0.00000635 ErrEst=0.00034527 s_Max=2097152
...
Pi~=3.14159082 Error= 0.00000183 ErrEst=0.00003052 s_Max=268435456
Pi~=3.14159228 Error= 0.00000037 ErrEst=0.00002158 s_Max=536870912
Pi~=3.14159176 Error= 0.00000090 ErrEst=0.00001526 s_Max=1073741824
Pi~=3.14159--- Error= 0.00000--- ErrEst=0.00001079 s_Max=2147418112
The last line above has not actually been computed yet, but shows
the limit of the program. The error estimate given is a safe upper
bound of the absolute value of the error, but it would appear that
a much better error estimate could be developed for large values
of s_Max. The error estimate is important because we cannot prove
that the process converges to pi unless we have an upper bound of
the absolute value of the error that approaches zero as s_Max goes
to infinity.
The proof is now simple when we realize that we are approximating
the area of a circle by the number of unit squares inside of it.
A unit square, with its center at integer coordinates other than
(0, 0), being counted if its center is inside or on the circle.
Sincerely,
Harry
times since October 20, 2004.