#include <math.h>
#include <stdio.h>
float fct(float x)
{
      return cos(x);
}
float integral(float a, int i, float h, int n);
int main(void)
{
/***********************************************************************
 *                                                                     *
 * This is a serial C-version of the MPI example on integration        *
 *                                                                     *
 ***********************************************************************/
  
      int n, p, i;
      float h, result, a, b, pi;

      pi = acos(-1.0);  /* = 3.14159... */
      a = 0.;           /* lower limit of integration */
      b = pi*1./2.;     /* upper limit of integration */
      p = 4;
      n = 500;          /* number of increment within each process */
      h = (b-a)/n/p;    /* length of increment */

      result = 0.0;
/* sum of integrals over all processes */
      for (i=0; i<p; i++)  {
        result += integral(a,i,h,n);
      }

      printf("The result =%f\n",result);
}
float integral(float a, int i, float h, int n)
{
      int j;
      float h2, aij, integ;

      integ = 0.0;                 /* initialize integral */
      h2 = h/2.;
      for (j=0;j<n;j++) {             /* sum over all "j" integrals */
        aij = a + (i*n + j)*h;        /* lower limit of "j" integral */
        integ += fct(aij+h2)*h;
      }  
      return integ;
}