#
# jointd.py
#
# Calculates the Cummulants of a Joint Gaussian
# using a random uniform distribution generator
#
# Carlos Rocha - February 2005
##############################

from math import *
from time import *

# THE Random Number Generator
##############################

class RandDouble:
    def __init__(self):
       self.seed =  int(time()) # Initial Random Seed
    
    def next(self):
        """Returns a random number from a uniform distribution [0,1]"""
        self.seed = (8121 * self.seed + 28411) % 134456
        return self.seed/134456.0
        

##############################

def main() :
    repeat = 100000  # Number of samples

    rand = RandDouble()

    sumM1 = 0 # The Sum of the raw moments Mx
    sumM2 = 0
    sumM3 = 0

    count = 0

    for i in range(repeat):

        x1 = rand.next()
        x2 = rand.next()

        if x1 != 0.0 and x2 != 0.0:

            y1 = sqrt(-2*log(x1))*sin(x2)
            y2 = sqrt(-2*log(x1))*cos(x2)

            p = exp(-(pow(y1,2)+pow(y2,2))/2)

            sumM1 += p
            sumM2 += pow(p,2)
            sumM3 += pow(p,3)
            count += 1

    M1 = sumM1 / count
    M2 = sumM2 / count
    M3 = sumM3 / count

    C1 = M1
    C2 = M2 - pow(M1,2)
    C3 = M3 - 3 * M1 * M2 + 2 * pow(M1,3)

    print "C1 = %f" % C1
    print "C2 = %f" % C2
    print "C3 = %f" % C3

if __name__ == "__main__" : main()