bellard's algorithm for pi calculation (slow version)

pi.c

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

#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 the prime number immediatly after n */
int next_prime(int n)
{
    do
	{
		n++;
    } while (!is_prime(n));
    return n;
}

void process(const char* input)
{
	int av, a, vmax, N, n, num, den, k, kq, kq2, t, v, s, i;
    double sum;
    
    n = atoi(input);
	if(n <= 0)
	{
		fprintf(stderr, "This program computes the n'th decimal digit of \\pi\nusage: pi n , where n is the digit you want\nthe digit to calculate must be greather than 0 but is [%s]\n", n);
		return;
	}
	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));
}

/* time macros */
#ifdef _WIN32
#	include <windows.h> 
#	pragma comment (lib, "WINMM.LIB")
#	define get_current_time() ( timeGetTime() )
#else
#	define get_current_time() ( time(NULL) )
#endif

#define render_duration(s,e) ( fprintf(stdout, "Process duration : %lld ms\n", (long long) (e-s)) )

/**
 * main program
 */
int main(int argc, char** argv)
{
	
	time_t start_time, end_time;
	
    if (argc < 2) {
		fprintf(stderr, "This program computes the n'th decimal digit of \\pi\n"
		       "usage: pi n , where n is the digit you want\n");
		exit(1);
    }

	start_time = get_current_time();
    process(argv[1]);
    end_time = get_current_time();
    
    render_duration(start_time, end_time);
	
	return 0;
}