import math,random,time
from cmath import *
import decimal as dc


g_gamma = 7
p_gamma = [0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]
eulers_e = 2.71828182845904523

def std_dev(list_temp):
	try:
		mean = sum(list_temp)/len(list_temp)
		
		sd = 0
		
		for v in list_temp:
			sd += (v**2 - mean**2)
		
		sd = sd/len(list_temp)
		return math.sqrt(sd)
	except:
		return -10
		
def gamma(z):
	z = complex(z)
	
	if z.real < 0.5:
		return pi / (sin(pi*z)*gamma(1-z))
	else:
		z -= 1
		x = p_gamma[0]
		for i in range(1, g_gamma+2):
			x += p_gamma[i]/(z+i)
		t = z + g_gamma + 0.5
		
		return sqrt(2*pi) * t**(z.real+0.5) * exp(-t.real) * x
		
def incomplete_gamma(s,x,alternate=False):
	#Described in computing the incomplete Gamma function to arbitrary precision  - Serge Winitzki
	summer = 0
	
	
	if alternate:
		fast = False
		if s > 0:
			for n in range(0,150):
				try:
					if fast:
						#using log - not as accurate
						num = log((x**(n+s))).real
						den = factorialRamanujan(n,logged=True)
						bit = (1.0/(s+n))*((-1)**n)*(e**(num - den))
					else:
						num = ((-1)**n)*(x**(n+s)).real
						den = ((s + n)*factorial(n))
						bit = num/den
		
		
					
					
					try:
						summer += bit
					except:
						print "Error", s,x,num,den
				except:
					pass
					
		return (gamma(s) - summer).real
	else:
		i = factorial(s-1)
		j = e**(-x)
		l = 0
		for k in range(0,s):
			l +=  float((x**k).real)/factorial(k)
		
		return  (i*j*l).real
		
def factorialGamma(m,s=0):
	if m == 0:
		return 1
	else:
		return gamma(m+1).real
	
def factorialRamanujan(m,s=0,logged=True):
	if m == 0:
		if logged:
			return log(1)
		else:
			return 1
		
	else:
		if logged:
			return (m*log(m)) - m + (log(m*(1+4*m*(1+2*m)))/6) + (log(pi)/2).real
		else:
			return e**((m*log(m)) - m + (log(m*(1+4*m*(1+2*m)))/6) + (log(pi)/2)).real
	
def factorial(m,s=0):
	value = long(1)
	if m != 0:
		while m != s:
			value = value*m
			m = m - 1
	return value 
	
def fishers(a,b,c,d):
	x = (factorialRamanujan(a+b) + factorialRamanujan(c+d) + factorialRamanujan(a+c) + factorialRamanujan(b+d)) 
	y = (factorialRamanujan(a) + factorialRamanujan(b) + factorialRamanujan(c) + factorialRamanujan(d) + factorialRamanujan(a+b+c+d))
	
	return e**(x - y).real
	
def poisson(expected,observed):
	expected = long(expected)
	
	try: 
		x = math.exp(-expected) 
		y =  pow(expected,observed)
		
		try:
			z =  math.log(factorial(observed))
		except:
			z =  factorialRamanujan(observed,logged=True).real
			
		if observed == 0:
			lnExpected = 1
		else:
			lnExpected = math.log(expected)

		return math.exp(observed*lnExpected - expected - z)
		
			
	except Exception,e:
		return 0.0
	
		
def cum_poisson(expected,observed):
	prob = 1.0
	for x in range(0,observed  +1):
		prob -= poisson(expected,x)
		
	
	if prob >= 0: 
		return prob
	else: 
		return 0.0
    
def binomial(k,n,p):
	if n > 140:
		return cum_poisson(int(n*p),k)
	else:
		c_bit = factorialRamanujan(n)/(factorialRamanujan(k)*factorialRamanujan(n-k))
		p_bit = (p**k)*((1-p)**(n-k))
		return c_bit * p_bit

def cum_binomial(k,n,p):
	summer = 0
	for i in range(k,n+1):
		summer += binomial(i,n,p)
	return summer

def uniform_product_density(p,n):
	a = (-1)**(n-1)
	b = factorialRamanujan(n-1)
	c = log(p).real**(n-1)

	return (a/b)*c
	
def cum_uniform_product(p,n,alternate = False):
	try:
		if p == 1:
			return 1
		else:
			if alternate:
				
				a = (-1)**(n)
				b = incomplete_gamma(n,-log(p))
				c = (log(p).real)**n
				d = factorialRamanujan(n-1,logged=False)
				e = (-log(p).real)**n
				#print p,"\t",n,"\t",a,"\t","%1.3g"%b,"\t","%1.3g"%c,"\t","%1.3g"%d,"\t","%1.5g"%e,"\t",
				return (a*b*c)/(d*e)
				
			else:
				a = incomplete_gamma(n,-log(p))
				b = factorial(n-1)
				#print p,"\t",n,"\t",a,"\t","%1.3g"%b,"\t",
				return a/b
	except:
		print "Error in Sig correction:",p,n
		raise
		
def erf(z):
	summer = 0
	for i in range(0,100):
	
		try:
			num = ((-1.0)**i) * z**(2*i+1)
		except OverflowError:
			print "Overflow"

		denum = factorial(i)*(2*i+1)
		summer += num/denum
		
	
	erf = summer*(2 /math.sqrt(math.pi)) 
	if erf > 1 or erf < -1:
		erf = 1
		
	return erf
		
def choose(k,n):
	return factorial(n)/(factorial(k)*factorial(n-k))

def product(values):
	#prod = 1
	prod = dc.Decimal(str(1))
	for v in values:
		prod *= dc.Decimal(str(v))
		
	return prod
	
	
def rlcProb(rlc,u,sd):
	if rlc > 6.2:
		return 0.00000000001
	if rlc < -6.2:
		return 1
	else:
		return 1 - (1 +erf((rlc - u)/(abs(sd)*math.sqrt(2))))/2
		
if __name__ == "__main__":
	st = time.time()	
	"""
	# test cum_uniform_product
	summer = 1
	
	
	binDict = {}
	samples = 100000
	for i in range(0,samples):
		randInt = int(log(random.random()*random.random()*random.random()*random.random()*random.random(),10).real)
		if randInt in binDict:
			binDict[randInt] += 1
		else:
			binDict[randInt] = 1
		
	sorter = binDict.keys()
	sorter.sort()
	#sorter.reverse()
	summer = 0
	
	for bin in sorter:
		summer += binDict[bin]
		print bin,binDict[bin],10**float(bin),float(samples - summer)/samples,1 - cum_uniform_product(10**float(bin),5), cum_uniform_product(10**float(bin),5)
	#"""
	
	#for i in range(0,10000):
	#	print i,factorial(i)
	
	print poisson(2,1.43)
	print cum_poisson(2,4)
	print "+++++"
	#print cum_poisson(231,229)
	"""
	for n in range(1,10):
		for  i in range(1,15):
			p = 0.1**i
			print n,"\t",
			print i,"\t",
			
			print poisson(n,i)
			#print "%1.5g"%(cum_uniform_product(p,n)),"\t",
			#print cum_uniform_product(p,n,alternate=True)
	
			#print p,-log(p),"%1.3g"%incomplete_gamma(n,-log(p)),"%1.3g"%incomplete_gamma(n,-log(p),True)
			
		print time.time() - st
		
		"""