/* integrate.c
 *
 * A program that performs numeric integration of polynomial functions
 * using various methods.  Note that we can integrate polynomials
 * analytically, so this program only serves to compare the calculated
 * integratal (using various methods) with the actual integral.
 *
 * Each of the integration routines can evaluate the polynomial
 * function, p, with a call to p(t).  See an example in the
 * integrateRiemann() function below.
 */


#include <stdio.h>
#include <stdlib.h>


typedef enum {			/* constants to define the integration method */
  RIEMANN, TRAPEZOID, SIMPSON
} IntegrationMethod;


/* global variables */

IntegrationMethod method;	/* the chosen method */
float rangeMin, rangeMax;	/* integration range */
int numIntervals;		/* number of intervals */
int numCoefficients;		/* number of polynomial coefficients */
float *coefficients;		/* array of coefficients, dynamically allocated (below) */


/* forward declarations of functions */

void getOptions( int argc, char **argv );
float integrate( float rangeMin, float rangeMax, int numIntervals );


/* Main program: Get the options; call the integrator; output the
 * result.
 */

int main( int argc, char **argv )

{
  float integral;

  getOptions( argc, argv );
  integral = integrate( rangeMin, rangeMax, numIntervals );
  printf( "computed integral = %g\n", integral );
}


/* Output a message about the usage of the program, then exit.
 */

void outputUsageAndQuit()

{
  fprintf( stderr, "usage: integrate method a b n coefficients\n" );
  fprintf( stderr, "     - `method' is one of: r (Riemann), t (trapezoid), or s (Simpson)\n" );
  fprintf( stderr, "     - `a' and `b' define the interval of integration\n" );
  fprintf( stderr, "     - `n' is the number of intervals\n" );
  fprintf( stderr, "     - `coefficients' define the polynomial to integrate, list in order\n" );
  fprintf( stderr, "     -                of decreasing power (e.g. t^2, then t^1, then t^0)\n" );
  exit(1);
}


/* Read the options from the command line
 */

void getOptions( int argc, char **argv )

{
  int i;

  if (argc < 6)
    outputUsageAndQuit();
  
  /* Get the method */

  switch (argv[1][0]) {
  case 'r':
    method = RIEMANN; break;
  case 't':
    method = TRAPEZOID; break;
  case 's':
    method = SIMPSON; break;
  default:
    outputUsageAndQuit();
  }

  /* Get the range */

  rangeMin = atof( argv[2] );
  rangeMax = atof( argv[3] );

  if (rangeMin >= rangeMax) {
    fprintf( stderr, "Error: Integration range [%g,%g] is reversed\n",
	     rangeMin, rangeMax );
    outputUsageAndQuit();
  }

  /* Get the number of intervals */

  numIntervals = atoi( argv[4] );
  if (numIntervals < 1) {
    fprintf( stderr, "Error: Number of intervals must be greater than zero\n" );
    outputUsageAndQuit();
  }

  /* Get the polynomial coefficients */

  numCoefficients = argc - 5;

  coefficients = (float *) malloc( numCoefficients * sizeof(float) ); /* allocate array */  

  for (i=0; i<numCoefficients; i++) /* fill array with coefficients */
    coefficients[numCoefficients-1-i] = atof( argv[i+5] );

  /* Output the options */

  printf( "method:     ");
  switch (method) {
  case RIEMANN:
    printf( "RIEMANN\n" ); break;
  case TRAPEZOID:
    printf( "TRAPEZOID\n" ); break;
  case SIMPSON:
    printf( "SIMPSON\n" ); break;
  }

  printf( "range:      [%g,%g]\n", rangeMin, rangeMax );
  printf( "intervals:  %d\n", numIntervals );

  printf( "polynomial: " );
  for (i=numCoefficients-1; i>=0; i--) {
    printf( "%g t^%d", coefficients[i], i );
    if (i > 0)
      printf( " + " );
  }
  printf( "\n" );
}


/* The polynomial function
 *
 * This function evaluates the polynomial at a particular point, t. 
 * It simply computes all the powers of t, multiplies them by the
 * coefficients, and sums them together.
 *
 * Note that there's another, more efficient way to do this.  It's
 * called "Horner's rule".  There's a small bonus mark if you look up
 * Horner's rule, implement it here, and report on the difference in
 * execution times (for several long computations) between Horner's
 * rule and the method below.  DON'T ATTEMPT THIS BONUS UNTIL YOU'VE
 * COMPLETED THE REST OF THE ASSIGNMENT!
 */

float p( float t )

{
  float power, total;
  int   i;

  power = 1;			/* the power is intially t^0 = 1 */
  total = 0;			/* sum of terms so far */

  /* sum up the terms */

  for (i=0; i<numCoefficients; i++) {
    total += coefficients[i] * power; /* add this term to total */
    power *= t;			/* increase power by one: t^i -> t^{i+1} */
  }

  return total;
}


/* The Riemann integrator
 *
 * YOUR COMMENTS HERE
 */

float integrateRiemann( float rangeMin, float rangeMax, int numIntervals )

{
  float integral;

  /* Here's an example of how to evaluate the polynomial function, p(t):
   * 
   * Delete this code and replace it with your own Riemann integrator.
   *
   * The following code evaluates the polynomial in the middle of the
   * interval and multiplies that value by the interval width.  This
   * is like a 1-interval-only Riemann integration.
   *
   * Note that we divide by 2.0, and not by 2.  If we used `2', C
   * would use integer division and would return an integer, rather
   * than a float.  This would (usually) *not* be the middle of the
   * interval!  Be careful to divide floats by floats and to divide
   * ints by ints!
   */

  integral = p( (rangeMin + rangeMax) / 2.0 ) * (rangeMax - rangeMin);

  return integral;
}


/* The trapezoid integrator
 *
 * YOUR COMMENTS HERE
 */

float integrateTrapezoid( float rangeMin, float rangeMax, int numIntervals )
     
{
  return 0;
}


/* The Simpson integrator
 *
 * YOUR COMMENTS HERE
 */

float integrateSimpson( float rangeMin, float rangeMax, int numIntervals )

{
  return 0;
}


/* The main integrator function, which calls the appropriate method
 * above.
 */


float integrate( float rangeMin, float rangeMax, int numIntervals )

{
  switch (method) {
  case RIEMANN:
    return integrateRiemann( rangeMin, rangeMax, numIntervals );
  case TRAPEZOID:
    return integrateTrapezoid( rangeMin, rangeMax, numIntervals );
  case SIMPSON:
    return integrateSimpson( rangeMin, rangeMax, numIntervals );
  }
}