/*
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));
}

inline REAL fprime(REAL x) { //derivatives
    //return 3*x*x;
    return (-1449 + 5*x*(4644 +25*x*(-1041 + 4*x*(646 + 25*x*(-29 + 12*x)))))/25000;
    //return exp(x);
    //return -10000.0*exp(-5000.0*(x-0.5)*(x-0.5))*(x-0.5)
}


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

int main(int argc,char *argv[])
{
    int N,i;
    REAL a,b,result,h,x,f0,f1,f2,fp,bdiff,fdiff,cdiff;
    
    //welcome info
    
    printf("# differentiation 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);
    
   
    x=a;
    h=(b-a)/N;
    f0=f(x);
    f1=f(x+h);
    for(i=1;i<N;i++) {
        x=a+i*h;
        f2=f(x+h);
        fp=fprime(x);
        
        bdiff=(f1-f0)/h;
        fdiff=(f2-f1)/h;
        cdiff=0.5*(f2-f0)/h;
        
        printf("%le %le %le %le %le %le\n",x,f1,fp,bdiff,fdiff,cdiff);
        
        
        f0=f1;
        f1=f2;
    }
    
    
    return 0;
    }

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