/* 
   The Taylor series for the exponential function e^x:
    e^z = lim n|->inf \Sum_1^inf z^n / n!
   See: http://en.wikipedia.org/wiki/Power_series

   This version uses arbitrary precision arithmetic from the GMP library.


   Basic GMP API usage is taken 
   from: http://www.cs.colorado.edu/~srirams/classes/doku.php/gmp_usage_tutorial
   Retrieved: 25/1/2013 
   GMP API: http://gmplib.org/manual/

   Compile: gcc -O2 -o p_gmp power_gmp.c -lgmp -lm
   Run:     time ./p_gmp 4000 20
*/

#include "gmp.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <assert.h>

/* types */
typedef unsigned long int  ui;

/* prototypes */
static inline void fact(mpq_t res, ui n);
static inline void pow_int(mpq_t res, ui z, ui n) ;
static inline void power_e(mpq_t res, ui z, ui d);

/* variables */
double startTime, stopTime;

/* aux functions */
// res = n!
static inline
void fact(mpq_t res, ui n){
  ui i;
  mpz_t p;
 
  mpz_init_set_ui(p,1); /* p = 1 */
  for (i=1l; i <= n ; ++i){
    mpz_mul_ui(p,p,i); /* p = p * i */
  }
  mpq_set_z(res, p);
  mpz_clear(p);
  // return res;
}

// res = z^n
static inline
void pow_int(mpq_t res, ui z, ui n) {
  mpz_t tmp, zq;
  mpz_init(tmp);
  mpz_init(zq);
  mpz_set_ui(zq,z);
  mpz_pow_ui(tmp, zq, n);
  mpq_set_z(res, tmp);
  mpz_clear(tmp);
  mpz_clear(zq);
  // return res;
}

/* exponential function, result will be stored in res */
void power_e(mpq_t res, ui z, ui d) {
  mpq_t sum, old_sum, eps;
  mpq_t tmp_num, tmp_den, tmp, tmpe, epsilon;
  ui n; 
  mpq_init(epsilon);
  mpq_init(eps);
  mpq_init(old_sum);
  mpq_init(sum);
  mpq_init(tmp_num);
  mpq_init(tmp_den);
  mpq_init(tmp);
  mpq_init(tmpe);
  // mpq_set_ui(epsilon,1l,1000000l); // 1e-6
  pow_int(tmpe, 10l, d);              // 10^d
  // mpq_set_z(epsilon, tmpe);        
  mpq_inv(epsilon, tmpe);             // 10^-d
#if defined(DEBUG)
  fprintf(stderr, "\nPRECISION %d ie EPSILON: ",d);
  mpq_out_str(stderr, 10, epsilon);   
#endif
  mpq_set_ui(sum, 0l, 1l);
  mpq_set_ui(eps, 1l, 1l);
  for (n=0; mpq_cmp(eps, epsilon)>0; n++) {
#if defined(DEBUG)
    fprintf(stderr, "\nIT %d\t", n);
#endif
    mpq_set(old_sum, sum);
    pow_int(tmp_num, z, n);
    fact(tmp_den, n);
    mpq_div(tmp, tmp_num, tmp_den);
#if defined(DEBUG)
    fprintf(stderr, "\tBONZO \t");
    mpq_out_str(stderr, 10, tmp);   // mpq_out_str(stdout, 10, old_sum);   mpq_out_str(stdout, 10, eps);
#endif
    mpq_add(sum, sum, tmp);
    mpq_sub(eps, sum, old_sum);
#if defined(DEBUG)
    fprintf(stderr, "\tSUM \t");
    mpq_out_str(stderr, 10, sum);   // mpq_out_str(stdout, 10, old_sum);   mpq_out_str(stdout, 10, eps);
#endif
  }
  mpq_clear(tmp_num);
  mpq_clear(tmp_den);
  mpq_clear(tmp);
  mpq_clear(tmpe);
  mpq_clear(old_sum);
  mpq_set(res,sum);
  // return sum;
}

int main (int argc, char **argv) {
  ui n, d;
  mpq_t res;
  mpf_t resf;
  double zz;
  if (argc<2) {
    fprintf(stderr, "Usage: power <n> <d> ... computes e^n to precision 10^-d, using a Taylor series");
    exit(1);
  }
  if (argc==2) { // just for testing
    n = atol(argv[1]);
    d = 6; // default precision
  } else {
    n = atol(argv[1]);
    d = atol(argv[2]);
  }
  if (d==0) { // special case: basic testing
    mpq_t tmp;
    ui i, z;
    z = atol(argv[1]);
    n = atol(argv[2]);
    mpq_init(tmp);
    fprintf(stderr, "Testing basic functions for z^n and n! only ...\n");
    for(i=0; i<n; i++) {
      pow_int(tmp, z, i);
      printf("\n%d^%d = ", z, i);
      mpq_out_str(stdout, 10, tmp);
    }
    for(i=0; i<n; i++) {
      fact(tmp, i);
      printf("\n%d! = ", i);
      mpq_out_str(stdout, 10, tmp);
    }
    exit(1);
  }
  /* MAIN ------------------- */
  fprintf(stderr, "Computing e^%d ...\n", n);
  mpq_init(res);
  mpf_init(resf);

  startTime = clock();
  power_e(res, n, d);
  stopTime = clock();

  printf("\ne^%d = ", n);
  mpq_out_str(stdout, 10, res);
  mpf_set_q(resf,res);
  printf("\ne^%d = ", n);
  mpf_out_str(stdout, 10, d, resf);
  if (n<100) {
    zz = exp((double)n);
    printf("\nExpected result: %.9f\n", zz);
  }
  printf("\nPrecision: %d digits", d);
  printf("\nElapsed time: %f secs\n", 
	 (stopTime-startTime)/CLOCKS_PER_SEC);
  exit(0);
}