#include <stdio.h>
#include "cquad.c"

/* Function pointer type. */
typedef double (*fptr) ( double , void * );


/* The test functions. */
double f1 ( double x , void *params ) {
    return exp(x);
    }
    
double f2 ( double x , void *params ) {
    return x >= 0.3;
    }

double f3 ( double x , void *params ) {
    return sqrt(x);
    }


double f4 ( double x , void *params ) {
    return (23.0/25) * cosh(x) - cos(x);
    }

double f5 ( double x , void *params ) {
    double x2 = x*x;
    return 1.0 / ( x2 * (x2 + 1) + 0.9);
    }

double f6 ( double x , void *params ) {
    return x * sqrt( x );
    }

double f7 ( double x , void *params ) {
    return 1.0 / sqrt(x);
    }

double f8 ( double x , void *params ) {
    double x2 = x*x;
    return 1.0 / (1 + x2*x2);
    }

double f9 ( double x , void *params ) {
    return 2.0 / (2 + sin(10*M_PI*x));
    }

double f10 ( double x , void *params ) {
    return 1.0 / (1 + x);
    }

double f11 ( double x , void *params ) {
    return 1.0 / (1 + exp(x));
    }

double f12 ( double x , void *params ) {
    return x / (exp(x) - 1.0);
    }

double f13 ( double x , void *params ) {
    return sin(100 * M_PI * x) / (M_PI * x);
    }

double f14 ( double x , void *params ) {
    return sqrt(50.0) * exp(-50*M_PI*x*x);
    }

double f15 ( double x , void *params ) {
    return 25.0 * exp(-25*x);
    }

double f16 ( double x , void *params ) {
    return 50 / M_PI * (2500 * x*x + 1);
    }

double f17 ( double x , void *params ) {
    double t1 = 50 * M_PI * x ,t2;
    t2 = sin(t1) / t1;
    return 50 * t2 * t2;
    }

double f18 ( double x , void *params ) {
    return cos( cos(x) + 3*sin(x) + 2*cos(2*x) + 3*sin(2*x) + 3*cos(3*x) );
    }

double f19 ( double x , void *params ) {
    return log(x);
    }

double f20 ( double x , void *params ) {
    return 1 / (x*x + 1.005);
    }

double f21 ( double x , void *params ) {
    return 1 / cosh( 10 * (x - 0.2) * 2 ) + 
          1 / cosh( 100 * (x - 0.4) * 4 ) + 
          1 / cosh( 1000 * (x - 0.6) * 8 );
    }

double f22 ( double x , void *params ) {
    return 4 * M_PI*M_PI * x * sin(20*M_PI*x) * cos(2*M_PI*x);
    }

double f23 ( double x , void *params ) {
    double t = 230*x - 30;
    return 1 / (1 + t*t);
    }

double f24 ( double x , void *params ) {
    return floor(exp(x));
    }

double f25 ( double x , void *params ) {
    return (x < 1) * (x + 1) + 
          (1 <= x && x <= 3) * (3 - x) + 
          (x > 3) * 2;
    }


/* The main function. */
int main ( int argc , char *argv[] ) {

    const static fptr funs[25] = { &f1 , &f2 , &f3 , &f4 , &f5 , &f6 , &f7 , 
        &f8 , &f9 , &f10 , &f11 , &f12 , &f13 , &f14 , &f15 , &f16 , &f17 ,
        &f18 , &f19 , &f20 , &f21 , &f22 , &f23 , &f24 , &f25 };
    const static double ranges[50] = { 0, 1 , 0, 1 , 0, 1 , -1, 1 , -1, 1 , 0, 1 , 0, 1 , 0, 1 , 0, 1 , 
        0, 1 , 0, 1 , 0, 1 , 0, 1 , 0, 10 , 0, 10 , 0, 10 , 0, 1 , 0, M_PI ,
        0, 1 , -1, 1 , 0, 1 , 0, 1 , 0, 1 , 0, 3 , 0, 5 };
    const static double f_exact[25] = { 1.7182818284590452354 , 0.7 , 2.0/3 , 0.4794282266888016674 , 
        1.5822329637296729331 , 0.4 , 2 , 0.86697298733991103757 , 
        1.1547005383792515290 , 0.69314718055994530942 , 0.3798854930417224753 , 
        0.77750463411224827640 , 0.49898680869304550249 , 
        0.5 , 1 , 0.13263071079267703209e+08 , 0.49898680869304550249 , 
        0.83867634269442961454 , -1 , 1.5643964440690497731 , 
        0.16349494301863722618 , -0.63466518254339257343 , 
        0.013492485649467772692 , 17.664383539246514971 , 7.5 };
    
    gsl_function F;
    double result, error;
    size_t neval;
    int fid, res;
    gsl_cquad_workspace *ws;
    
    
    /* Initialize the workspace */
    if ( (ws = gsl_cquad_workspace_alloc ( 200 )) == NULL ) {
        printf("error allocating workspace!\n");
        return 1;
        }
        
    /* Loop over the functions... */
    for ( fid = 0 ; fid < 25 ; fid++ ) {
    
        /* Initialize the function. */
        F.function = funs[fid];
        F.params = NULL;

        /* Call our quadrature routine. */
        res = gsl_cquad ( &F , ranges[2*fid] , ranges[2*fid+1] , 0.0 , 1.0e-12 , ws , &result , &error , &neval );

        /* Dump the output. */
        printf("f%i returned res=%i, result=%e, err=%e, neval=%i\n",
            fid+1, res, result, error, neval );
        printf("    relative error is %e (returned %e).\n",
            fabs( (result - f_exact[fid]) / f_exact[fid] ),
            fabs( error / f_exact[fid] ) );
        
        }
        
    /* All is well that ends well. */
    return 0;

    }