#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef double (*DFCN)(double *);
double li2(double);
double li3(double);
void vegas(DFCN fcn,double accuracy,int dimension,long ncalls,
	int itmax,int printflag);
double ranf(long dummy);
double ipow(double x,int y);
int iipow(int x,int y);
double f(double *x);

int main() {
	vegas(f,1.0E-08,1,100000,50,1);
	return(0);
}

double f(double *y)
{
	double x = y[0];
	return(log(1-x)*li2((1+x)*0.5)*li3((1-x)/(1+x))/(1+x));
}

static double a2[7] = {
	- 0.027777777777777777777,
	  0.000277777777777777778,
	- 0.000004724111866969010,
	  0.000000091857730746620,
	- 0.000000001897886998897,
	  0.000000000040647616451,
	- 0.000000000000892169102
/*
	  0.000000000000019939296,
	- 0.000000000000000451898
*/
};
static double li2of1 = 1.64493406684822643647;

double li2(double x)
/*
	The dilogarithm in the range -infty to +1
	Made by J.Vermaseren 26-feb-1992
	Rewritten into C 5-oct-1996
	a2 is generated with
		a2[0] = -1./36.;
		a2[1] = -1./100.*a2[0];
		a2[2] = -5./294.*a2[1];
		a2[3] = -7./360.*a2[2];
		a2[4] = -5./242.*a2[3];
		a2[5] = -7601./354900.*a2[4];
		a2[6] = -91./4146.*a2[5];
		a2[7] = -3617./161840.*a2[6];
		a2[8] = -3728695./164522862.*a2[7];
*/
{
	double y,z,u,ll;
	if ( x < -1. ) {
		y = log(1.-1./x);
		z = y*y;
		u = log(-x);
		ll = -(((((((a2[6]*z+a2[5])*z+a2[4])*z+a2[3])*z+a2[2])*z
			+a2[1])*z+a2[0])*z-1.)*y+0.25*z-0.5*u*u-li2of1;
	}
	else if ( x > 0.5 ) {
		if ( x > 1. ) {
			printf("li2 will not deal with arguments greater than 1\n");
			exit(-1);
		}
		if ( x == 1. ) ll = li2of1;
		else {
			u = log(1.-x);
			y = log(x);
			z = y*y;
			ll = li2of1-u*y
				-(((((((a2[6]*z+a2[5])*z+a2[4])*z+a2[3])*z+a2[2])*z
				+a2[1])*z+a2[0])*z-1.)*y+0.25*z;
		}
	}
	else {
		y = log(1.-x);
		z = y*y;
		ll = (((((((a2[6]*z+a2[5])*z+a2[4])*z+a2[3])*z+a2[2])*z
				+a2[1])*z+a2[0])*z-1.)*y-0.25*z;
	}
	return(ll);
}

static double a[13] = {
	 .99999999999999999912,  -.37500000000000040999,   .07870370370370388389,
	-.00868055555553163783,   .00012962962962362201,   .00008101851811970592,
	-.00000341935708570435,  -.00000132865361060450,   .00000008660826954690,
	 .00000002525094774351,  -.00000000214332111582,  -.00000000049736302346,
	 .00000000005039955318 };
static double b[10] = {
	1.00000000000000014921, .12500000000000063337, .03703703703629146347,
	 .01562499999873420964, .00800000059586145647, .00462963033767957726,
	 .00291528533212156926, .00195296346138406238,
	 .00139071443831508537, .00101609376397330399 };
static double li3of1 = 1.20205690315959428540;
static double pi6   = 1.64493406684822643647;

double li3(double x)
{
/*
	Trilog in the range -infty to +1
	Made by J.Vermaseren 15-mar-1992
	Rewritten to C: 5-oct-1996
	The series has been rewritten in terms of -ln(1-x) and then it
	has been improved with a Chebychev expansion. This gives 13 terms
	instead of 17. Result is good to 1.E-16
*/
	double y,z,l3;
	if ( x < -1. ) {
/*
		Use li3(x) = li3(1/x)-log(-x)*pi6-(log(-x))**3/6
*/
		z = log(-x);
		y = -log(1.-1./x);
		l3 = y*(a[0]+y*(a[1]+y*(a[2]+y*(a[3]+y*(a[4]+y*(a[5]+y*(a[6]
			+y*(a[7]+y*(a[8]+y*(a[9]+y*(a[10]+y*(a[11]+y*a[12]
			))))))))))))-z*(pi6+z*z/6.);
	}
	else if ( x > 0.5 ) {
		if ( x >= 1. ) {
			if ( x == 1. ) l3 = li3of1;
			else {
				printf("li3 cannot deal with x > 1\n");
				exit(-1);
			}
		}
		else {
/*
			Use Li3(x) = -Li3(1-x)-Li3(1-1/x)+Li3(1)
			             +log(x)*(pi6+(log(x))**2/6-.5*log(x)*log(1-x))
			The odd terms in the trilogs cancel!
*/
			y = log(x);
			z = y*y;
			l3 = -2.*z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))))
				+y*(pi6+z/6.-y*log(1.-x)*0.5)+li3of1;
		}
	}
	else if ( fabs(x) <= 0.1 ) {
		l3 = x*(b[0]+x*(b[1]+x*(b[2]+x*(b[3]+x*(b[4]+x*(b[5]+x*(b[6]
			+x*(b[7]+x*(b[8]+x*b[9])))))))));
	}
	else {
		y = -log(1.-x);
		l3 = y*(a[0]+y*(a[1]+y*(a[2]+y*(a[3]+y*(a[4]+y*(a[5]+y*(a[6]
			+y*(a[7]+y*(a[8]+y*(a[9]+y*(a[10]+y*(a[11]+y*a[12]))))))))))));
	}
	return(l3);
}

#define MAXDIM 10
#define MAXDIV 500

int ndo,it;
double si,si2,swgt,schi,scalls;
double xi[MAXDIV][MAXDIM],d[MAXDIV][MAXDIM],di[MAXDIV][MAXDIM];
int nxi[MAXDIV][MAXDIM];
double avgi,sd,chi2a;

void
vegas(
	DFCN fcn,          /* The function to be integrated */
	double accuracy,   /* The desired relative accuracy */
	int dimension,     /* Number of dimensions */
	long ncalls,       /* Number of points per iteration */
	int itmax,         /* Maximum number of iterations */
	int printflag)     /* How much output */
{
	int i,j,k,ndmx,nd,ng,npg,ndm,iaj,iaj1,mds;
	double calls,dxg,dv2g,xnd,xjac,rc,dr,xn,xo,ti,tsi,fb,f2b,wgt,f,f2,ti2;
	double xin[MAXDIV],r[MAXDIV],dx[MAXDIM],dt[MAXDIM];
	int ia[MAXDIM],kg[MAXDIM];
	double xl[MAXDIM],xu[MAXDIM],qran[MAXDIM],x[MAXDIM];
	double alpha = 1.5, one = 1.0;

	ndmx = MAXDIV;
	mds = 1;
	for ( i = 0; i < dimension; i++ ) { xl[i] = 0.; xu[i] = one; }
	ndo = 1;
	for ( j = 0; j < dimension; j++ ) xi[0][j] = one;

	it = 0;
	scalls = schi = swgt = si2 = si = 0;

	nd = ndmx;
	ng = 1;
	if ( mds != 0 ) {
	    ng = (int)pow( ncalls*0.5, 1.0/dimension );
    	mds = 1;
		if ( (2*ng-ndmx) >= 0 ) {
			mds = -1;
			npg = ng/ndmx+1;
			nd = ng/npg;
			ng = npg*nd;
		}
	}
	k = iipow(ng,dimension);
	npg = ncalls/k;
	if ( npg < 2 ) npg = 2;
	calls = npg*k;
	dxg = one/ng;
	dv2g = ipow(dxg,2*dimension)/(npg*npg*(npg-one));
	xnd = nd;
	ndm = nd-1;
	dxg *= xnd;
	xjac = one;
	for ( j = 0; j < dimension; j++ ) {
		dx[j] = xu[j]-xl[j];
		xjac *= dx[j];
	}
	if ( nd != ndo ) {
		rc = ndo/xnd;
		for ( j = 0; j < dimension; j++ ) {
			i = k = 0;
			dr = xn = xo = 0;
			while ( i < ndm ) {
				while ( rc > dr ) {
					dr += one;
					xo = xn;
					xn = xi[k++][j];
				}
				dr -= rc;
				xin[i++] = xn-(xn-xo)*dr;
			}
			for ( i = 0; i < ndm; i++ ) xi[i][j] = xin[i];
			xi[nd-1][j] = one;
		}
		ndo = nd;
	}
	if ( printflag ) {
		if ( printflag == 10 ) {
			printf("vegas: dimension = %d, ncalls = %8.0f, itmax = %d",dimension,calls,itmax);
			printf(", acc = %e, mds = %d, nd = %d\n",accuracy,mds,nd);
		}
		else {
			printf("input parameters for vegas:  dimension = %3d  ncalls = %8.0f\n",dimension,calls);
			printf("                             it = %5d  itmax =    %5d\n",it,itmax);
			printf("                             acc = %e\n",accuracy);
			printf("                             mds = %3d\n",mds);
			printf("                             nd = %4d\n",nd);
		}
	}
reiterate:
	it++;
	tsi = ti = 0.;
	for ( j = 0; j < dimension; j++ ) {
		kg[j] = 1;
		for ( i = 0; i < nd; i++ ) {
			nxi[i][j] = 0;
			d[i][j] = di[i][j] = ti;
		}
	}
reloop:
	f2b = fb = 0;
	for ( k = 0; k < npg; k++ ) {
		for ( j = 0; j < dimension; j++ ) qran[j] = ranf(0);
		wgt = xjac;
		for ( j = 0; j < dimension; j++ ) {
			xn = (kg[j]-qran[j])*dxg;
			ia[j] = (int)xn;
			iaj   = ia[j];
			iaj1  = iaj-1;
			if ( iaj > 0 ) {
				xo = xi[iaj][j]-xi[iaj1][j];
				rc = xi[iaj1][j]+(xn-iaj)*xo;
			}
			else {
				xo = xi[iaj][j];
				rc = (xn-iaj)*xo;
			}
			x[j] = xl[j]+rc*dx[j];
			wgt *= xo*xnd;
		}
		f    = fcn(x)*wgt;
		f2   = f*f;
		fb   += f;
		f2b  += f2;
		for ( j = 0; j < dimension; j++ ) {
			iaj  = ia[j];
			(nxi[iaj][j])++;
			di[iaj][j] += f/calls;
			if ( mds >= 0 ) d[iaj][j] += f2;
		}
	}
	f2b = f2b*npg;
	f2b -= fb*fb;
	ti  += fb;
	tsi += f2b;
	if ( mds < 0 ) {
		for ( j = 0; j < dimension; j++ ) d[ia[j]][j] += f2b;
	}
	for ( k = dimension-1; k >= 0; k-- ) {
		kg[k] = ( kg[k] % ng ) + 1;
		if ( kg[k] != 1 ) goto reloop;
	}
	ti /= calls;
    tsi *= dv2g;
	ti2 = ti*ti;
	wgt = ti2/tsi;
	si += ti*wgt;
	si2 += ti2;
	swgt += wgt;
	schi += ti2*wgt;
	scalls += calls;
	avgi = si/swgt;
	sd = swgt*it/si2;
	chi2a = 0;
	if ( it > 1 ) chi2a = sd*(schi/swgt-avgi*avgi)/(it-1);
	sd = sqrt(one/sd);
	if ( printflag ) {
		tsi = sqrt(tsi);
		if ( printflag == 10 ) {
			printf("%d%20.8g%12.4g%20.8g%12.4g%12.4g\n",it,ti,tsi,avgi,sd,chi2a);
		}
		else {
			printf("\nintegration by vegas.\n\n");
			printf("iteration no %3d.      integral =%16.10g\n",it,ti);
			printf("                       std dev  =%16.10g\n",tsi);
			printf("accumulated results.   integral =%16.10g\n",avgi);
			printf("                       std dev  =%16.10g\n",sd);
			printf("                 chi^2 per itn  =  %10.4g\n",chi2a);
		}
		if ( printflag < 0 ) {
			for ( j = 0; j < dimension; j++ ) {
				printf("data for axis%2d\n",j+1);
				printf("       x         delt i     convce        ");
/*				printf("       x         delt i     convce        "); */
				printf("       x         delt i     convce\n");
				for ( i = 0; i < nd; i++ ) {
					printf("%12.4g%12.4g%12.4g",xi[i][j],di[i][j],d[i][j]);
/*					if ( (i%3) == 2 ) printf("\n"); */
					if ( (i%2) == 1 ) printf("\n");
					else              printf("     ");
				}
				if ( (nd%3) != 0 ) printf("\n");
			}
		}
	}
	for ( j = 0; j < dimension; j++ ) {
		dt[j] = 0;
		for ( i = 0; i < nd; i++ ) {
			if ( nxi[i][j] > 0 ) d[i][j] /= nxi[i][j];
			dt[j] += d[i][j];
		}
	}
	for ( j = 0; j < dimension; j++ ) {
		rc = 0;
		for ( i = 0; i < nd; i++ ) {
			r[i] = 0;
			if ( d[i][j] > 0 ) {
				xo = dt[j]/d[i][j];
				r[i] = pow((xo-one)/(xo*log(xo)),alpha);
			}
			rc += r[i];
		}
		rc /= xnd;
		k = 0;
		dr = r[k];
		xo = 0;
		xn = xi[k][j];
		for ( i = 0; i < ndm; ) {
			while ( rc > dr ) {
				k++;
				dr += r[k];
				xo = xn;
				xn = xi[k][j];
			}
			dr -= rc;
			xin[i++] = xn-(xn-xo)*dr/r[k];
		}
		for ( i = 0; i < ndm; i++ ) xi[i][j] = xin[i];
		xi[nd-1][j] = one;
	}
	if ( ( it < itmax ) && ( fabs(accuracy) < fabs(sd/avgi) ) ) goto reiterate;
}

double rvec[100], seeds[24], carry = 0., twopm24;
static int mcall = 0;
static int ncall = 0;
void rcargo(long dummy);
void rcarry(int mcall);
static int notyet = 1;
static int i24 = 24, j24 = 10, iseeds[24];
static long icons = 2147483563;
static long ITWO24 = 0xFFFFFF;
static double TWOP12 = 4096.;

double
ranf(long dummy)
{
	if ( ncall == 0 ) {
		ncall = 1;
		if ( dummy ) rcargo(dummy);
	}
	if ( mcall == 0 ) {
		mcall = 100;
		rcarry(mcall);
	}
	return(rvec[--mcall]);
}

void rcargo(long inseed)
{
	long jseed = inseed, k;
	double twom24 = 1.0;
	int i;
	notyet = 0;
	for ( i = 0; i < 24; i++ ) {
		twom24 /= 2;
		k = jseed/53668;
		jseed = 40014*(jseed-k*53668) - k*12211;
		if ( jseed < 0 ) jseed += icons;
		iseeds[i] = jseed % ITWO24;
	}
	for ( i = 0; i < 24; i++ ) seeds[i] = iseeds[i]*twom24;
	i24 = 24;
	j24 = 10;
	carry = 0;
	if ( seeds[24] < seeds[14] ) carry = twom24;
	twopm24 = twom24;
}

void rcarry(int lenv)
{
	int ivec;
	double uni;
/*
	if not initialized -> default initialization with fixed number
*/
	if ( notyet ) rcargo((long)314159265);
	for ( ivec = 0; ivec < lenv; ivec++ ) {
		uni = seeds[--j24] - seeds[--i24] - carry;
		if ( uni < 0 ) { uni += 1.0; carry = twopm24; }
		else carry = 0.;
		seeds[i24] = uni;
		if ( i24 == 0 ) i24 = 24;
		if ( j24 == 0 ) j24 = 24;
		rvec[ivec] = uni;
	}
}

void rcarut(long *isdext)
{
	int i, icarry;
	for ( i = 0; i < 24; i++ ) isdext[i] = (int)(seeds[i]*TWOP12*TWOP12);
	icarry = ( carry > 0 );
	isdext[24] = 1000*j24 + 10*i24 + icarry;
}

void rcarin(long *isdext)
{
	int i;
	long isd;
	double twom24 = 1.0;
	for ( i = 0; i < 24; i++ ) twom24 /= 2;
	for ( i = 0; i < 24; i++ ) seeds[i] = isdext[i]*twom24;
	isd = isdext[24];
	carry = (isd % 10) * twom24;
	isd /= 10;
	i24 = isd % 100;
	j24 = isd / 100;
	twopm24 = twom24;
}

double ipow(double x,int y)
{
	int den;
	double z, u;
/*
	determines x to the power y
*/
	if ( y == 0 ) return(1.0);
	if ( y < 0 ) { den = 1; y = -y; }
	else den = 0;
	if ( y == 2 ) u = x*x;
	else {
		if ( ( y & 1 ) != 0 ) u = x;
		else u = 1.0;
		z = x;
		y >>= 1;
		while ( y ) {
			z = z*z;
			if ( ( y & 1 ) != 0 ) u *= z;
			y >>= 1;
		}
	}
	if ( den ) return(1.0/u);
	else return(u);
}

int iipow(int x,int y)
{
	int z, u;
/*
	determines x to the power y
*/
	if ( y == 0 ) return(1);
	if ( y < 0 ) {
		if ( x == 1 ) return(1);
		else if ( x == -1 ) {
			if ( ( y&1 ) != 0 ) return(-1);
			else return(1);
		}
		else if ( x == 0 ) {
			printf("Division by zero in iipow: %d^%d\n",x,y);
			exit(-1);
		}
		else return(0);
	}
	if ( y == 2 ) u = x*x;
	else {
		if ( ( y & 1 ) != 0 ) u = x;
		else u = 1;
		z = x;
		y >>= 1;
		while ( y ) {
			z = z*z;
			if ( ( y & 1 ) != 0 ) u *= z;
			y >>= 1;
		}
	}
	return(u);
}