Liam's Website


Digits of PI

Made on 3-14-2015 Written in C

On PI day (3-14-15) I wrote this program to calculate 10 million digits of PI. It makes use of the Gauss–Legendre algorithm to calculate PI, and the GMP Library for arbitrary precision math.

Be Careful: tenmillion.pi.txt is 10 megabytes and will take a long time to load!!

/*
Copyright 2015 William Ryan

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
#include <gmp.h>

#define ITTERATIONS 50
#define PRECISTION  100000000
#define GMPPREC     800000000

mpf_t a;
mpf_t b;
mpf_t t;
mpf_t p;
mpf_t a_n;
mpf_t b_n;
mpf_t t_n;
mpf_t p_n;
mpf_t pi;
mpf_t tmp;
mpf_t tmp2;
mp_exp_t n_exp;

int counter = 0;

int main(int argc, char *argv[]){
	mpf_set_default_prec(GMPPREC);

	// Set b
	mpf_init_set_ui(b, 0);
	mpf_sqrt_ui(b, 2);
	mpf_ui_div(b, 1, b);


	mpf_init_set_ui(a, 1);
	//mpf_init_set_d(b, 1/sqrt(2));
	mpf_init_set_d(t, 0.25);
	mpf_init_set_ui(p, 1);
	mpf_init_set_ui(a_n, 0);
	mpf_init_set_ui(b_n, 0);
	mpf_init_set_ui(t_n, 0);
	mpf_init_set_ui(p_n, 0);
	mpf_init_set_ui(pi , 0);
	mpf_init_set_ui(tmp, 0);
	mpf_init_set_ui(tmp2, 0);

	#ifdef ONLYLAST
	fprintf(stderr, "Compiled to only print pi at the end.\n");
	#endif


	while (counter < ITTERATIONS){
		/* (a+b)/2; */
		mpf_add(tmp, a, b);
		mpf_div_ui(a_n, tmp, 2);

		mpf_set_ui(tmp, 0);

		/* sqrt(a*b) */
		mpf_mul(tmp, a, b);
		mpf_sqrt(b_n, tmp);

		mpf_set_ui(tmp, 0);

		/* t - (p * powl(a - a_n, 2))     :( */
		mpf_sub(tmp, a, a_n);
		mpf_pow_ui(tmp, tmp, 2);
		mpf_mul(tmp, tmp, p);
		mpf_sub(t_n, t, tmp);

		mpf_set_ui(tmp, 0);

		/*  2 * p */
		mpf_mul_ui(p_n, p, 2);

		mpf_set_ui(tmp, 0);

		/*  PI CALC: powl(a_n + b_n, 2) / (4 * t_n)*/
		mpf_add(tmp, a_n, b_n);
		mpf_pow_ui(tmp, tmp, 2);

		mpf_mul_ui(tmp2, t_n, 4);

		mpf_div(pi, tmp, tmp2);

		mpf_set_ui(tmp, 0);
		mpf_set_ui(tmp2, 0);

		// Now, for printing
		#ifndef ONLYLAST
		char *str_pi = mpf_get_str(NULL, &n_exp, 10, PRECISTION, pi);
		printf("%3u: pi=%s\n", counter, str_pi);
		free(str_pi);
		#endif
		#ifdef ONLYLAST
		fprintf(stderr, "i=%u\n",counter);
		#endif

		mpf_set(a, a_n);
		mpf_set(b, b_n);
		mpf_set(t, t_n);
		mpf_set(p, p_n);

		counter++; 
	}
	#ifdef ONLYLAST
	char *str_pi = mpf_get_str(NULL, &n_exp, 10, PRECISTION, pi);
	fprintf(stderr, "pi=");
	printf("%s\n", str_pi);
	free(str_pi);
	#endif
	return 0;
}
	
all:
	gcc pi.c -o pi -lm -lgmp -I/opt/local/include -L/opt/local/lib -DONLYLAST