// function taking double variable to mpfr variable
void dobtomp(mpfr_t mp, double d) {
   int i;
   ostringstream strs;
   strs << scientific << setprecision(15) << d;
   mpfr_set_str(mp, strs.str().c_str(), 10, GMP_RNDN);
}

using namespace alglib;

// function to be integrated; bessel function of the first kind J_3(x)
void int_function_1_func(double x, double xminusa, double bminusx, double &y, void *ptr)
{
   mpfr_t my, mx;
   mpfr_inits2(200, mx, my, (mpfr_ptr)0);
   dobtomp(mx, x);
   mpfr_jn(my, 3, mx, GMP_RNDN);
   y=mpfr_get_d(my, GMP_RNDN);
}

// uses autogk to calculate integral
int main(int argc, char **argv)
{

    double a = 0;
    double b = 3;
    autogkstate s;
    double v;
    autogkreport rep;

    autogksmooth(a, b, s);
    alglib::autogkintegrate(s, int_function_1_func);
    autogkresults(s, v, rep);
   // prints result to screen
   printf("%.5f\n", double(v));
   
   // calculates analytical result and prints to screen
   mpfr_t m1, m2, m4;
   mpfr_inits2(200, m1, m2, m4, (mpfr_ptr)0);
   mpfr_set_str(m1, "3e0", 10, GMP_RNDN);
   mpfr_set_str(m4, "2e0", 10, GMP_RNDN);
   mpfr_jn(m2, 2, m1, GMP_RNDN);
   mpfr_mul(m2, m2, m4, GMP_RNDN);
   mpfr_j0(m1, m1, GMP_RNDN);
   mpfr_set_str(m4, "1e0", 10, GMP_RNDN);
   mpfr_sub(m4, m4, m1, GMP_RNDN);
   mpfr_sub(m4, m4, m2, GMP_RNDN);
   mpfr_out_str(stdout, 10, 10, m4, GMP_RNDN);
   cout << endl;
   mpfr_clears(m1, m2, m4, (mpfr_ptr)0);

    return 0;
}