#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 fun(real *x){return 1/(1-cos(x[0])*cos(x[1])*cos(x[2]))/pi/pi/pi;}
//real fun(real *x){return 4*pow(x[0],3);}

void main(){
int n=3,N=1000; real *a,*b,sum,sum2,res,err,vol=1;
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);

printf("--------------\nstratified, N=%d\n",N);
strata(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);
}