def bezier2polynomial(p, numpy_ordering=True, return_poly1d=False):
    """Converts a tuple of Bezier control points to a tuple of coefficients
    of the expanded polynomial.
    return_poly1d : returns a numpy.poly1d object.  This makes computations
    of derivatives/anti-derivatives and many other operations quite quick.
    numpy_ordering : By default (to accommodate numpy) the coefficients will
    be output in reverse standard order."""
    if len(p) == 4:
        coeffs = (-p[0] + 3*(p[1] - p[2]) + p[3],
                  3*(p[0] - 2*p[1] + p[2]),
                  3*(p[1]-p[0]),
                  p[0])
    elif len(p) == 3:
        coeffs = (p[0] - 2*p[1] + p[2],
                  2*(p[1] - p[0]),
                  p[0])
    elif len(p) == 2:
        coeffs = (p[1]-p[0],
                  p[0])
    elif len(p) == 1:
        coeffs = p
    else:
        # https://en.wikipedia.org/wiki/Bezier_curve#Polynomial_form
        n = len(p) - 1
        coeffs = [fac(n)//fac(n-j) * sum(
            (-1)**(i+j) * p[i] / (fac(i) * fac(j-i)) for i in range(j+1))
            for j in range(n+1)]
        coeffs.reverse()
    if not numpy_ordering:
        coeffs = coeffs[::-1]  # can't use .reverse() as might be tuple
    if return_poly1d:
        return poly1d(coeffs)
    return coeffs

def predict_errors(n,ct_main,freq_main,kmer_size,error_rate):
    '''
    predict kmer freq for kmers with exactly n errors
    '''
    
    #how many possible error kmers are there
    #k-choose-n combinations of bases can be mutated
    #for each combination there are 3^n combinations of bases to mutate to
    #pos = 3.0**n * sp.misc.comb(kmer_size,n)
    pos = 3.0**n * fac(kmer_size) / (fac(n)*fac(kmer_size-n))
    
    #probability of any one particular error kmer being generated
    #when the kmer is read
    #n bases misread as the particular error base
    #remainder of bases are not misread
    prob = (error_rate/3.0)**n * (1.0 - error_rate)**(kmer_size-n)

    freq = [0.0] * len(ct_main)

    #for every count in the error-free histogram
    for i,ct1 in enumerate(ct_main):
        mean = ct1 * prob
        val1 = freq_main[i] * pos * math.exp(-mean)
        
        #generate error predictions for every count in the error histogram
        for j,ct2 in enumerate(ct_main):
            #poisson for count ct2
            #scaled by val1 'trials'
            #freq[j] += val1 * poisson(mean,ct2)
            
            freq[j] += val1 * mean**ct2 / fac(ct2)
            if ct2 == 170: break #fac(171) is too big to be a float
    
    return freq

def solve_problem():
    counter = 0
    for n in range(1, 101):
        for r in range(0, n+1):
            if fac(n)/fac(r)/fac(n-r) > 1000000:
                counter += 1
    return counter

def part_p(n,k):
	result = 0
	for y in range(k, n+1):
		len_combs = fac(n)/(fac(n-y)*fac(y))
		result+=len_combs

	print result%1000000

 def a_function(l, r, i, l_1, l_2, pa, pb, pc, g):
     e = 1 / (4*g)
     out1 = Binomial.coefficient(l, l_1, l_2, pa, pb)
     out2 = (-1)**i * fac(l) * pc**(l - 2*r - 2*i) * e**(r + i)
     out3 = fac(r) * fac(i) * fac(l - 2*r - 2*i)
     ans = (-1)**l * out1 * (out2/out3)
     return ans

def _significance_direct_on_off(n_on, n_off, alpha):
    """Compute significance directly via Poisson probability.

    Use this method for small n_on < 10.
    In this case the Li & Ma formula isn't correct any more.

    * TODO: add reference
    * TODO: add large unit test coverage (where is it numerically precise enough)?
    * TODO: check coverage with MC simulation
    * TODO: implement in Cython and vectorize n_on (accept numpy  array n_on as input)
    """
    from math import factorial as fac
    from scipy.stats import norm

    # Compute tail probability to see n_on or more counts
    probability = 1
    for n in range(0, n_on):
        term_1 = alpha ** n / (1 + alpha) ** (n_off + n + 1)
        term_2 = fac(n_off + n) / (fac(n) * fac(n_off))
        probability -= term_1 * term_2

    # Convert probability to a significance
    significance = norm.isf(probability)

    return significance

def mend_prob(gen, num):
	prog = 2**gen  						# number of progeny in that generation
	result = 0
	for x in range(num, prog+1):		#sum the cumulative probability of exactly num or greater individuals occuring
		result += ((fac(prog)/(fac(x)*fac(prog-x)))*.25**x*.75**(prog-x))
										#formula for binomial probability 
	return ("%.3f" % result)			#round to 3 decimal places

def Hp_norm(p,k,l):
    s = 0

    for i in range(p+1):
        for j in range(i+1):
            s = s + (k*pi)**(2*(i-j))*(l*pi)**(2*j)*fac(i)/(fac(j)*fac(i-j)) 
        
    return sqrt(0.25*s)

def free_path(start, end):
    """Given start and end coords, calculate the total number of paths, assuming
    all entries are unblocked --> (m + n)! / (m! * n!)
    where m = num rows, n = num cols
    and start and end are tuples of the form (row pos, col pos)"""
    # num rows
    m = end[0] - start[0]
    # num cols
    n = end[1] - start[1]
    return fac(m + n) / (fac(m) * fac(n))

 def pn(self,n):
     """
     Probability in "n" state, holding n > 0 and n <= k
     """
     assert(n > 0 and n <= self.k)
     
     if n < self.c:
         return (pow(self.l, n) / (fac(n) * pow(self.m, n))) * self.p0()
     elif self.c <= n:
         return pow(self.l, n) / (pow(self.c, (n - self.c)) * fac(self.c) * pow(self.m, n)) * self.p0()

def per_match(seq):
	au_count = 0
	gc_count = 0
	for nuc in seq:
		if nuc == 'A':
			au_count += 1
		if nuc == 'G':
			gc_count += 1
	result = fac(au_count)*fac(gc_count)
	print result

def successorTerm(t):
    C = lambda n, k: round(fac(n) / (fac(k) * fac(n - k)))
    if t.exponent == 0:
        return (t,)
    if t.exponent == 1:
        return (t, Term(t.coefficient, 0))
    elif t.exponent > 1:
        return map(
            lambda z: Term(t.coefficient * C(z[0], z[1]), z[1]),
            zip(repeat(t.exponent, t.exponent + 1), range(t.exponent + 1)),  # n  # k
        )

def cool_sentence(n, a, v):
    from math import factorial as fac

    sent_count = 0

    if a < 8:
        for i in range(1, a + 1):
            sent_count += int(fac(a) / fac(a - i))
    else:
        for i in range(1, 8):
            sent_count += int(fac(a) / fac(a - i))

    return sent_count * n * v

 def b_function(self, l, ll, r, rr, i, l_1, l_2, a_x, b_x, p_x, g_1, l_3, l_4, c_x, d_x, q_x, g_2):
     sigma = self.sigma
     pa_x = p_x - a_x
     pb_x = p_x - b_x
     qc_x = q_x - c_x
     qd_x = q_x - d_x
     c_x = p_x - q_x
     delta = (1/(4*g_1)) + (1/(4*g_2))
     out1 = (-1)**l * sigma(l, l_1, l_2, pa_x, pb_x, r, g_1) * sigma(ll, l_3, l_4, qc_x, qd_x, rr, g_2)
     out2 = (-1)**i * (2 * delta)**(2 * (r + rr)) * fac(l + ll - 2*r - 2*rr) * delta**i * c_x**(l + ll - 2*(r + rr + i))
     out3 = (4 * delta)**(l + ll) * fac(i) * fac(l + ll - 2*(r + rr + i))
     ans = out1 * (out2 / out3)
     return ans

def lotto(use_random_coupon):
    weeks_to_play = 10000000
    max_val = 36
    all_numbers = [x for x in range(1, max_val + 1)]
    coupon = []

    if use_random_coupon:
        print("Using random generated coupon")
        coupon_rows = 10
        row_len = 7
        while len(coupon) < coupon_rows:
            row = sample(all_numbers, row_len)
            row.sort()
            coupon.append(row)
    else:
        print("using manually generated coupon")
        coupon = [
            [5, 10, 12, 19, 24, 26, 30],
            [7, 11, 12, 23, 25, 32, 34],
            [7, 11, 13, 24, 25, 29, 31],
            [4, 10, 14, 15, 21, 25, 29],
            [5,  6, 10, 11, 17, 30, 32],
            [4,  5,  8,  9, 10, 16, 18],
            [6, 12, 14, 15, 26, 33, 36],
            [1,  7, 10, 12, 16, 19, 34],
            [2, 10, 21, 22, 27, 30, 34],
            [6,  9, 17, 22, 25, 27, 34]]
        coupon_rows = len(coupon)
        row_len = len(coupon[0])

    wins = 0
    for week in range(1, weeks_to_play + 1):
        win_row = sample(all_numbers, row_len)
        win_row.sort()
        for row in range(0, coupon_rows):
            matches = 0
            for number in win_row:
                if number in coupon[row]:
                    matches += 1
                if matches == row_len:
                    wins += 1
                    print('Jackpot! Row {} in week {} with {}'
                          .format(row + 1, week, coupon[row]))
    else:
        print('Game over! You won {} times in {} years!'
              .format(wins, int(weeks_to_play * 7 / 365.25)))
        print('Your win ratio is 1 in {}'
              .format(int(weeks_to_play * coupon_rows / wins)))
        win_prob = int(fac(max_val) / (fac(row_len) * fac(max_val - row_len)))
        print('The mathematical win ratio is 1 in {}'
              .format(win_prob))

def prob(k,N):
    from fractions import Fraction as F
    from math import factorial as fac
    nCr = lambda n,r: int(fac(n)/(fac(r)*fac(n-r)))
    progeny_k = 2**k
    i = N
    prob_out = 0
    while i <= progeny_k:
        prob_ = (F(1/4)**i * F(3,4)**(progeny_k-i)) * nCr(progeny_k,i)
        dummy = (F(1/4)**i * F(3,4)**(progeny_k-i))
        prob_out += prob_
        print(i,dummy,nCr(progeny_k,i),prob_,sep='\t')
        i+=1
    return float(prob_out)

 def p0(self):
     """
     Probability in "0" state
     """
     if not self.p0_optimizer:
         self.p0_optimizer = True
         acum = 0
         for n in range(0, self.c):
             acum += pow(self.r, n) / fac(n)
         aux = ((pow(self.r, self.c))/(fac(self.c))) * (((1 - pow(self.rho, self.k - self.c + 1)) / (1 - self.rho)))
         self.p0_aux = 1 / (acum + aux)
         return self.p0_aux
     else:
         return  self.p0_aux

 def uniquePaths(self, m, n):
     """
     :type m: int
     :type n: int
     :rtype: int
     """
     if m>n:
         tmp=m;m=n;n=tmp;
     if m==0:
         return 0
     a=fac(n+m-2)
     b=fac(m-1)
     c=fac(n-1)
     return a/b/c

def anagrams(s, v=10000):
    anagram_list = []
    while len(anagram_list) < v and len(anagram_list) < fac(len(s)):
        candidate = ''.join(sample(s, len(s)))
        if candidate not in anagram_list:
            anagram_list.append(candidate)
    return anagram_list

 def lq(self):
     """
     Mean length og the tail
     """
     a = (self.p0() * pow(self.r, self.c) * self.rho) / (fac(self.c) * pow(1-self.rho, 2))
     b = 1 - pow(self.rho, self.k - self.c + 1) - ((1-self.rho) * (self.k - self.c + 1) * pow(self.rho, self.k - self.c))
     return a * b