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

   Compile: gcc -DSLOW -lm -o power power.c
            (add -DDEBUG for debugging info)
   Run:     ./power 2
            (to compute e^2)
*/

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

#define EPSILON 1e-6

/* prototypes */
double power_e(long z);       // compute e^z
long pow_int(long z, long n); // compute z^n
long fact(long n);            // compute n!

/* variables */
double startTime, stopTime;

#if defined(SLOW)

# warning "SLOW operations enabled"

/* a very slow algorithm, to make it more compute-intensive */

long add_slow(long m, long n) {
  long sum, i;
  for (sum=m, i=1; i<=n; sum+=1, i++) { }
  return sum;
}

long mult_slow(long m, long n) {
  long sum, i;
  for (sum=0, i=1; i<=n; sum=add_slow(sum,m), i++) { }
  return sum;
}

/* a really silly implementation of square, to make it more compute-intensive */
long sqr(long n) { return mult_slow(n,n); }

/* factorial: n! */
long fact(long n) {
  long i, prod; 
  for (i=1, prod=1; i<=n; prod=mult_slow(prod,i), i++) { };
  return prod;
}

/* power over integers: z^n */
long pow_int(long z, long n) {
  long i, prod; 
  for (i=1, prod=1; i<=n; prod=mult_slow(prod,z), i++) { };
  return prod;
}

#else
/* the proper algorithm */

/* square computation */
long sqr(long n) { return n*n; }

/* factorial: n! */
long fact(long n) {
  long i, prod; 
  for (i=1, prod=1; i<=n; prod=prod*i, i++) { };
  return prod;
}

/* power over integers: z^n */
long pow_int(long z, long n) {
  if (n==0)
    return 1;
  if (n==1)
    return z;
  if (n%2==0) 
    return sqr(pow_int(z,n/2));
  else 
    return z*sqr(pow_int(z,n/2));
}
#endif

/* exponential function */
double power_e(long z) {
  long n; 
  double sum, old_sum, eps;
  for (n=0, sum=0.0, eps=1.0; eps > EPSILON; n++) {
#if defined(DEBUG)
    fprintf(stderr, "IT %d\t", n);
#endif
    old_sum = sum;
    sum += ((double)pow_int(z,n)) / ((double) fact(n));
    eps = sum - old_sum;
#if defined(DEBUG)
    fprintf(stderr, "SUM %1.9f OLD_SUM %1.9f EPS %1.9f\n", sum, old_sum, eps);
#endif
  }
  return sum;
}

int main (int argc, char **argv) {
  long n;
  double z, zz;
  if (argc<2) {
    fprintf(stderr, "Usage: power <n> ... computes e^n, using a Taylor series");
    exit(1);
  }
  n = atol(argv[1]);
  fprintf(stderr, "Computing e^%d ...\n", n);
  startTime = clock();
  z = power_e(n);
  stopTime = clock();
  zz = exp((double)n);
  printf("e^%d = %1.9f (expected %1.9f)\nElapsed time: %f secs\n", 
	 n, z, zz, (stopTime-startTime)/CLOCKS_PER_SEC);
  exit(0);
}