/* code to integerate a function, even with
a singularity (point where integrand blows up).
gsl "qags" routine has its output modified slightly
to make it c++ like. read gsl manual for
more info.
*/



#include <iostream>
#include <iomanip>
#include <cmath>
#include <gsl/gsl_integration.h>

using namespace::std;

// don't touch next line:     
double f (double x, void * params) {
       double f = pow(x,x);
       return f;
     }
     
int main ()
     {

       // don't touch:
       gsl_integration_workspace * w 
         = gsl_integration_workspace_alloc (1000);
       
       double result, error;

     
       gsl_function F;
       F.function = &f;

       //don't touch:
       gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000,
                             w, &result, &error); 
     
       cout << setprecision(10) << "result = \t" << result << endl;
       cout << setprecision(10) << "estimated error = \t" << error << endl;

       // the "w->size" is magic for now:

       cout << setprecision(6) << "intervals = \t" << w->size << endl;

       // don't touch:
     
       gsl_integration_workspace_free (w);
     
       return 0;
     }