/* An implementation of a Serial Simpsons' Rule
 *
 * by John Weatherwax
 *
 * Input: 
 *    a, b: limits of integration.
 *    n: number of trapezoids.
 * Output:  Estimate of the integral from a to b of f(x) 
 *    using the trapezoidal rule and n trapezoids.
 *
 * Notes:  
 *    1.  f(x) is hardwired.
 *
 * See Chap. 4, pp. 60 & ff in PPMPI.
 */
#include <stdio.h>

main(int argc, char** argv) {
    float       a;         /* Left endpoint             */
    float       b;         /* Right endpoint            */
    int         n;         /* Number of trapezoids      */
    float       h;         /* Trapezoid base length     */

    float       integral;  /* Integral */

    /* Local declarations */ 
    float Simp(float local_a, float local_b, int local_n, float h);

    printf("Enter a, b, and n\n");
    scanf("%f %f %d", &a, &b, &n);

    h = (b-a)/n;    /* h is the same for all processes */
    integral = Simp(a,b,n,h);

    /* Print the result */
    printf("With n = %d Simpson's, our estimate\n", n);
    printf("of the integral from %f to %f = %18.10f\n", a, b, integral); 

} /*  main  */

/********************************************************************/
float Simp(
          float  local_a   /* in */, 
          float  local_b   /* in */, 
          int    local_n   /* in */, 
          float  h         /* in */) { 

    float integral;   /* Store result in integral  */ 
    float x; 
    int i; 

    float f(float x); /* function we're integrating */

    integral = f(local_a) + f(local_b);
    x = local_a; 
    for (i = 1; i <= local_n-1; i++) { 
        x = x + h; 
	/* METHOD #1: */ 
	if( (i % 2) == 1 ){ /* every ODD  term gets a multiplication by 4 */
	  integral = integral + 4*f(x); 
	}else{              /* every EVEN term gets a multiplication by 2 */
	  integral = integral + 2*f(x); 
	}
	/* VS METHOD #2. ... this should be faster ... add a second f(x) if needed */ 
	/* integral = integral + 2*f(x); 
	   if( (i % 2) == 1 ) integral = integral + 2*f(x); */ 
    } 
    integral = (integral*h)/3.0; 
    return integral;
} /*  Trap  */


/********************************************************************/
float f(float x) { 
    float return_val; 
    /* Calculate f(x). */
    /* Store calculation in return_val. */ 
    return_val = x*x;
    return return_val; 
} /* f */