/*
Demo program for numerical integration
 
Use "REAL" as floating point type instead of double or float
 
Compile with optional flag:
    -DSINGLE to use float instead of double
*/

#include <stdio.h>
#include <stdlib.h> /* for rand() */
#include <unistd.h> /* for getpid() */
#include <time.h> /* for time() */
#include <math.h>
#include <assert.h>

// construct REAL "type," depending on desired precision

#ifdef SINGLE
#define REAL float
#else
#define REAL double
#endif

// ******************************************************


inline REAL f(REAL x) {
    //return x*x*x;
    return (x - 0.1)*(x - 0.2)*(x - 0.3)*(x - 0.6)*(x - 0.8)*(x - 0.9);
    //return (exp(x)+1.0);
    //return exp(-5000.0*(x-0.5)*(x-0.5));
}


REAL rectangular_left(int N,REAL a,REAL b) {
    REAL h,sum;
    int i;
    
    //interval length
    h=(b-a)/N;
    
    //"left" sum
    sum=0.0;
    for(i=0;i<N;i++) {
        sum+=f(a+i*h);
    }
    
    //multiply the interval length
    return sum*h;
}

REAL rectangular_left_kahan(int N,REAL a,REAL b) {
    REAL h,f0,f1,sum;
    REAL c=0.0,y,t;
    int i;
    
    //interval length
    h=(b-a)/N;
    
    //"left" sum
    sum=0.0;
    for(i=0;i<N;i++) {
        y=f(a+i*h)-c;//function value - carryover
        t=sum+y;
        c=(t-sum)-y;
        sum=t;
    }
    
    //multiply the interval length
    return sum*h;
}

REAL rectangular_right(int N,REAL a,REAL b) {
    REAL h,sum;
    int i;
    
    //interval length
    h=(b-a)/N;
    
    //"left" sum
    sum=0.0;
    for(i=1;i<=N;i++) {
        sum+=f(a+i*h);
    }
    
    //multiply the interval length
    return sum*h;
}

REAL rectangular_right_kahan(int N,REAL a,REAL b) {
    REAL h,f0,f1,sum;
    REAL c=0.0,y,t;
    int i;
    
    //interval length
    h=(b-a)/N;
    
    //"left" sum
    sum=0.0;
    for(i=1;i<=N;i++) {
        y=f(a+i*h)-c;//function value - carryover
        t=sum+y;
        c=(t-sum)-y;
        sum=t;
    }
    
    //multiply the interval length
    return sum*h;
}

REAL rectangular_mid(int N,REAL a,REAL b) {
    REAL h,sum;
    int i;
    
    //interval length
    h=(b-a)/N;
    
    //"left" sum
    sum=0.0;
    for(i=0;i<N;i++) {
        sum+=f(a+(i+0.5)*h);
    }
    
    //multiply the interval length
    return sum*h;
}

REAL rectangular_mid_kahan(int N,REAL a,REAL b) {
    REAL h,f0,f1,sum;
    REAL c=0.0,y,t;
    int i;
    
    //interval length
    h=(b-a)/N;
    
    //"left" sum
    sum=0.0;
    for(i=0;i<N;i++) {
        y=f(a+(i+0.5)*h)-c;//function value - carryover
        t=sum+y;
        c=(t-sum)-y;
        sum=t;
    }
    
    //multiply the interval length
    return sum*h;
}


REAL trapezoid(int N,REAL a,REAL b) {
    REAL h,f0,f1,sum;
    int i;
    
    //interval length
    h=(b-a)/N;
    
    //trapezoid sum
    sum=0.0;
    f0=f(a);
    for(i=1;i<=N;i++) {
        f1=f(a+i*h);
        sum+=0.5*(f0+f1); //add midpoint value
        f0=f1;
    }
    
    //multiply the interval length
    return sum*h;
}

REAL trapezoid_kahan(int N,REAL a,REAL b) {
    REAL h,f0,f1,sum;
    REAL c=0.0,y,t;
    int i;
    
    //interval length
    h=(b-a)/N;
    
    //trapezoid sum
    sum=0.0;
    f0=f(a);
    for(i=1;i<=N;i++) {
        f1=f(a+i*h);
        y=0.5*(f0+f1)-c;//midpoint value - carryover
        t=sum+y;
        c=(t-sum)-y;
        sum=t;
        f0=f1;
    }
    
    //multiply the interval length
    return sum*h;
}


REAL simpson(int N,REAL a,REAL b) {
    REAL h,f0,f1,sum;
    int i;
    
    //interval length
    h=(b-a)/N;
    
    //simpson rule
    sum=0.0;
    f0=f(a);
    for(i=1;i<=N;i++) {
        f1=f(a+i*h);
        sum+=(f0+4*f(a+(i-0.5)*h)+f1); //simpson values
        f0=f1;
    }
    
    //multiply the interval length
    return sum*h/6.0;
}

//coordiantes & weights for Gauss-Legendre 5 point integration
static REAL GLx[] ={-0.90617984593866399280,
    -0.53846931010568309104,
    0,
    0.53846931010568309104,
    0.90617984593866399280};
static REAL GLw[] ={0.23692688505618908751,
    0.47862867049936646804,
    0.56888888888888888889,
    0.47862867049936646804,
    0.23692688505618908751};

REAL GL_5(int N,REAL a,REAL b) {
    REAL x,h,hh,sum;
    int i,j;
    
    //interval length
    h=(b-a)/N;
    hh=0.5*h;
    
    //Gauss-Legendre 5-point method
    sum=0.0;
    for(i=0;i<N;i++) {
        for(j=0;j<5;j++) {
			x=a+i*h+hh*(1.0+GLx[j]); //Gauss points
			sum+=GLw[j]*f(x); //multiply weights
		}
    }
    
    //multiply the half the interval length
    return sum*hh;
}

// ******************************************************

int main(int argc,char *argv[])
{
    int N;
    REAL a,b,result;
    
    //welcome info
    
    printf("integration demo using ");
#ifdef SINGLE
    printf("single");
#else
    printf("double");
#endif
    printf(" precision arithmetics.\n");
    
    printf("Optinal parameters are: discretization intervals (N), lower limit (a), upper limit (b)\n");
    
    //default parameters
    
    N=5;
    a=0.0;
    b=1.0;

    // check if arguments are present and read them
    
    if(argc > 1 ) N = atoi(argv[1]);
    if(argc > 2 ) a = atof(argv[2]);
    if(argc > 3 ) b = atof(argv[3]);

    //excute
   
    printf("using parameters: N=%d, a=%.16le, b=%.16le\n",N,a,b);
    
    result=rectangular_left(N,a,b);
    printf("left rectangular method:             %.16le\n",result);
    result=rectangular_left_kahan(N,a,b);
    printf("left rectangular method + kahan:     %.16le\n",result);
    
    result=rectangular_right(N,a,b);
    printf("right rectangular method:            %.16le\n",result);
    result=rectangular_right_kahan(N,a,b);
    printf("right rectangular method + kahan:    %.16le\n",result);
    
    result=rectangular_mid(N,a,b);
    printf("midpoint rectangular method:         %.16le\n",result);
    result=rectangular_mid_kahan(N,a,b);
    printf("midpoint rectangular method + kahan: %.16le\n",result);
    
    result=trapezoid(N,a,b);
    printf("trapezoid method:                    %.16le\n",result);
    result=trapezoid_kahan(N,a,b);
    printf("trapezoid method + kahan:            %.16le\n",result);
    
    result=simpson(N,a,b);
    printf("simpson method:                      %.16le\n",result);
    
    result=GL_5(N,a,b);
    printf("Gauss-Legendre 5-point method:       %.16le\n",result);
    
    return 0;
    }

// ******************************************************