#include <iostream>
#include <cmath>


extern"C"{
  int vegas_(double* region, int* ndim,double (*fcn)(double*) ,
	     int* init, int *ncall, int *itmx, int*nprn,
	     double *tgral, double *sd, double *chi2a);
}

double fcn(double p[10]);

int main(){

  FILE* outfile;
  outfile = fopen("vegas.dat","w");
  int dump;

  int ndim = 6, init = 0, ncall = 0, itmx = 10, nprn = -1;
  double result, sdr, chi, Pi = 3.1415926535, analytic = pow(Pi,3);
  double reg[12];

  reg[0] = -7;
  reg[1] = -7;
  reg[2] = -7;
  reg[3] = -7;
  reg[4] = -7;
  reg[5] = -7;
  reg[6] = 7;
  reg[7] = 7;
  reg[8] = 7;
  reg[9] = 7;
  reg[10] = 7;
  reg[11] = 7;

  for(int counter = 1; counter < 400; counter++){

    ncall = counter*100;

    dump = vegas_(reg,&ndim,fcn,&init,&ncall,&itmx,&nprn,&result,&sdr,&chi);
    fprintf(outfile,"%d \t%f\n",ncall,(result-analytic)/analytic);

  }

  return 0;
}
double fcn(double p[10]){

  double f;

  //Function definition valid for up to 10 dimensions
  f = exp(-pow(p[0],2))*exp(-pow(p[1],2))*exp(-pow(p[2],2))*exp(-pow(p[3],2))
    *exp(-pow(p[4],2))*exp(-pow(p[5],2))*exp(-pow(p[6],2))*exp(-pow(p[7],2))
    *exp(-pow(p[8],2))*exp(-pow(p[9],2));
  
  return f;
}