/*
** Pascal Sebah : September 1999
**
** Subject:
**
** A very easy program to compute Pi with many digits.
** No optimisations, no tricks, just a basic program to learn how
** to compute in multiprecision.
**
** Formulae:
**
** Pi/4 = arctan(1/2)+arctan(1/3) (Hutton 1)
** Pi/4 = 2*arctan(1/3)+arctan(1/7) (Hutton 2)
** Pi/4 = 4*arctan(1/5)-arctan(1/239) (Machin)
** Pi/4 = 12*arctan(1/18)+8*arctan(1/57)-5*arctan(1/239) (Gauss)
**
** with arctan(x) = x - x^3/3 + x^5/5 - ...
**
** The Lehmer's measure is the sum of the inverse of the decimal
** logarithm of the pk in the arctan(1/pk). The more the measure
** is small, the more the formula is efficient.
** For example, with Machin's formula:
**
** E = 1/log10(5)+1/log10(239) = 1.852
**
** Data:
**
** A big real (or multiprecision real) is defined in base B as:
** X = x(0) + x(1)/B^1 + ... + x(n-1)/B^(n-1)
** where 0<=x(i)<B
**
** Results: (PentiumII, 450Mhz)
**
** Formula : Hutton 1 Hutton 2 Machin Gauss
** Lehmer's measure: 5.418 3.280 1.852 1.786
**
** 1000 decimals: 0.2s 0.1s 0.06s 0.06s
** 10000 decimals: 19.0s 11.4s 6.7s 6.4s
** 100000 decimals: 1891.0s 1144.0s 785.0s 622.0s
**
** With a little work it's possible to reduce those computation
** times by a factor 3 and more:
**
** => Work with double instead of long and the base B can
** be choosen as 10^8
** => During the iterations the numbers you add are smaller
** and smaller, take this in account in the +, *, /
** => In the division of y=x/d, you may precompute 1/d and
** avoid multiplications in the loop (only with doubles)
** => MaxDiv may be increased to more than 3000 with doubles
** => ...
*/
#include <time.h>
#include <stdio.h>
#include <malloc.h>
#include <math.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 = 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);
}
