from numpy import *
from pylab import *

# number of bits in shift register
M = 24
# initialize register with all zeros
bits = [0]*M
# populate the last element of the register
bits[-1] = 1

# load lags data from file
lags = loadtxt("lag_values.txt", dtype = {'names': ('M', 'i'), 'formats': ('i2', 'S9')})

# extract relevant row
trow = lags[M-2]
# extract string of tap indices
data = trow[1]
# initialize taps with all zeros
taps = [0]*M

# set taps
start = 0
while True:
	# find the next comma
	end = data.find(",", start)
	if end == -1:
		# extract the last tap index
		taps[int16(data[start:len(data)]) - 1] = 1
		break
	else:
		# extract tap index
		taps[int16(data[start:end]) - 1] = 1
		# set the new start search value
		start = end + 1

# initialize powers of 2
def f(x): return 2**x
twos = map(f, range(M))

# define transformations
def y1(x1, x2): return sqrt(-2*log(x1))*sin(x2)
def y2(x1, x2): return sqrt(-2*log(x1))*cos(x2)		

# initialize output
all_x1 = []
all_x2 = []
all_y1 = []
all_y2 = []

# number of randoms generated
N = 2000
# to reduce correlation of LFSR, leap forward M iterations
leap = M

for i in range(N):
	for j in range(leap):
		# increment the register
		bits.insert(0, mod(dot(bits, taps), 2))
		# push all bits to the right
		del bits[-1]
	# convert from binary to decimal, uniform distribution from 0 to 1
	x1 = float(dot(bits, twos))/(2**M - 1)
	# save x1
	all_x1.append(x1)
	
	for j in range(leap):
		# increment the register
		bits.insert(0, mod(dot(bits, taps), 2))
		# push all bits to the right
		del bits[-1]
	# convert from binary to decimal, uniform distribution from 0 to 2*pi
	x2 = float(dot(bits, twos))*2*pi/(2**M - 1)
	# save x2
	all_x2.append(x2)
	
	# compute and save y1, y2
	all_y1.append(y1(x1, x2))
	all_y2.append(y2(x1, x2))

# compute first, second, and third cumulants
C1_1 = float(sum(all_y1))/N
print 'The first cumulant of y1 is ' + str(C1_1)
C1_2 = float(sum(all_y2))/N
print 'The first cumulant of y2 is ' + str(C1_2)
C2_1 = float(sum([x**2 - C1_1**2 for x in all_y1]))/N
print 'The second cumulant of y1 is ' + str(C2_1)
C2_2 = float(sum([x**2 - C1_1**2 for x in all_y2]))/N
print 'The second cumulant of y2 is ' + str(C2_2)
C3_1 = float(sum([x**3 - 3*C2_1*C1_1 + 2*C1_1**3 for x in all_y1]))/N
print 'The third cumulant of y1 is ' + str(C3_1)
C3_2 = float(sum([x**3 - 3*C2_1*C1_1 + 2*C1_1**3 for x in all_y2]))/N
print 'The third cumulant of y2 is ' + str(C3_2)

# plot and save distribution
fig = figure()
ax = fig.add_subplot(111)
ax.scatter(all_y1, all_y2)
ax.set_aspect('equal')
title('LFSR with $M = 24$ mapped to Gaussian')
xlim([-3, 3])
ylim([-3, 3])
savefig('Gaussian_LFSR_M24_ff24.pdf')
show()