/* 
   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!

   This is NOT a proper parallel solution to the above problem.
   It only sketches the communication for the specific input
   below, ie. input n=2 on 3 processors.

   Compile: mpicc -DDEBUG -DSLOW -lm -o power_par power_par.c
   Run:     mpirun -n 3 power_par 2 
*/

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

#define EPSILON 1e-6

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

#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 from, long to) {
  long n; 
  double sum, old_sum, eps;
  for (n=from, sum=0.0, eps=1.0; n<=to ; 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) {
  int p, id;
  long n;
  double z, zz;
  double elapsed_time;

  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &p);
  MPI_Comm_rank(MPI_COMM_WORLD, &id);

  if (id == 0) { /* master processor */
    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);

    // start the timer
    MPI_Barrier(MPI_COMM_WORLD);
    elapsed_time = - MPI_Wtime();

    /* SPECIAL CASE: INPUT n=2, and 2 workers */

    /* Q: how shall we split up the interval in general? */
    long from1=0, to1=7, from2=8, to2=14;
    double res, res1, res2, res_check;
    MPI_Status status;

    MPI_Send(&n, 1, MPI_LONG, 1, 0, MPI_COMM_WORLD);
    MPI_Send(&from1, 1, MPI_LONG, 1, 0, MPI_COMM_WORLD);
    MPI_Send(&to1, 1, MPI_LONG, 1, 0, MPI_COMM_WORLD);

    MPI_Send(&n, 1, MPI_LONG, 2, 0, MPI_COMM_WORLD);
    MPI_Send(&from2, 1, MPI_LONG, 2, 0, MPI_COMM_WORLD);
    MPI_Send(&to2, 1, MPI_LONG, 2, 0, MPI_COMM_WORLD);

    MPI_Recv(&res1, 1, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
    MPI_Recv(&res2, 1, MPI_DOUBLE, 2, 0, MPI_COMM_WORLD, &status);

    fprintf(stderr, "FROM PE1 res1=%1.9f\n", res1);
    fprintf(stderr, "FROM PE2 res2=%1.9f\n", res2);

    /* Q: how can we determine in general that this is close enough to the solution? */
    res = res1+res2;

    /* Q: what if this is not close enough? */

    // stop the timer
    elapsed_time += MPI_Wtime();

    res_check = exp((double)n);
    printf("e^%d = %1.9f (expected %1.9f)\nElapsed time: %f secs\n", 
	   n, res, res_check, elapsed_time);
  } else { // worker
    MPI_Status status;
    long from, to;
    double res;

    // start the timer
    MPI_Barrier(MPI_COMM_WORLD);
    // elapsed_time = - MPI_Wtime();

    MPI_Recv(&n, 1, MPI_LONG, 0, 0, MPI_COMM_WORLD, &status);
    MPI_Recv(&from, 1, MPI_LONG, 0, 0, MPI_COMM_WORLD, &status);
    MPI_Recv(&to, 1, MPI_LONG, 0, 0, MPI_COMM_WORLD, &status);

    fprintf(stderr, "[%d] from=%d, to=%d\t", id, from, to);
    res = power_e(n, from, to);
    fprintf(stderr, "res=%1.9f\n", res);

    MPI_Send(&res, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
  }

  MPI_Finalize();
  return 0;
}