pi.c

计算PI(圆周率)的程序

/* +++Date last modified: 05-Jul-1997 */

/*
** This method implements the Salamin / Brent / Gauss Arithmetic-
** Geometric Mean pi formula.
**
** Let A[0] = 1, B[0] = 1/Sqrt(2)
**
** Then iterate from 1 to 'n'.
** A[n] = (A[n-1] + B[n-1])/2
** B[n] = Sqrt(A[n-1]*B[n-1])
** C[n]^2 = A[n]^2 - B[n]^2  (or) C[n] = (A[n-1]-B[n-1])/2
**                        n
** PI[n] = 4A[n+1]^2 / (1-(Sum (2^(j+1))*C[j]^2))
**                       j = 1
**
** There is an actual error calculation, but it comes out  to  slightly
** more  than double on each iteration.  I think it results in about 17
** million  correct  digits,  instead  of  16  million  if  it actually
** doubled.  PI16 generates 178,000 digits.  PI19 to  over  a  million.
** PI22 is 10 million, and PI26 to 200 million.
**
** For what little it's worth, this program is placed into the public
** domain by its author, Carey Bloodworth, on September 21, 1996.
**
** The math routines in this program are not general purpose routines.
** They have been optimized and written for this specific use.  The
** concepts are of course portable, but not the implementation.
**
** The program run time is about O(n^1.7).  That's fairly good growth,
** compared to things such as the classic arctangents.  But not as good
** as it should be, if I used a FFT based multiplication.  Also, the 'O'
** is fairly large.  In fact, I'd guess that this program could compute
** one million digits of pi in about the same time as my previously
** posted arctangent method, in spite of this one having less than n^2
** growth.
**
** The program has not been cleaned up.  It's written rather crudely
** and dirty.  The RSqrt() in particular is rather dirty, having gone
** through numerous changes and attempts at speeding it up.
** But I'm not planning on doing any more with it.  The growth isn't as
** low as I'd hoped, and until I find a faster multiplication, the
** method isn't any better than simpler arctangents.
**
** I currently use a base of 10,000 simply because it made debugging
** easier.  A base of 65,536 and modifying the FastMul() to handle sizes
** that aren't powers of two would allow a bit better efficiency.
*/

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

#ifdef __TURBOC__
typedef short int Indexer;
#else
typedef long int Indexer;
#endif

typedef short int Short;
typedef long  int Long;

Short *a, *b, *c, *Work1, *MulWork, *Work2, *Work3;
int SIZE;

int ispow2(int x){      return! ((~(~0U>>1)|x)&x -1) ;}
/*
** Work1 is explicitly used in Reciprocal() and RSqrt()
** Work1 is implicitly used in Divide() and Sqrt()
** Work2 is explicitly used in Divide() and Sqrt()
** Work3 is used only in the AGM and holds the previous reciprocal
**  square root, allowing us to save some time in RSqrt()
*/


void Copy(Short *a, Short *b, Indexer Len)
{
      while (Len--) a[Len] = b[Len];
}

/*
** This rounds and copies a 2x mul result into a normal result
** Our number format will never have more than one unit of integer,
** and after a mul, we have two, so we need to fix that.
*/

void Round2x(Short *a, Short *b, Indexer Len)
{
      Indexer x;
      Short carry;

      carry = 0;
      if (b[Len+1] >= 5000)
            carry = 1;
      for (x = Len; x > 0; x--)
      {
            carry += b[x];
            a[x-1] = carry % 10000;
            carry /= 10000;
      }
}

void DivBy2(Short *n, Indexer Len)
{
      Indexer x;
      Long temp;

      temp = 0;
      for (x = 0; x < Len; x++)
      {
            temp = (Long)n[x]+temp*10000;
            n[x] = (Short)(temp/2);
            temp = temp%2;
      }
}

void DoCarries(Short *limit, Short *cur, Short carry)
{
      Long temp;

      while ((cur >= limit) && (carry != 0))
      {
            temp = *cur+carry;
            carry = 0;
            if (temp >= 10000)
            {
                  carry = 1;
                  temp -= 10000;
            }
            *cur = temp;
            --cur;
      }
}

void DoBorrows(Short *limit, Short *cur, Short borrow)
{
      Long temp;
      while ((cur >= limit) && (borrow != 0))
      {
            temp = *cur-borrow;
            borrow = 0;
            if (temp < 0)
		{
			borrow = 1;
			temp += 10000;
		}
            *cur = temp;
            --cur;
      };
}

void PrintShort2(char *str, Short *num, Indexer Len)
{
      Indexer x;

      printf("%s ", str);
      printf("%u.", num[0]);
      for (x = 1; x < Len; x++)
            printf("%04u", num[x]);
      printf("\n");
}

void PrintShort(char *str, Short *num, Indexer Len)
{
      Indexer x;
      int printed = 0;

      printf("%s ", str);

      printf("%u.\n", num[0]);

      for (x = 1; x < Len; x++)
      {
            printf("%02d", num[x]/100);
            printed += 2;
            if ((printed % 1000) == 0)
            {
                  printf("\n\n\n\n");
                  printed = 0;
            }
            else if ((printed % 50) == 0)
            {
                  printf("\n");
            }
            else if ((printed % 10) == 0)
            {
                  printf(" ");
            }
            printf("%02d", num[x] % 100);
            printed += 2;
            if ((printed % 1000) == 0)
            {
                  printf("\n\n\n\n");
                  printed = 0;
            }
            else if ((printed % 50) == 0)
            {
                  printf("\n");
            }
            else if ((printed % 10) == 0)
            {
                  printf(" ");
            }
      }
      printf("\n");

}

/* sum = a + b */

Short Add(Short *sum, Short *a, Short *b, Indexer Len)
{
      Long s, carry;
      Indexer x;

      carry = 0;
      for (x = Len - 1; x >= 0; x--)
      {
            s = (Long)a[x] + (Long)b[x] + carry;
            sum[x] = (Short)(s%10000);
            carry = s/10000;
      }
      return carry;
}

/* dif = a-b */

Short Sub(Short *dif, Short *a, Short *b, Indexer Len)
{
      Long d, borrow;
      Indexer x;

      borrow = 0;
      for (x = Len - 1; x >= 0; x--)
      {
            d = (Long)a[x] - (Long)b[x] - borrow;
            borrow = 0;
            if (d < 0)
            {
                  borrow = 1;
                  d += 10000;
            }
            dif[x] = (Short)d;
      }
      return borrow;
}

void Negate(Short *num, Indexer Len)
{
      Indexer x;Long d, borrow;

      borrow = 0;
      for (x = Len - 1; x >= 0; x--)
      {
            d = 0 - num[x] - borrow;
            borrow = 0;
            if (d < 0)
            {
                  borrow = 1;
                  d += 10000;
            }
            num[x] = (Short)d;
      }
}

/* prod = a*b.  prod should be twice normal length */

void Mul(Short *prod, Short *a, Short *b, Indexer Len)
{
      Long p;
      Indexer ia, ib, ip;

      if ((prod == a) || (b == prod))
      {
            printf("MUL product can't be one of the other arguments\n");
            exit(1);
      }

      for (ip = 0; ip < Len * 2; ip++)
            prod[ip] = 0;
      for (ib = Len - 1; ib >= 0; ib--)
      {
            if (b[ib] == 0)
                  continue;
            ip = ib + Len;
            p = 0;
            for (ia = Len - 1; ia >= 0; ia--)
            {
                  p = (Long)a[ia]*(Long)b[ib] + p + prod[ip];
                  prod[ip] = p%10000;
                  p = p/10000;
                  ip--;
            }
            while ((p) && (ip >= 0))
            {
                  p += prod[ip];
                  prod[ip] = p%10000;
                  p = p/10000;
                  ip--;
            }
      }
}

/*
** This is based on the simple O(n^1.585) method, although my
** growth seems to be closer to O(n^1.7)
**
** It's fairly simple.  a*b is: a2b2(B^2+B)+(a2-a1)(b1-b2)B+a1b1(B+1)
**
** For a = 4711 and b = 6397, a2 = 47 a1 = 11 b2 = 63 b1 = 97 Base = 100
**
** If we did that the normal way, we'd do
** a2b2 = 47*63 = 2961
** a2b1 = 47*97 = 4559
** a1b2 = 11*63 =  693
** a1b1 = 11*97 = 1067
**
** 29 61
**    45 59
**     6 93
**       10 67
** -----------
** 30 13 62 67
**
** Or, we'd need N*N multiplications.
**
** With the 'fractal' method, we compute:
** a2b2 = 47*63 = 2961
** (a2-b1)(b1-b2) = (47-11)(97-63) = 36*34 = 1224
** a1b1 = 11*97 = 1067
**
** 29 61
**    29 61
**    12 24
**    10 67
**       10 67