Pi

pi.c

/*!
 *  \file   pi.c
 *  \brief  Calculates the nth digit of pi
 */

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

/*!
 * \brief Calculate prime numbers up to a given limit
 *
 * Returns an array of all primes up to a given value.  Works by using the
 * sieve of Erastothenes.  Start with the assumption that 2 is prime.  Cross
 * out all multiples of this value.  The next number not crossed out is a prime,
 * and so it is added to the list.  Multiples of this number are crossed out,
 * and the cycle repeats.
 *
 * This method maintains a static list of primes and a static sieve, along with
 * a high water mark of the \a to parameter.  If it is called again with \a to
 * less than the high water mark, cached results are returned.  If \a to is
 * higher than this mark, then all values between the high water mark and \a to
 * are sieved using the process above (with the modification that it must 
 * restart sieving at the beginning).  
 *
 * \param   to      All primes \f$ p \le to \f$ are returned
 * \param   count   Out parameter for the number of primes found
 *
 * \return  Array of prime integers
 */
static long* calc_primes(long to, size_t *count)
{
	long                   i, j;

	const int              primes_growth_rate = 2;      /* How fast to grow the primes array */
	const int              primes_start_size  = 100;    /* How many primes to allocate initially */
	
	static long            high_water_mark    = 0;      /* Maximum value of the \a to arg */
	static size_t          num_primes         = 0;      /* Number of primes found */
	static int             primes_alloced     = 0;      /* Amount of space in the primes array */
	static long           *primes             = NULL;   /* Array of prime numbers */
	static char           *sieve              = NULL;   /* Sieve; 1 => composite */
	
	/* First call */
	if(high_water_mark == 0) {
		/* Alloc the memory */
		sieve = malloc(to * sizeof(char));
		memset(sieve, 0, to * sizeof(char));
		primes = malloc(primes_start_size * sizeof(long));
		primes_alloced = primes_start_size;
		
		/* Prime the system (hehe) */
		sieve[0]        = 1;
		primes[0]       = 2;
		high_water_mark = 2;
		num_primes      = 1;
	}
	
	if(to <= high_water_mark) {
		/* Find the number of primes <= to */
		/* TODO: Better search algorithm */
		if(count) {
			*count = 0;
			
			for(i = num_primes - 1; i >= 0; i--) {
				if(primes[i] <= to) {
					*count = i + 1;
					break;
				}
			}
		}
		
		return primes;
	} else {
		/* Allocate some extra memory for the sieve */
		sieve = realloc(sieve, to * sizeof(char));
		memset(sieve + high_water_mark, 0, (to - high_water_mark) * sizeof(char));
		
		/* Divide all the values above the high water mark by all primes below
		   the high water mark */
		for(i = high_water_mark; i < to; i++) {
			for(j = 0; j < num_primes; j++) {
				if((i + 1) % primes[j] == 0) {
					sieve[i] = 1;
				}
			}
		}
		
		/* Now continue as a normal sieve--find the first 0, mark it prime, divide
		   the rest of the sieve by it, and continue until we hit to */
		for(i = high_water_mark; i < to; i++) {
			if(!sieve[i]) {
				if(primes_alloced <= num_primes) {
					primes_alloced *= primes_growth_rate;
					primes = realloc(primes, primes_alloced * sizeof(long));
				}
				primes[num_primes] = i + 1;
				num_primes++;

				for(j = i + 1; j < to; j++) {
					if(((j + 1) % (i + 1)) == 0) {
						sieve[j] = 1;
					}
				}
			}
		}
		
		high_water_mark = to;
		
		if(count) {
			*count = num_primes;
		}
		return primes;
	}
}

/*!
 *  \brief Factor of an integer
 *
 *  A single factor of an integer.  Stored as a linked list so that factors can
 *  easily be added as found and also be easily iterated
 */
struct factor {
	long   value;  /*!< Value of this factor         */
	struct factor  *next;   /*!< Next factor, or NULL if none */
};

/*!
 *  Frees the memory used by a list of factors
 *
 *  \param   factors   Factors to free
 */
static void free_factors(struct factor *factors)
{
	struct factor *prev = NULL;
	struct factor *cur  = NULL;
	
	cur = factors;
	while(cur) {
		prev = cur;
		cur = cur->next;
		free(prev);
	}
}

/*!
 *  Calculates all the prime factors of the given value.  Does this in a 
 *  simplistic manner, using trial division on all values below the sqrt(n).
 *
 *  \param   n   Value to factorize
 *
 *  \return  List of factors
 */
static struct factor *calc_prime_factors(long n, long max_factor, size_t *num_factors)
{
	long i;
	
	struct factor   *ret = NULL;
	struct factor   *tail = NULL;
	long            *primes;
	size_t           num_primes;
	
	primes = calc_primes(max_factor, &num_primes);
	*num_factors = 0;
	
	if(max_factor == 0) {
		return ret;
	}
	
	for(i = 0; i < num_primes; i++) {
		if(n % primes[i] == 0) {
			*num_factors = *num_factors + 1;
			if(!tail) {
				ret = malloc(sizeof(struct factor));
				tail = ret;
			} else {
				tail->next = malloc(sizeof(struct factor));
				tail = tail->next;
			}
			tail->value = primes[i];
			tail->next = NULL;
		}
	}
	
	if(!ret && n <= max_factor) {
		*num_factors = 1;
		ret = malloc(sizeof(struct factor));
		ret->value = n;
		ret->next = NULL;
	}
	
	return ret;
}

/*!
 *  Calculates \f$ b^e \bmod m \f$
 *
 *  \param   b   Base
 *  \param   e   Exponent
 *  \param   m   Modulus
 */
__attribute__((const))
static long expt_mod(long b, long e, long m)
{
	long ret = 1;
	
	while(e > 0) {
		if(e & 1) {
			ret = (ret * b) % m;
		}
		
		e >>= 1;
		b = (b * b) % m;
	}
	
	return ret;
}

/*!
 *  Calculates \f$ n \bmod m \f$.  Works properly for negative values.
 *
 *  \param   n   Value
 *  \param   m   Modulus
 */
__attribute__((const))
static inline long int_modulus(long n, long m)
{
	long res;
	res = n - m * (n / m);
	if(res < 0) {
		res += m;
	}
	
	return res;
}

/*!
 *  Calculates \f$ (a \cdot b) \bmod m \f$
 *
 *  \param   a   First multiplicand
 *  \param   b   Second multiplicand
 *  \param   m   Modulus
 */
__attribute__((const))
static inline long mod_mul(long a, long b, long m)
{
	return a * b - (long) ((1.0 / (double) m) * (double) a * (double) b) * m;
}

/*!
 *  Calculates the modular multiplicative of \a n with respect to modulo \a m.
 *  That is to say, it calculates a number \a x such that
 *  \f$ nx = 1 \pmod {m} \f$.
 *
 *  Uses the extended Euclidean algorithm to do the calculation.
 *
 *  \param   n   Number to invert
 *  \param   m   Modulus
 *
 *  \return  Multiplicative inverse of n mod m
 */
__attribute__((const))
static long mod_inv(long n, long m)
{
	long      x = 0;
	long prev_x = 1;
	long      y = 1;
	long prev_y = 0;
	
	long q = 0;
	long tmp;
	
	while(m) {
		q = n / m;
		
		tmp = m;
		m = int_modulus(n, m);
		n = tmp;
		
		tmp = x;
		x = prev_x - q * x;
		prev_x = tmp;
		
		tmp = y;
		y = prev_y - q * y;
		prev_y = tmp;
	}
	
	return prev_x;
}

/*!
 *  Calculates \f$ \sum_{j=0}^k{binom(n, j)} \bmod m \f$ where binom is the
 *  binomial coefficient of n and j.
 *
 *  \param   k   Summation upper limit
 *  \param   n   Upper term of binomial coefficient
 *  \param   m   Modulus
 */
__attribute__((const, hot))
static long mod_sum_binom(long k, long n, long m)
{
	long j;
	long A, B, C, C_acc;
	long a, b, a_star, b_star;
	
	size_t         num_factors_m;
	struct factor *prime_factors_m;
	struct factor *cur_fact;
	
	struct factor *r      = NULL;
	struct factor *cur_r  = NULL;
	
	if(k > n / 2) {
		return int_modulus(expt_mod(2, n, m) - mod_sum_binom(n - k - 1, n, m), m);
	}
	
	/* Step 1 */
	prime_factors_m = calc_prime_factors(m, k, &num_factors_m);
	
	/* Step 2 */
	A = 1; B = 1; C = 1;
	
	for(j = 0; j < num_factors_m; j++) {
		if(!cur_r) {
			r = malloc(sizeof(struct factor));
			cur_r = r;
		} else {
			cur_r->next = malloc(sizeof(struct factor));
			cur_r = cur_r->next;
		}
		cur_r->value = 1;
		cur_r->next = NULL;
	}
	
	for(j = 1; j <= k; j++) {
		
		/* Step 3a */
		a = n - j + 1;
		b = j;
		
		/* Steps 3b and 3c */
		cur_fact = prime_factors_m;
		cur_r = r;
		a_star = a;
		b_star = b;
		while(cur_fact) {
			while(a_star % cur_fact->value == 0) {
				a_star /= cur_fact->value;
				cur_r->value *= cur_fact->value;
			}
			
			while(b_star % cur_fact->value == 0) {
				b_star /= cur_fact->value;
				cur_r->value /= cur_fact->value;
			}
			
			cur_fact = cur_fact->next;
			cur_r = cur_r->next;
		}
		
		/* Step 3d */
		A = mod_mul(A, a_star, m);
		B = mod_mul(B, b_star, m);
	
		C = mod_mul(C, b_star, m);
	
		C_acc = A;
		for(cur_r = r; cur_r; cur_r = cur_r ->next) {
			C_acc = mod_mul(C_acc, cur_r->value, m);
		}
		
		C = int_modulus(C + C_acc, m);
	}
	
	free_factors(prime_factors_m);
	free_factors(r);
	
	/* Step 4 */
	return mod_mul(C, mod_inv(B, m), m);
}

/*!
 *  Gets \a n0 digits of \f$ \pi \f$ starting with digit \a n.  The first digit after the decimal point is
 *  digit 0.  Because the digits are returned in a double, the number of digits in the result cannot be higher
 *  than the precision of a double, regardless of n0.
 *
 *  This method only returns accurate results when \f$ n \ge 4 \cdot n0 \f$.
 *
 *  \param   n   Offset of digit to retrieve
 *  \param   n0  Number of digits to retrieve; also the precision of the result
 *
 *  \return  A decimal number whose integer part is 0 and whose decimal part corresponds to the nth digit of
 *           pi.  Another way to put it is, this is the decimal part of \f$ \pi \cdot 10^n \f$ , accurate to
 *           \a n0 decimal places.
 */
double get_pi_digits(long n, long n0)
{
	long k;
	long M, N, m, s;
	double b, c, x, t;
	int sign;
	double log_n;
	
	log_n = log(n);
	
	M = 2 * (long) (3 * n / log_n / log_n / log_n);
	N = (long) ((n + n0 + 1) * (log(10) / (log(2 * M_E * M)))) + 1;
	N += N % 2;
	
	b = 0;	
	for(k = 0; k < (M + 1) * N; k += 2) {
		x = 0;
		
		m = 2 * k + 1;
		s = expt_mod(10, n, m);
		s = mod_mul(4, s, m);
		x += (double) s / (double) m;
		
		m = 2 * k + 3;
		s = expt_mod(10, n, m);
		s = mod_mul(4, s, m);
		x -= (double) s / (double) m;
		
		b += x;
		if(b <= -0.5) {
			b += 1;
		} else if(b >= 1) {
			b -= 2;
		}
	}
	
	c = 0;
	sign = -1;
	for(k = 0; k < N; k++) {
		m = 2 * M * N + 2 * k + 1;
		s = mod_sum_binom(k, N, m);
		s = mod_mul(s, expt_mod(5, N, m), m);
		s = mod_mul(s, expt_mod(10, n - N, m), m);
		s = mod_mul(4, s, m);
		b += sign * (double) s / (double) m;
		b = b - floor(b);

		sign = -sign;
	}
	
	return modf(b, &t);
}

/*!
 *  Returns the nth decimal digit of \f$ \pi \f$.  Uses get_pi_digits() under the cover, but provides a nicer 
 *  API for getting a single digit.
 *
 *  \param   n  Offset of digit to retrieve
 *  \return  The given digit
 *  
 */
int get_pi_digit(long n) {
	const int initial_digits[] = { 
		1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 
		8, 9, 7, 9, 3, 2, 3, 8, 4, 6, 
		2, 6, 4, 3, 3, 8, 3, 2, 7, 9, 
		5, 0, 2, 8, 8, 4, 1, 9, 7, 1, 
		6, 9, 3, 9, 9, 3, 7, 5, 1, 0 };
	
	if(n < 50) {
		return initial_digits[n];
	} else {
		return (int) (10 * get_pi_digits(n, 14));
	}
}

static int print_digits(FILE *stream, long n, int num_digits) {
	int i;
	double t, val;
	
	val = get_pi_digits(n, num_digits);
	
	for(i = 0; i < num_digits; i++) {
		val *= 10;
		fprintf(stream, "%d", (int) val);
		val = modf(val, &t);
	}
}

int main(int argc, char **argv)
{
	long n;
	long N = 2000;
	double pct = 0;
	int last_pct = -1;
	FILE *out;
	
	out = fopen("pi.out", "w");

	fprintf(out, "3.\n");
	for(n = 0; n < 50; n++) {
		fprintf(out, "%d", get_pi_digit(n));
	}
	fprintf(out, "\n");
	
	for(n = 50; n < N; n += 10) {
		pct = (double) n / (double) N;
		if(floor(pct * 100) > last_pct) {
			fprintf(stderr, "\r%d%%", (int) (pct * 100));
			last_pct++;
		}
		print_digits(out, n, 10);
		if(n != 0 && (n + 10) % 50 == 0) {
			fprintf(out, "\n");
		}
	}
	fprintf(stderr, "\r100%%\n");
	fprintf(out, "\n");
	
	fclose(out);
}