code-appendix.txt

#include <math.h>
#include <gmp.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>

void leibniz();
void circle();
void wallis();
void ramanujan(unsigned long long times);
void ramanujanNW(int times);
void ramd(int times);
void factz();
void circleD(unsigned long long size);
int inCircle(unsigned long long a, unsigned long long b, unsigned long long size);
void ramanjuanRat(unsigned long long times);
void monteCarlo(unsigned long long times);
void chudnovsky(unsigned long long times);
void newton(unsigned long long times);
void doubFactz(unsigned long long x);

mpz_t z,f;
mpf_t times, n;
mpz_t radius;

int main(){
  //mpf_inp_str(
  unsigned long long clo = clock();
  
  //printf("%lu\n",ULONG_MAX);
  //for things using mpf
  /*mpf_init2(times,1024);
  mpf_init2(n,1024);
  mpf_set_ui(times,1);//enter terms/1000
  mpf_mul_ui(times,times,100);*/

  //mpf_set_ui(times,10);//enter terms/1000
  //mpf_mul_ui(times,times,1000);
  //mpf_init2(n,1024);

  //for things using mpz
  mpz_init(radius);
  mpz_set_ui(radius,128);//enter terms/1000
  mpz_mul_ui(radius,radius,1000);
  

  //ramanujan(6400);
  //newton(16000);
  chudnovsky(6400);
  //leibniz();
  //wallis();
  //circle();
  //unsigned long long numterms = 128000;
  //circle(numterms);
  //mpz_init(z); 
  //chudnovsky(1);
  //mpz_init(z);
  //mpz_set_ui(z,1);
  //doubFactz(5);
  //printf("%d\n", mpz_get_ui(z));
  //ramd(5);
  //ramanujan(5);
  //ramanjuanRat(5);
  //mpf_init2(z,1024);
  //mpf_set_ui(z,1);
  //factz(5);
  
  //leibniz(numterms);
  /*mpf_t x;
  mpf_init2(x,1024);
  mpf_set_d(x,M_PI);
  mpf_out_str(stdout, 10, 0, x);
  printf("\n");*/
  //wallis(100000);
  
  //mpf_t x;
  //mpf_init2(x, 1024);
  //mpf_set_ui(x,100000);
  //mpf_add_ui(x, x, 3);
  //mpf_out_str(stdout,10,0,z);
  //printf("\n");
  //circle(4800);
  //circleD(100000);
  //monteCarlo(100000);
  //newton(1000);
  printf("%ld\n",(clock()-clo)/1000);
  //printf("Terms:%d\n", numterms);
  //sprintf("%e\n", M_PI);
  
  
  return 0;
}
void leibniz(){
  //  printf("runtohere\n");
  mpf_t pi, term,numer;
  mpf_init2(pi,1024);
  mpf_init2(term,1024);
  mpf_init2(numer,1024);
  mpf_set_ui(n,0);

  //int a;
  //for (a=0; a<times; a++){
  while(mpf_cmp(n,times)){
    //mpf_set_si(numer,(a%2)?(-1):(1));
    if (mpf_cmp_ui(numer,0)<=0)
      mpf_set_ui(numer,1);
    else
      mpf_set_si(numer,-1);
    //printf("run\n");
    //mpf_set_ui(term,(2*a+1));
    mpf_set(term,n);
    mpf_mul_ui(term,term,2);
    mpf_add_ui(term,term,1);
    mpf_div(term,numer,term);
    mpf_add(pi,pi,term);
    /*if (!(a%(times/100))) {
      mpf_t pi4;
      mpf_init2(pi4, 1024);
      mpf_mul_ui(pi4,pi,4);
      mpf_out_str(stdout, 10, 0, pi4);
      printf("\n");
      }*/
    //if (!(a%(times/100)))
    //printf("%d\n",a/(times/100));
    mpf_add_ui(n,n,1);
    }
  //return pi*4;
  mpf_mul_ui(pi, pi, 4);
  mpf_out_str(stdout, 10, 0, pi);
  printf("\n");
}
void circle(){
  mpz_t in, temp,temp2,size2,x,y;
  //in=in circle
  //temp=x^2
  //temp2=y^2
  //size2=size^2
  mpz_init2(in,1024);
  mpz_init2(temp,1024);
  mpz_init2(temp2,1024);
  mpz_init2(size2,1024);
  mpz_init2(x,1024);
  mpz_init2(y,1024);
  mpz_set(size2,radius);
  mpz_mul(size2,size2,radius);


  
  //unsigned long long x, y;
  //int f=0;//checking that the denom is right
  //x=0;
  //y=0;
  mpz_set_ui(x,0);
  mpz_set_ui(y,0);
  while (1){
    //printf("loop run\n");
    //f++;
    mpz_set(temp,x);
    mpz_mul(temp,temp,x);
    mpz_set(temp2,y);
    mpz_mul(temp2,temp2,y);
    mpz_add(temp,temp,temp2);
    
    if (mpz_cmp(temp,size2)<=0){
      mpz_add_ui(in,in,1);
      //x++;
      mpz_add_ui(x,x,1);
    }
    else{
      mpz_add_ui(y,y,1);
      //f+=(size-x);
      //y++;
      //if (size<y)
      if (mpz_cmp(y,radius)>0)
	break;
      //x=0;
      mpz_set_ui(x,0);
    }
  }
  mpf_t pi,tot;
  mpf_init2(pi,1024);
  mpf_set_z(pi,in);
  mpf_init2(tot,1024);
  mpf_set_z(tot,radius);
  mpf_add_ui(tot,tot,1);
  mpf_mul(tot,tot,tot);
  //mpz_out_str(stdout, 10, in);
  //printf("\n");
  //mpf_out_str(stdout, 10, 0, tot);
  //printf("\n");
  mpf_div(pi,pi,tot);
  mpf_mul_ui(pi,pi,4);
  //printf("Size:%d\n",size);
  mpf_out_str(stdout, 10, 0, pi);
  printf("\n");
  //printf("%d,%d\n",f,(size+1)*(size+1));

}
void wallis(){
  mpf_t pi, numer, term,temp;
  //mpf_div_ui(times,times,2);
  mpf_init2(pi,1024);
  mpf_init2(numer,1024);
  mpf_init2(term, 1024);
  mpf_init2(temp,1024);
  mpf_set_ui(pi,1);
  mpf_set_ui(n,1);
  //int a;
  //for (a=1; a<=times/2; a++){
  while(mpf_cmp(times,n)){
    //printf("run\n");
    mpf_set_ui(numer,2);
    mpf_mul(numer,numer,n);
    mpf_mul(numer,numer,numer);
    mpf_set(term,n);
    mpf_mul_ui(term,term,2);
    mpf_sub_ui(term,term,1);
    //mpf_set(pi,

    //mpf_mul_ui(term,term,2*a+1);
    mpf_set_ui(temp,2);
    mpf_mul(temp,temp,n);
    mpf_add_ui(temp,temp,1);
    mpf_mul(term,term,temp);
    mpf_div(term,numer,term);
    mpf_mul(pi,pi,term);
    //if (!(a%(times/200)))
    //printf("%d\n",a*2/(times/100));
    mpf_add_ui(n,n,1);
    }
  mpf_mul_ui(pi, pi, 2);
  mpf_out_str(stdout, 10, 0, pi);
  printf("\n");
}
void ramanujan(unsigned long long times){
  long long prec = 256*1024;
  mpf_t pi,term,tmp;
  mpf_init(pi);
  mpf_init(term);
  //mpf_init(mult);
  mpf_init(tmp);
  mpz_init(z);
  mpz_set_ui(z,1);
  //mpf_init2(rpi,16384);
  mpf_set_prec(pi,prec);
  mpf_set_prec(term,prec);
  //mpf_set_prec(mult,16384);
  mpf_set_prec(tmp,prec);
  //mpf_set_ui(mult,1);
  long long n;
  mpf_set_ui(pi,0);
  //mpf_out_str(stdout, 10, 0, tmp);
  //printf("\n");
  for (n=0; n<times; n++){
    mpf_set_ui(term,1103);
    mpf_set_ui(tmp,26390);
    mpf_mul_ui(tmp,tmp,n);
    mpf_add(term,term,tmp);
    mpz_set_ui(z,1);
    factz(4*n);
    mpf_set_z(tmp,z);
    mpf_mul(term,term,tmp);
    mpf_set_ui(tmp,396);
    mpf_pow_ui(tmp,tmp,4*n);
    mpf_div(term,term,tmp);
    mpz_set_ui(z,1);
    factz(n);
    mpf_set_z(tmp,z);
    mpf_pow_ui(tmp,tmp,4);
    mpf_div(term,term,tmp);
    mpf_add(pi,pi,term);
    //mpf_out_str(stdout, 10, 0, term);
    //printf("\n");
  }//while (mpf_cmp_d(rpi,M_PI)>0);
  mpf_sqrt_ui(tmp,2);
  mpf_mul_ui(tmp,tmp,2);
  mpf_div_ui(tmp,tmp,9801);
  //mpf_set_d(tmp,(2*sqrt(2))/9801);
  mpf_mul(pi,pi,tmp);
  mpf_ui_div(pi,1,pi);
  mpf_out_str(stdout, 10, 0, pi);
  printf("\n");
  //mpf_set_d(tmp,M_PI);
  //mpf_sub(tmp,pi,tmp);
  //printf("delta=\n");
  //mpf_out_str(stdout, 10, 0, tmp);
  //printf("\n");


  mpz_set_ui(z,1);
}
void factz(unsigned long long x){
  if (x==1||!x)
    return;
  mpz_mul_ui(z,z,x);
  factz(x-1);
}
void chudnovsky(unsigned long long times){
  long long prec = 524288;
  mpf_t tmp,pi,term;
  mpz_init(z);
  mpf_init2(tmp,prec);
  //mpf_init2(tmp2,1024);
  mpf_init2(pi,prec);
  mpf_init2(term,prec);
  mpf_set_ui(pi,0);
  long long a;
  for (a=0; a<times; a++){
    mpf_set_ui(term,13591409);
    mpf_set_ui(tmp,545140134);
    mpf_mul_ui(tmp,tmp,a);
    mpf_add(term,term,tmp);
    mpz_set_ui(z,1);
    factz(6*a);
    mpf_set_z(tmp,z);
    mpf_mul(term,term,tmp);
    //if (a%2)
      //mpf_neg(term,term);
    mpz_set_ui(z,1);
    factz(3*a);
    mpf_set_z(tmp,z);
    mpf_div(term,term,tmp);
    mpz_set_ui(z,1);
    factz(a);
    mpf_set_z(tmp,z);
    mpf_pow_ui(tmp,tmp,3);
    mpf_div(term,term,tmp);
    mpf_set_si(tmp,-640320);
    mpf_pow_ui(tmp,tmp,3*a);
    mpf_div(term,term,tmp);
    mpf_add(pi,pi,term);
    //mpf_out_str(stdout,10,0,term);
    //printf("\n");
  }
  mpf_sqrt_ui(tmp,10005);
  mpf_mul_ui(tmp,tmp,426880);
  mpf_div(pi,tmp,pi);
  mpf_out_str(stdout, 10, 0, pi);
  printf("\n");
}
void doubFactz(unsigned long long x){
  if (x==1||!x)
    return;
  mpz_mul_ui(z,z,x);
  doubFactz(x-2);
}
void newton(unsigned long long times){
  mpf_t pi, term, tmp;
  mpf_init2(pi,16384);
  mpf_init2(term,16384);
  mpf_init2(tmp,16384);
  mpz_init(z);
  mpf_set_ui(pi,0);
  long long a;
  for (a=0; a<times; a++){
    //printf("started\n");
    mpz_set_ui(z,1);
    factz(a);
    mpf_set_z(term,z);
    mpz_set_ui(z,1);
    doubFactz(2*a+1);
    mpf_set_z(tmp,z);
    mpf_div(term,term,tmp);
    mpf_add(pi,pi,term);
    //printf("ended\n");
    //mpf_out_str(stdout,10,0,term);
    //printf("\n");
  }
  mpf_mul_ui(pi,pi,2);
  mpf_out_str(stdout, 10, 0, pi);
  printf("\n");
}