eeg@ski.UUCP (eeg systems ) (07/23/85)
I am posting this as an interesting approach to solving
pi to many decimal places. Originally compiled and run on an Apple
II, it took 40 minutes to solve pi to 1000 places. On a Vax, it took
1.5 seconds! I limited its precision to 1000 places because the only
listing for pi I had available only listed that many places.
To use, cut below and compile.
Enjoy. Bryan Costales ( ...!dual!ptsfa!ski!eeg!bcx )
--------------------------< cut here >-------------------------------
/*
** NAME
** pi -- a program to calculate the value of pi
**
** SYNOPSIS
** pi {-}digits {t}
**
** DESCRIPTION
** Pi calculates the value of pi to the number of digits
** specified. Output is to the standard output, if the
** digits specifed is negative, output is formatted
** for attractive display on an 80col printer.
**
** The optional last argument "t" enables timing of the
** computational loop. The pitime() function must be supplied
** by the individual site implimentation.
**
** This is an adaptation in C of the BASIC program "Apple Pi"
** by Robert J. Bishop, published in Micro Apple (c) 1981
** Micro Ink, Inc. That program, when run on an Apple II,
** took 40 hours to solve pi to 1000 places. Arrays of digits
** are used to impliment extended precision arithmatic. The
** program solves pi using the series expansion:
**
** infinity infinity
** ____ 16*(-1e(k+1)) ____ 4*(-1e(k+1))
** \ \
** pi = > ------------- - > ------------
** / /
** ---- (2k-1)*5e(2k-1) ---- (2k-1)*239e(2k-1)
** k = 1 k = 1
**
** AUTHOR
** Bryan Costales ( ...!dual!ptsfa!ski!eeg!bcx )
** pi.c (c) 1985 Bryan Costales
**
** BUGS
**
*/
#include <stdio.h>
/*
* The maximum number of digits (plus 5 for rounding).
*/
#define MAX 10005
/*
* Global variables
*/
int Sign,
Zero,
Size,
Pass, /* two passes */
Exp, /* exponent for divide */
Divide;
/*
* These character arrays represent the digits for extended
* precision arithmatic. (These may need to be int on some machines).
*/
#define DIGIT char
DIGIT Power[MAX], Term[MAX], Result[MAX];
main(argc, argv)
int argc;
char *argv[];
{
static int constant[3]={0,25,239};
register int i;
int charcount, printer = 0;
void exit(); /* for lint */
if(argc == 1)
{
(void)fprintf(stderr,"usage: pi {-}digits {t}\n");
exit(1);
}
Size = atoi(argv[1]);
if(Size < 0)
{
Size = -Size;
printer++;
}
if(Size >= MAX-4 || Size <= 0)
{
(void)fprintf(stderr,"'digits' must be 1 thru %d.\n",MAX-5);
exit(1);
}
if( (argc==3) && (*argv[2]=='t') )
{
(void)fprintf(stderr,"start computation: ");
pitime();
}
for(Pass=1; Pass < 3; Pass++)
{
init();
do
{
copy();
Divide = Exp;
div(Term);
if ( Sign > 0 )
add();
if ( Sign < 0 )
sub();
Exp = Exp + 2;
Sign *= -1;
Divide = constant[Pass];
div(Power);
if (Pass == 2)
div(Power);
} while( Zero != 0 );
}
if( (argc == 3) && (*argv[2] == 't') )
{
(void)fprintf(stderr,"end computation: ");
pitime();
}
if( printer )
{
(void)printf("\t\tThe value of pi" );
(void)printf(" to %d decimal places\n\n\t\t",Size);
}
(void)printf("%d.",Result[0]);
charcount = 0;
for(i = 1; i <= Size; i++)
{
(void)printf( "%d", Result[i]);
charcount++;
if ( (charcount == 50) && printer )
{
(void)printf( "\n\t\t " );
charcount = 0;
}
}
(void)printf( "\n" );
return( 0 ); /* for lint */
}
add()
{
register DIGIT *r, *t;
register int sum, carry = 0;
r = Result + Size;
t = Term + Size;
while( r >= Result )
{
sum = *r + *t + carry;
carry = 0;
if( sum >= 10 )
{
sum -= 10;
carry++;
}
*r-- = sum;
--t;
}
}
sub()
{
register DIGIT *r, *t;
register int diff, loan = 0;
r = Result + Size;
t = Term + Size;
while( r >= Result )
{
diff = *r - *t - loan;
loan =0;
if( diff < 0 )
{
diff += 10;
loan++;
}
*r-- = diff;
--t;
}
}
div(array)
register DIGIT *array;
{
register DIGIT *end;
register int quotient, residue, digit = 0;
Zero = 0;
for( end = array + Size; array <= end; )
{
digit += *array;
quotient = digit / Divide;
residue = digit % Divide;
if((Zero !=0) || ((quotient+residue) != 0))
Zero = 1;
else
Zero = 0;
*array++ = quotient;
digit = 10 * residue;
}
}
init()
{
register DIGIT *p, *t, *r, *end;
p = Power;
t = Term;
r = Result;
end = Power+Size;
while( p <= end )
{
*p++ = 0;
*t++ = 0;
if (Pass == 1)
*r++ = 0;
}
*Power = 16 / (Pass * Pass);
if( Pass == 1 )
{
Divide = 5;
}
else
{
Divide = 239;
}
div(Power);
Exp = 1;
Sign = 3 - (2 * Pass);
}
copy()
{
register DIGIT *t, *p, *end;
t = Term;
p = Power;
end = Term + Size;
while (t <= end)
{
*t++ = *p++;
}
}
pitime()
{
/* Print the current time \n */
/* Here we just beep */
(void)fprintf( stderr, "%c\n", 0x07 );
}
cagordon@watnot.UUCP (Chris A. Gordon) (07/25/85)
In article <187@ski.UUCP> eeg@ski.UUCP (eeg systems (bcx) writes: [ article deleted - program to calculate pi to X digits ] >** infinity infinity >** ____ 16*(-1e(k+1)) ____ 4*(-1e(k+1)) >** \ \ >** pi = > ------------- - > ------------ >** / / >** ---- (2k-1)*5e(2k-1) ---- (2k-1)*239e(2k-1) >** k = 1 k = 1 >** Here is a more simple sum-evaluation of pi (thought I don't know if it will work with the original program): oo ---- k-1 | | \ (-1) pi = +--+ > ---------- | / 2k-1 ---- k=1
ken@boring.UUCP (07/27/85)
I decided to rewrite the posted program in "bc". (If you think having
only 286 distinct names in BASIC is bad, try using "bc"!) To get more
decimal places I had to use scaling, see program for details.
Fortunately, we know 1 < pi < 10 :-).
It took 2 hours to compute pi to a 1000 places on a 780 but it is
limited only by the amount of memory "dc" can get. And at last I have
something to soak up all those spare CPU cycles.
Anybody know a faster converging series for pi?
Enjoy, Ken
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# pi.b
# pi.sh
# This archive created: Fri Jul 26 15:37:57 1985
# By: Ken Yap ()
export PATH; PATH=/bin:$PATH
if test -f 'pi.b'
then
echo shar: will not over-write existing file "'pi.b'"
else
cat << \SHAR_EOF > 'pi.b'
/*
** NAME
** pi -- a program to calculate the value of pi
**
** SYNOPSIS
** bc pi.b
**
** DESCRIPTION
** Pi calculates the value of pi to the number of digits
** specified.
**
** This is an adaptation in BC of the BASIC program "Apple Pi"
** by Robert J. Bishop, published in Micro Apple (c) 1981
** Micro Ink, Inc. That program, when run on an Apple II,
** took 40 hours to solve pi to 1000 places. The
** program solves pi using the series expansion:
**
** infinity infinity
** ____ 16*(-1e(k+1)) ____ 4*(-1e(k+1))
** \ \
** pi = > ------------- - > ------------
** / /
** ---- (2k-1)*5e(2k-1) ---- (2k-1)*239e(2k-1)
** k = 1 k = 1
**
** AUTHORS
** Bryan Costales ( ...!dual!ptsfa!ski!eeg!bcx )
** pi.c (c) 1985 Bryan Costales
**
** Modified for bc by Ken Yap (..!seismo!{rochester,mcvax}!ken)
** pi.b
**
** BUGS
**
*/
scale = 99;
d/*en*/ = 1;
define i(p/*laces*/) {
auto x;
for (x = 0; x < p/*laces*/; x++) d/*en*/ *= 10;
"Result will be pi times 10 to the power of ";
return (p/*laces*/);
}
define p/*i*/(k) {
k *= 2;
s/*ign*/ = 1;
x = 0;
y = 0;
m = 5;
n = 239;
for(e/*xp*/ = 1; e/*xp*/ < k; e/*xp*/ += 2) {
f = d/*en*/ / (e/*xp*/ * m);
g = d/*en*/ / (e/*xp*/ * n);
if ( s/*ign*/ > 0 ) {
x += f;
y += g;
}
if ( s/*ign*/ < 0 ) {
x -= f;
y -= g;
}
m *= 25;
n *= 57121;
s/*ign*/ *= -1;
}
return(x * 16 - y * 4);
}
SHAR_EOF
if test 1529 -ne "`wc -c < 'pi.b'`"
then
echo shar: error transmitting "'pi.b'" '(should have been 1529 characters)'
fi
fi # end of overwriting check
if test -f 'pi.sh'
then
echo shar: will not over-write existing file "'pi.sh'"
else
cat << \SHAR_EOF > 'pi.sh'
#! /bin/sh
bc -c pi.b > pi.d << PROG
i(1000)
p(1000)
PROG
nice dc < pi.d > pi.1000 &
SHAR_EOF
if test 85 -ne "`wc -c < 'pi.sh'`"
then
echo shar: error transmitting "'pi.sh'" '(should have been 85 characters)'
fi
chmod +x 'pi.sh'
fi # end of overwriting check
# End of shell archive
exit 0
--
UUCP: ..!{seismo,okstate,garfield,decvax,philabs}!mcvax!ken Voice: Ken!
Mail: Centrum voor Wiskunde en Informatica, Kruislaan 413, 1098 SJ, Amsterdam.mjs@eagle.UUCP (M.J.Shannon) (07/28/85)
> In article <187@ski.UUCP> eeg@ski.UUCP (eeg systems (bcx) writes: > [ article deleted - program to calculate pi to X digits ] > >** infinity infinity > >** ____ 16*(-1e(k+1)) ____ 4*(-1e(k+1)) > >** \ \ > >** pi = > ------------- - > ------------ > >** / / > >** ---- (2k-1)*5e(2k-1) ---- (2k-1)*239e(2k-1) > >** k = 1 k = 1 > >** > > Here is a more simple sum-evaluation of pi (thought I don't know if it will work > with the original program): > > oo > ---- k-1 > \ (-1) > pi = 4 > ---------- > / 2k-1 > ---- > k=1 The reason why the 2 infinite sums are used is that the single infinite sum converges *VERY* slowly. -- Marty Shannon UUCP: ihnp4!eagle!mjs Phone: +1 201 522 6063
ctk@ecsvax.UUCP (Tim Kelley) (07/29/85)
In article <11307@watnot.UUCP> cagordon@watnot.UUCP (Chris A. Gordon) writes: >In article <187@ski.UUCP> eeg@ski.UUCP (eeg systems (bcx) writes: >[ article deleted - program to calculate pi to X digits ] >>** infinity infinity >>** ____ 16*(-1e(k+1)) ____ 4*(-1e(k+1)) >>** \ \ >>** pi = > ------------- - > ------------ (Expression 1) >>** / / >>** ---- (2k-1)*5e(2k-1) ---- (2k-1)*239e(2k-1) >>** k = 1 k = 1 >>** > >Here is a more simple sum-evaluation of pi (thought I don't know if it will work >with the original program): > > oo > ---- k-1 > | | \ (-1) > pi = +--+ > ---------- (Expression 2) > | / 2k-1 > ---- > k=1 Friends, I hate to include the complete text of the article I'm responding to but I have no alternative. The reason expression 1 is better than expression 2 is that the series in 1 converges faster than that in 2. One can check this out with a programmable calculator. Expression 2 will run all night and give you maybe 3 figures. If you want 3 digits of accuracy from expression 2 you'd need about 1000 terms; 4 digits would require 10,000. -- C.T. Kelley decvax!mcnc!ecsvax!ctk Dept. of Math. N.C. State U. Box 8205 Raleigh, N.C. 27695-8205, 919-737-7895
kl@unido.UUCP (07/29/85)
# This is another program to calculate a few decimals of pi, using the
# "standard" series for computing. There are two versions given:
# pi, which uses a normal C-program (div.c) to do multi-precision division,
# and pix, which uses an assembler routine (exdiv.s) instead. Exdiv.s contains
# the special "extended division" opcode for a VAX, thereby allowing twice
# the number of digits per integer (see def.h).
# To use:
# 1. Unpack this file
# 2. Edit def.h
# It should only be necessary to define the number of DIGITS
# to compute, or the default LOG file to use.
# 3. If you don't have rusage-stuff (4.2BSD), edit status() in pi.c
# 4. do make
# 5. run pi or pix, an argument is the logfile to use.
#
# Note: a status() is given at the beginning and end of execution.
# Since pi (pix) will probably be a long running program :-), a
# HUP (#1) signal sent to pi (pix) will cause a status entry in
# the LOG file. For really long running progs, see ./clock.
#
# Some sample (cpu) times (on a VAX/750):
# pix pi
# DIGITS
# 1000 about 12 sec. about 25 sec.
# 10000 about 20 min. about 44 min.
# 100050 about 35 hrs. about 78 hrs.
#
# This is a shell archive. Unpack with "sh <file".
echo x - Makefile
cat > "Makefile" << '!Funky!Stuff!'
CFLAGS = -O
all: pi pix print
pi: def.h pi.o div.o print
cc $(CFLAGS) pi.o div.o -o pi -lm
pix: def.h pix.o exdiv.o print
cc $(CFLAGS) pix.o exdiv.o -o pix -lm
pix.o: pi.c def.h
cc -DEXDIV -S pi.c
as -o pix.o pi.s
rm pi.s
exdiv.o: exdiv.s def.h
/lib/cpp exdiv.s tmp.s -DEXDIV
as -o exdiv.o tmp.s
rm tmp.s
print: print.o
cc $(CFLAGS) -o print print.o
clean:
-rm *.o core pi pix print
!Funky!Stuff!
echo x - clock
cat > "clock" << '!Funky!Stuff!'
if [ x$1 = x ]
then
echo "Usage: $0 <pid of pi-program>"
exit
fi
while true
do
kill -1 $1
sleep 3600
done
!Funky!Stuff!
echo x - def.h
cat > "def.h" << '!Funky!Stuff!'
#define DIGITS 10000 /* no. of digits of pi we wish to compute */
#define LOG "log" /* the logfile */
#ifdef EXDIV
# define SIZE 100000000 /* max. power of ten <= maxint */
# define NSIZE 8 /* no. of digits in SIZE, ie 10^NSIZE = SIZE */
#else
# define SIZE 10000 /* max. power of ten <= sqrt(maxint) */
# define NSIZE 4 /* no. of digits in SIZE, ie 10^NSIZE = SIZE */
#endif EXDIV
#define N (DIGITS/NSIZE+2) /* +2 for rounding */
!Funky!Stuff!
echo x - div.c
cat > "div.c" << '!Funky!Stuff!'
#include "def.h"
div (a,n)
int *a, n;
{
register int *ra, *ul, q;
ul = &a[N];
q = 0;
for (ra=a; ra<ul;) {
q = q*SIZE + *ra;
*ra++ = q/n;
q = q % n;
}
}
!Funky!Stuff!
echo x - exdiv.s
cat > "exdiv.s" << '!Funky!Stuff!'
#include "def.h"
LL0:
.data
.comm _quad,8
.text
.align 1
.globl _div
_div:
.word L12
jbr L14
L15:
clrq _quad
clrl r11
L18:
cmpl r11,$N
jgeq L17
emul $SIZE,_quad,*4(ap)[r11],_quad
ediv 8(ap),_quad,*4(ap)[r11],_quad
clrl _quad+4
L16:
incl r11
jbr L18
L17:
ret
.set L12,0xc00
L14:
jbr L15
.data
!Funky!Stuff!
echo x - pi.c
cat > "pi.c" << '!Funky!Stuff!'
#include "def.h"
add (a,b)
int *a, *b;
{
register int *ap, *bp;
for (ap = &a[N], bp = &b[N]; ap>a; ) {
*--bp += *--ap;
if (*bp > SIZE-1) {
*bp -= SIZE;
bp[-1]++;
}
}
}
sub (a,b)
int *a, *b;
{
register int *ap, *bp;
for (ap = &a[N], bp = &b[N]; ap>a; ) {
*--bp -= *--ap;
if (*bp < 0) {
*bp += SIZE;
bp[-1]--;
}
}
}
copy(x,y)
int *x, *y;
{
register int *rx, *ry;
for (rx = &x[N], ry = &y[N]; rx>x; )
*--ry = *--rx;
}
#include <sys/time.h>
#include <sys/resource.h>
#include <stdio.h>
#define PRINT(x) fprintf(Log,"x: %d\t", ru.x)
char *ctime();
struct rusage ru;
FILE *Log;
int cnt5, cnt239;
int lim5, lim239;
status()
{
long now;
time(&now);
fprintf (Log,"at %s", ctime(&now));
fprintf(Log,"cnt5= %d(%d)\tcnt239= %d(%d)\n",4*cnt5,lim5,
4*cnt239,lim239);
getrusage(RUSAGE_SELF,&ru);
fprintf(Log,"usr time: %d.%06d [sec]\n", ru.ru_utime.tv_sec,
ru.ru_utime.tv_usec);
fprintf(Log,"sys time: %d.%06d [sec]\n", ru.ru_stime.tv_sec,
ru.ru_stime.tv_usec);
/* edit these PRINT's according to taste */
PRINT(ru_maxrss); PRINT(ru_minflt); PRINT(ru_majflt);
fprintf(Log,"\n");
PRINT(ru_nswap); PRINT(ru_inblock); PRINT(ru_oublock);
fprintf(Log,"\n");
PRINT(ru_nsignals); PRINT(ru_nvcsw); PRINT(ru_nivcsw);
fprintf(Log,"\n\n");
fflush(Log);
};
#include <math.h>
#include <signal.h>
FILE *popen();
main(argc, argv)
int argc;
char *argv[];
{
int term[N], pterm[N], sum[N];
register int i, j;
char aux[20];
FILE *Print;
cnt5 = cnt239 = 0;
signal(SIGHUP,status);
lim5 = (int)((double)DIGITS/log10(5.0)) + 10; /* +10 for safety */
lim239 = (int)((double)DIGITS/log10(239.0)) + 10;
if (argv[1] != NULL)
Log = fopen(argv[1],"w");
else
Log = fopen(LOG,"w");
if (Log == NULL) {
perror("fopen");
exit(1);
}
fprintf (Log,"DIGITS = %d\nlim5 = %d\nlim239 = %d\n",
DIGITS,lim5,lim239);
fprintf (Log,"start of exec. ");
status();
term[0] = 32*(SIZE/100);
for (j=1;j<N;j++)
term[j] = 0;
copy (term,sum);
i = 1;
while (i < lim5) {
i += 2;
div(term,25);
copy (term,pterm);
div(pterm,i);
sub(pterm,sum);
i += 2;
div(term,25);
copy(term,pterm);
div(pterm,i);
add(pterm,sum);
cnt5++;
}
term[0] = 4*(SIZE/10);
for (j=1;j<N;j++) term[j] = 0;
div(term,239);
sub(term,sum);
i = 1;
while (i < lim239) {
i += 2;
div(term,57121);
copy (term,pterm);
div(pterm,i);
add(pterm,sum);
i += 2;
div(term,57121);
copy(term,pterm);
div(pterm,i);
sub(pterm,sum);
cnt239++;
}
Print = popen("print","w");
for (j=0;j<N-2;j++) {
sprintf(aux,"%0*d",NSIZE,sum[j]);
fprintf(Print,"%s",aux);
}
putchar('\n');
fclose(Print);
fprintf (Log,"end of exec. ");
status();
fclose(Log);
}
!Funky!Stuff!
echo x - print.c
cat > "print.c" << '!Funky!Stuff!'
#include <stdio.h>
main() {
char c;
int i;
i = 0;
putchar('\n');
c = getchar();
putchar (c);
putchar(',');
putchar('\n');
while ( (c=getchar()) != EOF) {
putchar(c);
i++;
if (i%5 == 0 ) putchar(' ');
if (i%50 == 0) putchar('\n');
}
putchar('\n');
}
!Funky!Stuff!
#
#
# Karlheinz Lehner kl@unido.UUCP
# Univ. of Dortmund hpfcla!hpbbn!unido!kl
# West Germany mcvax!unido!kl
#grayson@uiucuxc.Uiuc.ARPA (08/02/85)
Here is the fastest way ever, and it is not a series. It is called
the Gauss-Legendre method, and was discovered independently by Salamin and
Brent. It is based on the arithmetic-geometric mean, which is:
agm(a,b):
loop until a - b is very small
amean = (a+b)/2; gmean = sqrt(a*b);
a = amean; b = gmean;
return a;
pi
a = 1; b = 1/sqrt(2); t = 1/4; x = 1;
loop until a - b is very small
y = a; a = (a+b)/2; b = sqrt(b*y); t = t - x*(a-y)^2; x = 2*x;
return pi = (a+b)^2 / (4*t);
With this method they've gotten 10000000 digits in 4 hours.
Sources:
R. P. Brent, Multiple-precision zero-finding methods and the complexity
of elementary funciton evaluation, in "Analytic Computation
Complexity", ed. J. F. Traub, Academic Press 1976
R. P. Brent, ACM Transactions on mathematical Software, Algorithm 524,
4 (1978) 57-70, 71-81.
R. P. Brent, MP source code, in fortran on the CYBER.
from:
uucp: {ihnp4,pur-ee}!uiucdcs!uiucuxc!grayson
Dan Grayson, Math Dept, Univ of Ill, Urbana 61801jmc@inset.UUCP (John Collins) (08/09/85)
In article <11307@watnot.UUCP> cagordon@watnot.UUCP (Chris A. Gordon) writes: >Here is a more simple sum-evaluation of pi (thought I don't know if it will >work with the original program): > > oo > ---- k-1 > | | \ (-1) > pi = +--+ > ---------- > | / 2k-1 > ---- > k=1 The convergence of this series is *TERRIBLE* You'd have to evaluate many thousands of terms to get anywhere (you are still messing around with the third decimal place after 500 terms!). Moreover the alternation of the signs makes life even more difficult - each term cancels most of the work of the previous term. Use of Aitken acceleration has a dramatic improvement though - but the resulting series can be shown to be equivalent to the series for 6 * arcsin(.5) which was given in a previous message. -- John M Collins ....mcvax!ist!inset!jmc Phone: +44 727 57267 Snail Mail: 47 Cedarwood Drive, St Albans, Herts, AL4 0DN, England.
gabriel@anl-mcs (John Gabriel) (08/11/85)
I am entering this conversation not having seen previous correspondence, but here is a trick from my high school days (in England during WWII). arctan(x+h) = arctan(x) + arctan(h/(1 + h*x + x**2)) Therefore arctan(1) = arctan(1/2) + arctan(1/3) = arctan(1/4) + arctan(2/9) + arctan(1/6) + arctan(3/19) and so on etc. If we stop at the above 4 way split, the error after N terms is roughly 1 bit in 2*N. Thus for a 64 bit floating point number with a 56 bit fractional part we need only 28 terms. A further split reduces the problem to 8 series with the worst being a power series in 1/8. Algebraic ingeneuity can combine terms so that your may not be doing quite as much work as summing all the series. After 8 terms you lose as much as you gain from further splits. Better solutions for use in subroutine libraries exist and are given in "A Software Manual for the Elementary Functions" by W.J. Cody and W.M. Waite, Prentice Hall. However if you want to compute Pi to a few thousand decimal places other considerations arise like Do you want to use rational fraction arithmetic?, should you be using methods like continued fraction expansions?, and on and on. I seem to recollect results sort of like pi/10 = arctan(r) where r was some reasonable rational fraction, but I can't for the life of me remember how it was obtained. There are of course the results from the fact that a polygon with 2**n + 1 sides can be solved using nothing worse than square roots, and so once you have that solution you also have a rapidly convergent series, e.g. for arctan(2*pi/17). But you still must compute the square roots. On the whole methods using repeated fractions are probably best. For details about that, get in touch with Henry Thacher c/o Dept. of Computer Science, U. of Kentucky, Lexington KY. by snail mail or telephone. Henry has retired and lives somewhere in CA now, but the Dept should still know where to find him. John Gabriel
bjchrist@convex.UUCP (08/19/85)
As an aside: If you know PI you can check a random number generator by the accuracy with which it calculates pi by the following means: |a square with coordinates: (0,0) (1,0) (1,1) (0,1) has an area |of one (1) | given: |a circle with radius 1 has as its area PI. Plot the circle on graph |paper with is center at (0,0). The area of the circle falling |in quadrant 1 is PI/4. So the ratio of points within the square's |area AND within the circle's area is PI/4. write a program to generate two random numbers 0<= x <= 1.0 keep count of the number of times sqrt(x**2 + y**2) <=1.0 the ratio of the hits (sqrt(x^2 + y^2) <= 1.0) to misses will give you PI/4. The rate of convergence is horrible but it's interesting to compare random number generators and their accuracy. If a REAL *random* number generator is used. After an infinite number of shots, you will really have PI. Assuming you have an infinitly accurate calculator.