#include "real.h"
//#define pi 3.1415927
#define pi 1
#include <stdio.h>
#include <math.h>

void plainmc(int n,real a[],real b[],real (*func)(real*),
      int N,real *sum,real *sum2,void (*rnd)(int,real*,real*,real*));
void pseudo(int n,real *a,real *b, real *x);
void quasi(int n,real *a,real *b, real *x);
real strata(int n,real *a,real *b,real (*func)(real*),
      int M,real acc,void (*rnd)(int,real*,real*,real*));

int K;
//real fun(real *x){K++;return 1/(1-cos(x[0])*cos(x[1])*cos(x[2]))/pi/pi/pi;}
real fun(real *x){K++;return 4*pow(x[0],3)*4*pow(x[1],3)*4*pow(x[2],3);}
//real fun(real *x){K++;return 1;}

void main(){
int n=3,N=pow(20,n); real *a,*b,sum,sum2,res,err,vol=1,var;
double exact=1.3932039296856768591842462603255;

a=(real*)malloc(n*sizeof(real)); b=(real*)malloc(n*sizeof(real));
a[0]=0; a[1]=0; a[2]=0; b[0]=pi; b[1]=pi; b[2]=pi;

{int i; for(i=0;i<n;i++) vol*=(b[i]-a[i]);}

printf("--------------\npseudo-random, N=%d\n",N);
plainmc(n,a,b,&fun,N,&sum,&sum2,&pseudo);
res=sum/N*vol; err=sqrt(sum2/N/N-sum*sum/N/N/N)*vol;
printf("result\t=%g\n",res);
printf("exact\t=%g\n",exact);
printf("estimated error\t=%g\n",err);
printf("actual error\t=%g\n",res-exact);

printf("--------------\nquasi-random, N=%d\n",N);
plainmc(n,a,b,&fun,N,&sum,&sum2,&quasi);
res=sum/N*vol; err=sqrt(sum2/N/N-sum*sum/N/N/N)*vol;
printf("result\t=%g\n",res);
printf("exact\t=%g\n",exact);
printf("estimated error\t=%g\n",err);
printf("actual error\t=%g\n",res-exact);

N=100 ; err=0.01; K=0;
printf("--------------\nstratified, N0=%d, vol=%g\n",N,vol);
res=strata(n,a,b,&fun,N,err*err,&pseudo);
printf("K\t=%i\n",K);
printf("result\t=%g\n",res);
printf("exact\t=%g\n",exact);
printf("estimated error\t=%g\n",err);
printf("actual error\t=%g\n",res-exact);
}