#include <cmath>
#include <iostream>

#include "cubaint.hpp"

int test_integrand (const int *ndim, const double x[],
		    const int *ncomp, double f[], void *userdata)
{
  for (int icomp=0; icomp<*ncomp; icomp++) {  
    f[icomp] = 1;
    for (int idim=0; idim<*ndim; idim++)
      f[icomp] *= exp(-pow(x[idim],2));
  }

  return 0;
}

int test_integrand_c (const int *ndim, const double x[],
		      const int *ncomp, std::complex<double> f[], void *userdata)
{
  for (int icomp=0; icomp<*ncomp; icomp++) {  
    f[icomp] = 1.;
    for (int idim=0; idim<*ndim; idim++)
      f[icomp] *= exp(-pow(x[idim],2));
  }
  
  return 0;
}

int main ()
{
  int nregions, neval, fail;
  double error, prob, integral;

  std::cout << "int_0^1 dx e-x^2: " << std::endl;
  for (int itype=0; itype<4; itype++) {
    cubaint::integrate (&test_integrand, cubaint::intmethod(itype),
			{{0,1}}, 1,
			nregions,
			neval, fail, &integral, 
			&error, &prob, 0);
    std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl;
  }

  std::cout << "int_0^INF dx e-x^2: " << std::endl;
  for (int itype=0; itype<4; itype++) {
    cubaint::integrate (&test_integrand, cubaint::intmethod(itype),
			{{0,cubaint::INF}}, 1,
			nregions,
			neval, fail, &integral, 
			&error, &prob, 0);
    std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl;
  }

  std::cout << "int_-INF^0 dx e-x^2: " << std::endl;
  for (int itype=0; itype<4; itype++) {
    cubaint::integrate (&test_integrand, cubaint::intmethod(itype),
			{{cubaint::INF,0}}, 1,
			nregions,
			neval, fail, &integral, 
			&error, &prob, 0);
    std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl;
  }

  std::cout << "int_-INF^INF dx e-x^2: " << std::endl;
  for (int itype=0; itype<4; itype++) {
    cubaint::integrate (&test_integrand, cubaint::intmethod(itype),
			{{cubaint::INF,cubaint::INF}}, 1,
			nregions,
			neval, fail, &integral, 
			&error, &prob, 0);
    std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl;
  }

  std::cout << "int_1^INF dx e-x^2: " << std::endl;
  for (int itype=0; itype<4; itype++) {
    cubaint::integrate (&test_integrand, cubaint::intmethod(itype),
			{{1,cubaint::INF}}, 1,
			nregions,
			neval, fail, &integral, 
			&error, &prob, 0);
    std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl;
  }

  std::cout << "int_-INF^1 dx e-x^2: " << std::endl;
  for (int itype=0; itype<4; itype++) {
    cubaint::integrate (&test_integrand, cubaint::intmethod(itype),
			{{cubaint::INF,1}}, 1,
			nregions,
			neval, fail, &integral, 
			&error, &prob, 0);
    std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl;
  }

  std::complex<double> error_c, prob_c, integral_c;
  std::cout << "int_0^1 dx (e-x^2+i0): " << std::endl;
  for (int itype=0; itype<4; itype++) {
    cubaint::integrate (&test_integrand_c, cubaint::intmethod(itype),
			{{0,1}}, 1,
			nregions,
			neval, fail, &integral_c, 
			&error_c, &prob_c, 0);
    std::cout << integral_c << " +- " << error_c << " (" << itype << ")" << std::endl;
  }

}