/* calculation of  Inverse Transform Method */
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;

unsigned int rnd()
{
    static unsigned long long int s1 = 1234UL;
    static unsigned long long int s2 = 5678UL;

    s1 = s1 ^ (s1 >> 17);
    s1 = s1 ^ (s1 << 31);
    s1 = s1 ^ (s1 >>  8);
    s2 = (s2 & 0xffffffff)*4294957665UL + (s2>>32);

    return (unsigned int)(s1 ^ s2);
}

double drnd() { return (double)rnd()/4294967296.; }

double func(double a,double b,double c,double d,double e,double f)
{
    return sin(a)+sin(b*2.)+sin(c*3.)+cos(d)+cos(e*2.)+cos(f*3.);
}

double intfunc(double a,double b,double c,double d,double e,double f)
{
    return -cos(a)-cos(b*2.)/2.-cos(c*3.)/3.+sin(d)+sin(e*2.)/2.+sin(f*3.)/3.;
}

int main()
{
    double intf_exact = intfunc(1.,1.,1.,1.,1.,1.)-intfunc(0.,0.,0.,0.,0.,0.);
	cout << "Exact  = "  << intf_exact << endl;

//    printf("Exact  = %.16f\n",intf_exact);

    double intf_rand = 0.;
    for(int i=0;i<1000000;i++) {
        intf_rand += func(drnd(),drnd(),drnd(),drnd(),drnd(),drnd());
    }
    intf_rand /= 1000000;

    cout << "Random = "  << intf_rand << endl;

//    printf("Random = %.16f\n",intf_rand);

	double intf_boxes = 0.;
	double h = 0.05;
	double a,b,c,d,e,f;

	a = h/2.; while(a<1.) {
	b = h/2.; while(b<1.) {
	c = h/2.; while(c<1.) {
	d = h/2.; while(d<1.) {
	e = h/2.; while(e<1.) {
	f = h/2.; while(f<1.) {

		intf_boxes += func(a,b,c,d,e,f);

	f += h;}
	e += h;}
	d += h;}
	c += h;}
	b += h;}
	a += h;}

	intf_boxes *= h*h*h*h*h*h;

	cout << "Boxes  = "  << intf_boxes << endl;
//	printf("Boxes  = %.16f\n",intf_boxes);
}