Вычисление числа П - C++ - Ответ 2264588

#include <time.h> 
#include <stdio.h> 
#include <malloc.h> 
#include <math.h> 
#include <conio.h> 
long B=10000; /* Working base */ 
long LB=4;    /* Log10(base)  */ 
long MaxDiv=450;  /* about sqrt(2^31/B) */ 
/* 
** Set the big real x to the small integer Integer  
*/ 
void SetToInteger (long n, long *x, long Integer) { 
  long i; 
  for (i=1; i<n; i++) x[i] = 0; 
  x[0] = Integer; 
} 
/* 
** Is the big real x equal to zero ? 
*/ 
long IsZero (long n, long *x) { 
  long i; 
  for (i=0; i<n; i++)   
    if (x[i])   return 0; 
        return 1; 
} 
/* 
** Addition of big reals : x += y 
**  Like school addition with carry management 
*/ 
void Add (long n, long *x, long *y) { 
  long carry=0, i; 
  for (i=n-1; i>=0; i--) { 
    x[i] += y[i]+carry; 
    if (x[i]<B) carry = 0; 
    else { 
      carry = 1; 
      x[i] -= B; 
    } 
  }   
} 
/* 
** Substraction of big reals : x -= y 
**  Like school substraction with carry management 
**  x must be greater than y 
*/ 
void Sub (long n, long *x, long *y) { 
  long i; 
  for (i=n-1; i>=0; i--) { 
    x[i] -= y[i]; 
                if (x[i]<0) { 
                  if (i) {       
        x[i] += B; 
        x[i-1]--; 
      } 
                } 
  }   
} 
/* 
** Multiplication of the big real x by the integer q  
** x = x*q. 
**  Like school multiplication with carry management 
*/ 
void Mul (long n, long *x, long q) { 
  long carry=0, xi, i; 
  for (i=n-1; i>=0; i--) { 
    xi  = x[i]*q;                
    xi += carry;                 
    if (xi>=B) { 
      carry = xi/B; 
      xi -= (carry*B); 
    } 
    else  
      carry = 0; 
    x[i] = xi; 
        }   
} 
/* 
** Division of the big real x by the integer d  
** The result is y=x/d. 
**  Like school division with carry management 
**  d is limited to MaxDiv*MaxDiv. 
*/ 
void Div (long n, long *x, long d, long *y) { 
  long carry=0, xi, q, i; 
  for (i=0; i<n; i++) { 
    xi    = x[i]+carry*B; 
    q     = xi/d; 
    carry = xi-q*d;    
    y[i]  = q;         
  }   
} 
/* 
** Find the arc cotangent of the integer p (that is arctan (1/p)) 
**  Result in the big real x (size n) 
**  buf1 and buf2 are two buffers of size n 
*/ 
void arccot (long p, long n, long *x, long *buf1, long *buf2) { 
  long p2=p*p, k=3, sign=0; 
  long *uk=buf1, *vk=buf2; 
  SetToInteger (n, x, 0); 
  SetToInteger (n, uk, 1);      /* uk = 1/p */ 
  Div (n, uk, p, uk); 
  Add (n, x, uk);               /* x  = uk */ 
 
  while (!IsZero(n, uk)) { 
    if (p<MaxDiv) 
      Div (n, uk, p2, uk);  /* One step for small p */ 
    else { 
      Div (n, uk, p, uk);   /* Two steps for large p (see division) */ 
      Div (n, uk, p, uk);   
    } 
    /* uk = u(k-1)/(p^2) */ 
    Div (n, uk, k, vk);       /* vk = uk/k  */ 
    if (sign) Add (n, x, vk); /* x = x+vk   */ 
    else Sub (n, x, vk);      /* x = x-vk   */ 
    k+=2; 
    sign = 1-sign; 
  } 
} 
/* 
** Print the big real x 
*/ 
void Print (long n, long *x) { 
  long i;  
  printf ("%d.", x[0]); 
  for (i=1; i<n; i++) { 
    printf ("%.4d", x[i]); 
    if (i%25==0) printf ("%8d\n", i*4); 
  } 
  printf ("\n"); 
} 
/* 
** Computation of the constant Pi with arctan relations 
*/ 
void main () {   
  clock_t endclock, startclock;  
  long NbDigits=10000, NbArctan; 
  long p[10], m[10]; 
  long size=1+NbDigits/LB, i; 
  long *Pi      = (long *)malloc(size*sizeof(long)); 
  long *arctan  = (long *)malloc(size*sizeof(long)); 
  long *buffer1 = (long *)malloc(size*sizeof(long)); 
  long *buffer2 = (long *)malloc(size*sizeof(long));  
  startclock = clock();     
  /* 
  ** Formula used:  
  **    
  **   Pi/4 = 12*arctan(1/18)+8*arctan(1/57)-5*arctan(1/239) (Gauss) 
  */ 
  NbArctan = 3; 
  m[0] = 12; m[1] = 8;  m[2] = -5; 
  p[0] = 18; p[1] = 57; p[2] = 239;  
  SetToInteger (size, Pi, 0); 
  /* 
  ** Computation of Pi/4 = Sum(i) [m[i]*arctan(1/p[i])]  
  */ 
  for (i=0; i<NbArctan; i++) { 
    arccot (p[i], size, arctan, buffer1, buffer2); 
    Mul (size, arctan, abs(m[i])); 
    if (m[i]>0) Add (size, Pi, arctan);   
    else        Sub (size, Pi, arctan);   
  } 
  Mul (size, Pi, 4); 
  endclock = clock (); 
  Print (size, Pi);  /* Print out of Pi */ 
  printf ("Computation time is : %9.2f seconds\n",   
         (float)(endclock-startclock)/(float)CLOCKS_PER_SEC );  
  free (Pi); 
  free (arctan); 
        free (buffer1); 
        free (buffer2); 
        getch(); 
         
}