import operator
import math
import time
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import numpy as np


class LFSR:
    def __init__(self, n, seed):
        self.tapsDict = {2 : (1,2),3: (1,3),4: (1,4), 5: (2,5),6: (1,6),
            7: (3,7), 8: (2,3,4,8),9: (4,9),10: (3,10),11: (2,11),12: (1,4,6,12),
            13: (1,3,4,13),14: (1,6,10,14),15: (1,15),16: (1,3,12,16),17: (3,17),
            18: (7,18),19: (1,2,5,19),20: (3,20),21: (2,21),22: (1,22),23: (5,23),
            24: (1,2,7,24),25: (3,25),26: (1,2,6,26),27: (1,2,5,27),28: (3,28),
            29: (2,29),30: (1,2,23,30),31: (3,31),32: (1,2,22,32),33: (13,33),
            34: (1,2,27,34)}
        self.n = n

        self.taps = self.tapsDict[self.n]
        self.shift_register = self.get_bits(seed, n)

    def shift(self):
        new_bit = sum([self.shift_register[tap-1] for tap in self.taps])%2
        self.shift_register = self.shift_register[:-1]
        self.shift_register.insert(0, new_bit)
        return new_bit

    def random(self, length):
        bitstring=""
        for i in range(length):
            bitstring += str(self.shift())
        return float(int(bitstring, 2))/pow(2,self.n)

    def get_bits(self, k, bits_to_get):
        """returns the first bits_to_get bits of k starting with the lowest"""
        bits = []
        for i in xrange(bits_to_get):
            bits.append(k%2)
            k = k >> 1
        return bits

if __name__== "__main__":
    num_samples = 10000
    my_LFSR = LFSR(16, int(1000 * time.time()))

    samples1 = [my_LFSR.random(16) for i in range(num_samples)]
    samples2 = [my_LFSR.random(16) for i in range(num_samples)]

    g_pairs = [[math.sqrt(-2*math.log(x1))*math.sin(2*math.pi*x2), math.sqrt(-2*math.log(x1))*math.cos(2*math.pi*x2)] for x1, x2 in zip(samples1, samples2)]
    g_samples = reduce(operator.__add__, g_pairs) 


    plt.hist(g_samples)
    plt.show()


    
    uniform_m1 = sum(samples1)/len(samples1)
    uniform_m2 = sum([x**2 for x in samples1])/len(samples1)
    uniform_m3 = sum([x**3 for x in samples1])/len(samples1)

    gauss_m1 = sum(g_samples)/len(g_samples)
    gauss_m2 = sum([x**2 for x in g_samples])/len(g_samples)
    gauss_m3 = sum([x**3 for x in g_samples])/len(g_samples)

    uniform_C1 = uniform_m1
    uniform_C2 = uniform_m2 - uniform_m1**2
    uniform_C3 = uniform_m3 - 3*uniform_m1*uniform_m2 + 2*uniform_m1**3

    gauss_C1 = gauss_m1
    gauss_C2 = gauss_m2 - gauss_m1**2
    gauss_C3 = gauss_m3 - 3*gauss_m1*gauss_m2 + 2*gauss_m1**3

    print "Uniform cumulants: "
    print "C1: ", uniform_C1
    print "C2: ", uniform_C2
    print "C3: ", uniform_C3
        
    print "Gaussian cumulants: "
    print "C1: ", gauss_C1
    print "C2: ", gauss_C2
    print "C3: ", gauss_C3