Bellard的基于BBP公式的Pi计算程序

bbp.pi.c

/*
 * Computation of the n'th decimal digit of \pi with very little memory.
 * Written by Fabrice Bellard on January 8, 1997.
 * 
 * We use a slightly modified version of the method described by Simon
 * Plouffe in "On the Computation of the n'th decimal digit of various
 * transcendental numbers" (November 1996). We have modified the algorithm
 * to get a running time of O(n^2) instead of O(n^3log(n)^3).
 * 
 * This program uses mostly integer arithmetic. It may be slow on some
 * hardwares where integer multiplications and divisons must be done
 * by software. We have supposed that 'int' has a size of 32 bits. If
 * your compiler supports 'long long' integers of 64 bits, you may use
 * the integer version of 'mul_mod' (see HAS_LONG_LONG).  
 */

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

/* uncomment the following line to use 'long long' integers */
/* #define HAS_LONG_LONG */

#ifdef HAS_LONG_LONG
#define mul_mod(a,b,m) (( (long long) (a) * (long long) (b) ) % (m))
#else
#define mul_mod(a,b,m) fmod( (double) a * (double) b, m)
#endif

/* return the inverse of x mod y */
int inv_mod(int x,int y) {
  int q,u,v,a,c,t;

  u=x;
  v=y;
  c=1;
  a=0;
  do {
    q=v/u;
    
    t=c;
    c=a-q*c;
    a=t;
    
    t=u;
    u=v-q*u;
    v=t;
  } while (u!=0);
  a=a%y;
  if (a<0) a=y+a;
  return a;
}

/* return (a^b) mod m */
int pow_mod(int a,int b,int m)
{
  int r,aa;
   
  r=1;
  aa=a;
  while (1) {
    if (b&1) r=mul_mod(r,aa,m);
    b=b>>1;
    if (b == 0) break;
    aa=mul_mod(aa,aa,m);
  }
  return r;
}
      
/* return true if n is prime */
int is_prime(int n)
{
   int r,i;
   if ((n % 2) == 0) return 0;

   r=(int)(sqrt(n));
   for(i=3;i<=r;i+=2) if ((n % i) == 0) return 0;
   return 1;
}

/* return the prime number immediatly after n */
int next_prime(int n)
{
   do {
      n++;
   } while (!is_prime(n));
   return n;
}

int main(int argc,char *argv[])
{
  int av,a,vmax,N,n,num,den,k,kq,kq2,t,v,s,i;
  double sum;

  if (argc<2 || (n=atoi(argv[1])) <= 0) {
    printf("This program computes the n'th decimal digit of \\pi\n"
     "usage: pi n , where n is the digit you want\n"
	   );
    exit(1);
  }
  
  N=(int)((n+20)*log(10)/log(2));

  sum=0;

  for(a=3;a<=(2*N);a=next_prime(a)) {

    vmax=(int)(log(2*N)/log(a));
    av=1;
    for(i=0;i<vmax;i++) av=av*a;

    s=0;
    num=1;
    den=1;
    v=0;
    kq=1;
    kq2=1;

    for(k=1;k<=N;k++) {

      t=k;
      if (kq >= a) {
	 do {
	    t=t/a;
	    v--;
	 } while ((t % a) == 0);
	 kq=0;
      }
      kq++;
      num=mul_mod(num,t,av);

      t=(2*k-1);
      if (kq2 >= a) {
	 if (kq2 == a) {
	    do {
	       t=t/a;
	       v++;
	    } while ((t % a) == 0);
	 }
	 kq2-=a;
      }
      den=mul_mod(den,t,av);
      kq2+=2;
      
      if (v > 0) {
	t=inv_mod(den,av);
	t=mul_mod(t,num,av);
	t=mul_mod(t,k,av);
	for(i=v;i<vmax;i++) t=mul_mod(t,a,av);
	s+=t;
	if (s>=av) s-=av;
      }
      
    }

    t=pow_mod(10,n-1,av);
    s=mul_mod(s,t,av);
    sum=fmod(sum+(double) s/ (double) av,1.0);
  }
  printf("Decimal digits of pi at position %d: %09d\n",n,(int)(sum*1e9));
  return 0;
}