Newsgroups: comp.lang.c
Subject: followup to fooling with pi
From: pmor1@cs.aukuni.ac.nz (Peter William       Morrison       )
Date: Fri, 16 Oct 1992 00:43:43 GMT
Organization: Computer Science Dept. University of Auckland
Keywords: pi
Summary: Well, it's the summary, as requested!

Here's the followup to my original article about fooling with pi.
For those of you not in the know, this was a request posted by
myself (Peter Morrison) for the C wizards of the world to come up
with programs to calculate pi to 1,000,000 d.p.
There was a great response. Thanks all!

People to mention in my summary post to usenet for
"fooling with pi"

Marc_North@mindlink.bc.ca (Marc North)
Date: 30 Sep 1992 21:07 -0700 (PDT)
Wants a summary

William Kaufman <wkaufman@us.oracle.com>
Date: 30 Sep 1992 22:29:42 -0700 (PDT)
completely missed the point with "double pi = atan(1) * 4;"

auj@aberystwyth.ac.uk (Alun Jones)
Date: 01 Oct 1992 11:29:51 +0100
Sent a program in ml, one in perl, and one in C (finally!), mentioned other
bignum packages which would do it for me, and hinted at a few other sequences
for pi - also suggested which identity to use.

kajihara@pysgjb.physics.ncsu.edu (Scott Kajihara)
Date: 01 Oct 1992 09:05:53 -0400
Replied that the Taylor expansion method was VERY slow... yup!

feg@clyde.att.com (Forrest E Gehrke)
Date: 01 Oct 1992 09:05:19 -0400 (EDT)
Mentioned Microsoft's demo prog, provided with MS-C V4.0, that gave pi
to the requested accuracy up to 4 from the last d.p.

pjames@hubcap.clemson.edu (perry r james)
Date: 01 Oct 1992 10:10:25 -0400
Said to start with the smallest _x_ terms in the Taylor expansion, and
work backwards, since the errors would cancel out. Mentioned sci.math's FAQ
having routines to handle arbitrary precision arithmetic. Reminisced about
wanting to do the same but using completely rational arithmetic until the end

Ian.Collier@prg.oxford.ac.uk
Date: 01 Oct 1992 15:24:58 -0300 (BST)
Told me "You could always try (shock horror) writing it yourself..."
And gave me some other series:
	  [ 1    1^2     (1*3)^2   (1*3*5)^2   (1*3*5*7)^2       ]
 pi = 6 * [ - + ------ + ------- + --------- + ----------- + ... ]
	  [ 2   3!*2^3   5!*2^5     7!*2^7       9!*2^9          ]

 pi^2/6 = 1 + 1/4 + 1/9 + 1/25 + 1/36 + 1/49 + ...

	  [      1!    2!       3!          4!          ]
 pi = 2 * [ 1 + --- + ----- + ------- + --------- + ... ]
	  [     1*3   1*3*5   1*3*5*7   1*3*5*7*9       ]
and also "the most efficient calculation"...

a=1; b=sqrt(.5); c=.25; x=1; pi=0;
Repeat the following until the change in pi is sufficiently small:
   tempa=a
   c=c-x*(b-a)*(b-a)/4
   a=(a+b)/2
   b=sqrt(tempa*b)
   x=x+x
   pi=(a+b)*(a+b)/(4*c)
gave me the signature that implements this in 2 lines of C (cool!) to 800 d.p.

stephen@mont.cs.missouri.edu (Stephen Montgomery-Smith)
Date: 01 Oct 1992 11:41:33 -0500 (CDT)
Gave me a program to find e, with the comment that a simple modification should
produce pi.

stephen@mont.cs.missouri.edu (Stephen Montgomery-Smith)
Date: 01 Oct 1992 11:43:12 -0500 (CDT)
Also sent me a program by dik@cwi.nl (Dik T. Winter) that prints the
first 800 digits of pi.

stephen@mont.cs.missouri.edu (Stephen Montgomery-Smith)
Date: 01 Oct 1992 11:44:07 -0500 (CDT)
And again!!! What a guy!
Sent me another program to find 1000 digits of pi using the
'Spigot' algorithm, originally due to Stanly Rabinowitz

Jeff.Parker@East.Sun.COM (Jeff Parker - SunSoft)
Date: 01 Oct 1992 19:51:49 +0000 (GMT)
Did some very advanced algebra (ha ha) on my first identity - 
pi/4 = atn(1) = 1 - 1/3 + 1/5 - 1/7 + 1/9 ...
=> pi = 4(1 - 1/3 + 1/5 - 1/7 + 1/9 ...)
And mentioned that it converges REAL slow. Thanks...

dalessio@marble.rtsg.mot.com (Mario D'Alessio)
Date: 01 Oct 1992 15:58:45 -0500 (CDT)
Sent me the same program as Stephen Montgomery-Smith's 3rd email - 
that 'Spigot' algorithm again (giving due credit to Stanly Rabinowitz)
And a really neat signature...
************************************************************
* ******************************************************** *
* *  __________________           __                     * *
* *  \_________________|)____.---'--`---.____            * *
* *                ||    \----.________.----/            * *
* *  __/\__      __||____/ /_    `--'                    * *
* *  \RUSH/     |___         \   Mario D'Alessio         * *
* *   |/\|          `--------'   dalessio@rtsg.mot.com   * *
* *                                                      * *
* *    Work: (708) 632-2323  9am to 6pm Central Time     * *
* ******************************************************** *
************************************************************

CFKARLSO@apsc.com (CHARLES F KARLSON)
Date: 05 Oct 1992 16:51:40 -0700 (MST)
Sent me a another 'Spigot algorithm' program, but omitted any
credits, or explanation as to the workings...

Hmmm... there's a total lack of programs in this.
Here they all are:

#include <stdio.h>
unsigned short line[10000];
main(argc,argv)
int argc;
char *argv[];
{
   unsigned short n, a, i, q, v, dummy, h=0, w=0;

   if (argc==2)
       n=atoi(argv[1]);
   else
       n=20;

   n+=2;

   for (i=(n*10)/3 ; i > 0 ; i--) {
      line[i-1] = 2;
   }

   for (dummy=0; dummy < n; dummy++) {
       q=0;
       for (i=(n*10)/3 ; i > 0 ; i--) {
	   a = line[i-1] * 10 + (q*i);
	   q=a/(i*2-1);
	   line[i-1] = a%(i*2-1);
       }
       v = w;
       w = h+(q/10);
       h = q%10;
       if (w==10) {
	   w=0;
	   v++;
       }
       if (dummy!=0)
	   printf("%d", v);
   }
   printf("\n");
														     }

That one's not very accurate (i.e. not at all above 1000 d.p.)

This one is good:

/*
 * From: dik@cwi.nl (Dik T. Winter) Date: 30 Jun 91 12:39:11 GMT
 * Organization: CWI, Amsterdam -- dik t. winter, cwi, amsterdam, nederland
 * dik@cwi.nl
 *
 * The program below prints the first 800 digits of pi:
 */

int a = 10000, b, c = 28000, d, e, f[28001], g;

main()
{
    while ( b - c )
	f[b++] = a / 5;
    for (; d = 0, g = c * 2; c
      -= 14, printf("%.4d", e + d / a), e = d % a)
	for (b = c; d += f[b] * a, f[b] = d % --g,
	d /= g--, --b; d *= b);
}
/*  1000 digits of PI                                       */
/*  'Spigot' algorithm origionally due to Stanly Rabinowitz */

#include <stdio.h>
#include <math.h>

main()
{
  int d = 4, r = 10000, n = 251, m = 3.322*n*d;
  int i, j, k, q;
  static int a[3340];

  for (i = 0; i <= m; i++) a[i] = 2;
    a[m] = 4;

  for (i = 1; i <= n; i++) {
    q = 0;
    for (k = m; k > 0; k--) {
      a[k] = a[k]*r+q;
      q = a[k]/(2*k+1);
      a[k] -= (2*k+1)*q;
      q *= k;
    }
    a[0] = a[0]*r+q;
    q = a[0]/r;
    a[0] -= q*r;
    printf("%04d", q);
  }
  fflush(stdout);
}

This is the one I used to get pi to 50000 d.p. (It needs over 13MB to
do 1,000,000 d.p., and I've only got 6)

Ok, happy programming all you c buffs out there, and I'll be
back next year