#include <math.h>
#include <stdio.h>
int num_dim = 1;
extern long random();
double acc = -1, alph = 1.5, xl[10], xu[10];
int bins_max = 50, mode, prn = 1, max_it_num = 5, calls = 1000;
static int  it_start, bins_prev, sub_calls, it_num, bins, boxes;
static double jac, dxg, wtd_int_sum, sum_wgts, chi_sum, vol;
static double delx[10], grid_sum[51][10], bin_sum[51][10], y_bin[51][10];
 
/* The following variables are needed for the print functions. */
 
//FILE *o_file, *fopen();
FILE *o_file;
 
/*  VEGAS: Performs N-dimensional Monte-Carlo Integration
Created by G. P. Lepage in Sept 1976/(rev)Aug 1979
Algorithm described in J. Comp. Phys. 27, 192(1978)
Rewritten into C by D. R. Yennie, September 1984     */

//void vegas1(), vegas2(), vegas3(); 
 
/*  Description of coordinates:  x[j] is the j-th variable of integration,
with upper and lower limits xu[j] and xl[j].  The integration length in the
j-th direction is delx[j].  In case the actual range of integration is
infinite, it must first be mapped into a finite range by an appropriate
transformation.  Each coordinate x[j] is mapped linearly to a variable y[j] in
the range 0 to 1. This range is divided into bins with boundaries y_bin[i][j],
where i=0 corresponds to y=0 and i=bins to y=1.  The integral contribution
from the bin preceding y_bin[i][j] is called bin_sum[i][j].  Another similar
sum which is used to refine the grid is called grid_sum[i][j].  A third
parameter used in defining the real coordinate using random numbers is called
z.  It ranges from 0 to bins.  Its integer part gives the lower index of the
bin into which a call is to be placed, and the  remainder gives the location
inside the bin.  */
 
void refine ( double y_bin[51][10], double r[51], double rc, int j, int n1, int n2 );
void change_kg ( int kg[10], int ng, int j );
void init_array( double array[51][10], int n1,int n2 );
void prn_lim( double a[],double b[],int m );
void prn_head( int a,int b,int c,int d,double e,int f,double g,int h,int i ,int j );
void prn_res ( int a, double b, double c, double d, double e, double f );
void prn_grid ( double y[25][10],double s[25][10],int m,int n,int p );
void vegas1( double (*fxn)(double*), double* tot_int_ptr, double* tot_sig_ptr, double* chi_sq_ptr );
void vegas2( double (*fxn)(double*), double* tot_int_ptr, double* tot_sig_ptr, double* chi_sq_ptr );
void vegas3( double (*fxn)(double*), double* tot_int_ptr, double* tot_sig_ptr, double* chi_sq_ptr );
void vegas ( double (*fxn)(double*), double* tot_int_ptr, double* tot_sig_ptr, double* chi_sq_ptr )
  {
    int j;
    o_file = fopen("vegas.out", "w");

    init_array( y_bin, 50, num_dim );
    for ( j = 1 ; j <= num_dim ; ++j )
        y_bin[1][j] = 1.;
    bins_prev= 1;
    vol = 1;
    for ( j = 1 ; j <= num_dim ; ++j )
      {
        delx[j] = xu[j]-xl[j];
        vol *= delx[j];
      }
    prn_lim ( xl,xu,num_dim );
    vegas1( fxn, tot_int_ptr, tot_sig_ptr, chi_sq_ptr );
  }
 
void vegas1( double (*fxn)(double*), double* tot_int_ptr, double* tot_sig_ptr, double* chi_sq_ptr )
  {
    wtd_int_sum= 0;
    sum_wgts= 0;
    chi_sum= 0;
    it_num= 1;
    vegas2( fxn, tot_int_ptr, tot_sig_ptr, chi_sq_ptr );
  }
 
void vegas2( double (*fxn)(double*), double* tot_int_ptr, double* tot_sig_ptr, double* chi_sq_ptr )
  {
 
    /*  At this point, the calculation is divided into three possible cases:
    mode = 0 : obtained only by choice at the beginning.  VEGAS uses
    importance sampling: the points are distributed over the whole grid
    randomly, but the bins adjust to emphasize the regions where the
    function is largest.
    mode = 1 or -1 as decided by a test (based on the number of calls,
    dimension, and maximum number of intervals).  The +1 choice
    concentrates increments where the integrand is largest and the -1
    choice where the error is largest.  In addition to the bin structure,
    the volume is divided into boxes, with a certain number of calls
    ("sub_calls") going into each box.  In both cases, the number of calls
    is given by calls = sub_calls * boxes ** num_dim and is usually a bit
    smaller than the nominal number specified.  For +1, the number of
    bins in each dimension ("bins") is "bins_max", while the number of
    boxes in each dimension ("boxes") is smaller.  Thus, each box contains
    several (not necessarily an integral number!) bins.  For -1, "bins" is
    an integer multiple of "boxes", and "bins" is smaller than "bins_max".
    */
    int i, j, k;
    double one = 1., two = 2.;
    bins= bins_max;
    boxes= 1;
    if ( mode != 0 )
      {
        boxes = pow( calls/two , one/num_dim )+.00001;
        mode = 1;
        if ( (2*boxes-bins_max) >= 0 )
          {
            mode = -1;
            i = boxes/bins_max + 1;
            bins = boxes/i;
            boxes = i*bins;
          }
      }
    k = pow( one*boxes , one*num_dim )+.00001;
    sub_calls = calls/k;
    if ( sub_calls<2 )
        sub_calls= 2;
    calls = sub_calls*k;
    dxg = one*bins/boxes;
    jac = one/calls;
    jac = jac * vol * pow ( one*bins , one*num_dim );
    /*         jac is the total volume of x-space divided by the average number
    of calls per bin.   */
 
    /* If the number of bins changes, bins are expanded or contracted
    accordingly, while preserving bin density */
 
    if ( bins != bins_prev )
      {
        double r[51], rc;
        rc = one*bins_prev/bins;                   /*ratio of bin sizes*/
        for ( j = 1 ; j <= num_dim ; ++j )
          {
            for ( i= 1 ; i <= bins_prev ; ++i )
                r[i] = 1;
            refine ( y_bin, r, rc, j, bins_prev, bins );
          }
        bins_prev = bins;
      }
    if ( prn >= 0 )
      {
        int a, b, c, d;
        a = num_dim;
        b = calls;
        c = it_num;
        d = max_it_num;
        prn_head ( a, b, c, d, acc, prn, alph, mode, bins, boxes );
      }
    vegas3( fxn, tot_int_ptr, tot_sig_ptr, chi_sq_ptr );
  }
 
void vegas3( double (*fxn)(double*), double* tot_int_ptr, double* tot_sig_ptr, double* chi_sq_ptr )
{
 
    /* Start of iteration  */
 
    int i, j, k, q, p;
    double cum_int, cum_sig, chi_sq;
    it_start = it_num;
    cum_int= 1.;
    cum_sig= 1.;
    q = cum_int;
    p = max_it_num;
    //for ( it_num = it_start ; it_num <= p && acc*fabs( q ) < cum_sig ; ++it_num )
    for ( it_num = it_start ; it_num <= p && acc*abs( q ) < cum_sig ; ++it_num )
      {
        double intgrl, intgrl_sq, sig;
        double wgt, x[10];
        int kg[10];
        for ( j= 0 ; j <= num_dim ; ++j )
            kg[j] = 0;
        init_array( grid_sum, bins, num_dim );
        init_array( bin_sum, bins, num_dim );
        intgrl= 0.;
        sig= 0.;
        for ( ; kg[0] == 0 ; change_kg ( kg, boxes, num_dim ) )
          {
            int ia[10];
            double  f, f_sq, f_sum, f_sq_sum;
            /*Note: rand returns an integer between 0 and 2**31-1=2147483647 ! */
            for ( k= 1, f_sum = 0., f_sq_sum = 0. ; k <= sub_calls ; ++k )
              {
                double jacbin, y, z;
                jacbin = jac;
                for ( j= 1 ; j <= num_dim ; ++j )
                  {
                    double binwdth;
                    z = ( kg[j]+random()/2147483647. )*dxg;
                    ia[j] = z;
                    if ( ia[j]==0 )
                      {
                        binwdth= y_bin[1][j];
                        y = z*binwdth;
                      }
                    else
                      {
                        binwdth = y_bin[ia[j]+1][j] - y_bin[ia[j]][j];
                        y = y_bin[ia[j]][j] + (z - ia[j])*binwdth;
                      }
                    x[j] = xl[j] + y*delx[j];
                    jacbin = jacbin*binwdth;
                  }
                /* At this point, jacbin is the real volume of the bin/(ave. no. calls/bin) */
                f = jacbin;
                f *= fxn( x );
                f_sq = f * f;
                f_sum += f;
                f_sq_sum += f_sq;
                for ( j= 1 ; j<= num_dim ; ++j )
                  {
                    bin_sum[ia[j]+1][j] += f;
                    if ( mode != -1 )
                        grid_sum[ia[j]+1][j] += f_sq;
                  }
              }                                   /* end of k loop */
            f_sq_sum = f_sq_sum*sub_calls - f_sum*f_sum;
            intgrl += f_sum;
            sig += f_sq_sum;
            if ( mode == -1 )
              {
                for ( j= 1 ; j <= num_dim ; ++j )
                    grid_sum[ia[j]+1][j] += f_sq_sum;
              }
          }                                      /* end of kg loop */
        /* Compute final results for this iteration   */
 
        sig = sig/(sub_calls-1);
        intgrl_sq = intgrl*intgrl;
        wgt = 1/sig;
        wtd_int_sum += intgrl*wgt;
        sum_wgts += wgt;
        chi_sum += intgrl_sq*wgt;
        cum_int = wtd_int_sum/sum_wgts;
        chi_sq = 0.;
        if ( it_num != 1 )
            chi_sq = (chi_sum-wtd_int_sum*cum_int)/(it_num-1.);
        cum_sig = sqrt( 1/sum_wgts );
        if ( prn >= 0 )
          {
            sig = sqrt( sig );
 
            /*  Print out the answer    */
 
            prn_res ( it_num, intgrl, sig, cum_int, cum_sig, chi_sq );
            if ( it_num == max_it_num && prn > 0)
                prn_grid ( y_bin, bin_sum, num_dim, bins, prn );
          }
 
        /*   Refine grid  */
 
        for ( j = 1 ; j <= num_dim ; ++j )
          {
            double old, xnew, grid_tot, r[51], rc;
            old = grid_sum[1][j];
            xnew = grid_sum[2][j];
            grid_sum[1][j] = (old+xnew)/2;
            grid_tot = grid_sum[1][j];
            for ( i = 2 ; i < bins ; ++i )
              {
                grid_sum[i][j] = old+xnew;
                old = xnew;
                xnew = grid_sum[i+1][j];
                grid_sum[i][j] = (grid_sum[i][j]+xnew)/3;
                grid_tot += grid_sum[i][j];
              }
            grid_sum[bins][j] = (xnew+old)/2;
            grid_tot += grid_sum[bins][j];
            for ( i = 1, rc = 0 ; i<= bins ; ++i )
              {
                r[i] = 0;
                if ( grid_sum[i][j] > 0 )
                  {
                    old = grid_tot/grid_sum[i][j];
                    r[i] = pow ( ( (old-1)/old/log( old ) ), alph );
                  }
                rc += r[i];
              }
            rc /= bins;
            refine ( y_bin, r, rc, j, bins, bins );
          }
      }
    *tot_int_ptr = cum_int;
    *tot_sig_ptr = cum_sig;
    *chi_sq_ptr  = chi_sq;
  }
 
void refine ( double y_bin[51][10], double r[51], double rc, int j, int n1, int n2 )
  {
    int i, k;
    double xold, xnew, xin[51], dr;
    xnew= 0;
    dr= 0;
    i= 1;
    for ( k= 1 ; k <= n1 ; ++k )
      {
        dr += r[k];
        xold = xnew;
        xnew = y_bin[k][j];
        for ( ; dr > rc ; ++i )
          {
            dr = dr-rc;
            xin[i] = xnew-(xnew-xold)*dr/r[k];
          }
      }
    for ( i = 1 ; i < n2 ; ++i )
        y_bin[i][j] = xin[i];
    y_bin[n2][j] = 1;
    return;
  }
 
void change_kg ( int kg[10], int ng, int j )
  {
    ++kg[j];
    if ( j == 0 )
        return;
    if ( kg[j] % ng == 0 )
      {
        kg[j] = 0;
        change_kg ( kg, ng, j-1 );
      }
  }
 
void init_array( double array[51][10], int n1,int n2 )
  {
    int i ,j;
    for ( j= 1 ; j <= n2 ; ++j )
      {
        for ( i = 1 ; i <= n1 ; ++i )
            array[i][j] = 0;
      }
  }
/*   Set of programs for printing vegas results   */
 
void prn_lim( double a[],double b[],int m )
  {
    int j;
    fprintf ( o_file, "The limits of integration are:\n" );
    for ( j = 1 ; j <= m ; ++j )
        fprintf ( o_file,"\nxl[%d]=%f    xu[%d]=%f", j, a[j], j, b[j] );
    fprintf ( o_file,"\n" );
  }
 
void prn_head( int a,int b,int c,int d,double e,int f,double g,int h,int i ,int j )
  {
    fprintf(o_file, "\nnum_dim=%d,  calls=%d,  it_num=%d,  max_it_num=%d,  acc=%f,\
            \n", a, b, c, d, e );
    fprintf(o_file,  "   prn=%d,  alph=%f,  mode=%d,  bins=%d  boxes=%d\n",f, g, h, i, j );
    fprintf(o_file,  "\n\
            single.......iteration             accumulated......results   \n\
iteration    integral      sigma       integral     sigma     chi-sq            /it\n\
            number \n" );
 
    fprintf( o_file, "\nnum_dim=%d,  calls=%d, it_num=%d,\
             max_it_num=%d, acc=%f,\n\
             prn=%d,  alph=%f,  mode=%d\
             bins=%d  boxes=%d\n\
             \n       single.......iteration             \
             accumulated......results   \n\
             iteration   integral     sigma       integral       sigma\
             chi-sq/it\n    number \n",
             a, b, c, d, e, f, g, h, i, j );
  }
 
void prn_res ( int a, double b, double c, double d, double e, double f )
  {
    fprintf( o_file, "%6d%17.6e%10.2e%16.6e%10.2e%12.2e\n", a, b, c, d, e,
             f );
  }
 
void prn_grid ( double y[25][10],double s[25][10],int m,int n,int p )
  {
    int mod, i, j;
    for ( j = 1 ; j <= m ; ++j )
      {
        fprintf( o_file, "\n axis %d \n", j );
        fprintf( o_file, "      x          delta         x     " );
        fprintf( o_file, "    delta         x          delta\n" );
        for ( i= 1+p/2, mod = 1 ; i <= n ; i += p, ++mod )
          {
            fprintf ( o_file, "%11.2e%13.2e ", y[i][j], s[i][j] );
            if ( mod%3 == 0 )
                fprintf ( o_file, "\n" );
          }
        fprintf ( o_file, "\n" );
      }
  }