#include <stdio.h>
#include <string.h>
#include <math.h>
#include "rng.h"

#define NSAMPLES    10000


float sc1=0.0;
float sc2=0.0;
float sc3=0.0;

unsigned int		SEED[] = {
	0xffffffffU,	0xffffffffU, 0xffffffffU, 0xffffffffU,
	0xffffffffU,	0xffffffffU, 0xffffffffU, 0xffffffffU,
	0xffffffffU,	0xffffffffU, 0xffffffffU, 0xffffffffU,
	0xffffffffU,	0xffffffffU, 0xffffffffU, 0xffffffffU};

float fy1(float x1, float x2) {
    return sqrt(-2*log(x1))*sin(x2);
}

float fy2(float x1, float x2) {
    return sqrt(-2*log(x1))*cos(x2);
}

int main(int argc, char **argv) {
    float sc1=0.0,sc2=0.0,sc3=0.0;
    float x1,x2,y1,y2,pr;
    float e1,e2,e3;
    int i;

    RNG_init(SEED);
    for (i=0;i<NSAMPLES;i++) {
        x1 = RNG_rand();
        x2 = RNG_rand();

        if ((x1 == 0.0) || (x2 == 0.0)) {
            i--;
            continue;
        }

        y1 = fy1(x1,x2);
        y2 = fy2(x1,x2);

        pr = exp(-((y1*y1)+(y2*y2))/2.0);

        sc1 += pr;
        sc2 += pr*pr;
        sc3 += pr*pr*pr;
    }

    e1 = sc1/NSAMPLES;
    e2 = sc2/NSAMPLES;
    e3 = sc3/NSAMPLES;

    printf("C1=%f\nC2=%f\nC3=%f\n",e1,e2-e1*e1,e3-3*e1*e2+2*e1*e1*e1);

    return 0;
}