def f_f(j):
    n=2*j
    r=n//2
    num=math.factorial(n)
    den=math.factorial(n-r) * math.factorial(r)
    return num//den
    

def n_choose_k(n, k):
    """Computes n choose k."""

    if n < k:
        return 0
    else:
        return factorial(n) / factorial(k) / factorial(n - k)

def keyspace(n,m,k):

    lastterm = factorial(k) / (factorial(k - m))

    result = int(n * comb(n,m,exact=True) * lastterm)

    return result

def probabilityOfKeepSet(keep_set):
	probability = math.factorial(len(keep_set))
	for value in range(1,7):
		num_value_in_set = sum([1 for x in keep_set if x == value])
		probability /= math.factorial(num_value_in_set)
	probability = float(probability)/6**len(keep_set) # Divide by all possible choices
	return probability

def choose(a, b):
    if b == 0 or a == b:
        return 1
    else:
        numer = factorial(a)
        denom = factorial(b) * factorial(a-b)
        return numer//denom

 def solve(self, cipher):
     """
     Labeled permutation variants
     :param cipher: the cipher
     """
     N, M = cipher
     return math.factorial(N+M-1)/math.factorial(N)/math.factorial(M-1)%MOD

def multCoeff(x):
    a = 0
    b = 0
    c = 0
    d = 0
    e = 0
    f = 0
    g = 0
    h = 0
    i = 0
    for j in str(x):
        if j == '1':
            a += 1
        elif j == '2':
            b += 1
        elif j == '3':
            c += 1
        elif j == '4':
            d += 1
        elif j == '5':
            e += 1
        elif j == '6':
            f += 1
        elif j == '7':
            g += 1
        elif j == '8':
            h += 1
        elif j == '9':
            i += 1
    num = factorial(a+b+c+d+e+f+g+h+i)
    den = factorial(a)*factorial(b)*factorial(c)*factorial(d)*factorial(e)\
    *factorial(f)*factorial(g)*factorial(h)*factorial(i)
    return num/den

def test_factorial():
    """ Tests the   factorial   function. """
    # DONE: 3a. Implement this function, using it to test the NEXT
    #   function. Write the two functions in whichever order you prefer.
    #   Include at least 4 tests.
    #
    # ** Use the  math.factorial  function as an ORACLE for testing. **
    print()
    print('--------------------------------------------------')
    print('Testing the   factorial   function:')
    print('--------------------------------------------------')
    actual_answer = factorial(5)
    oracle_answer = math.factorial(5)
    test_case = 'factorial(5).  Actual, Oracle answers:'
    print('   Called ', test_case, actual_answer, oracle_answer)
    actual_answer = factorial(8)
    oracle_answer = math.factorial(8)
    test_case = 'factorial(8).  Actual, Oracle answers:'
    print('   Called ', test_case, actual_answer, oracle_answer)
    actual_answer = factorial(11)
    oracle_answer = math.factorial(11)
    test_case = 'factorial(11).  Actual, Oracle answers:'
    print('   Called ', test_case, actual_answer, oracle_answer)
    actual_answer = factorial(14)
    oracle_answer = math.factorial(14)
    test_case = 'factorial(14).  Actual, Oracle answers:'
    print('   Called ', test_case, actual_answer, oracle_answer)

def dig_fact():
    
    # import timeit and start timer to measure time of function
    import timeit
    start = timeit.default_timer()
    
    import math # import math package for factorial function
    
    sum_total = 0 # initialize sum variable to hold sum of digit factorials
    
    list = [] # initalize list to hold all numbers which are equal to the sum of the factorial of their digits
    
    for i in range(3,math.factorial(9)*7): # iterate through range and check conditions per problem
   
        sum_total = 0 # reset sum variable for each number in range above
   
        for j in range(0,len(str(i))): # iterate through each digit of number in range above
                
            sum_total += math.factorial(int(str(i)[j])) # add factorial of each digit of number to "sum" variable for checking problem condition below
            
        if i == sum_total: # append number to list if it is equal to the sum of the factorial of its digits
            
            list.append(i)
            
    print sum(list) # return solution        
            
    # stop timer and print total elapsed time
    stop = timeit.default_timer()
    print "Total time:", stop - start, "seconds"

 def check_topological_embedding_brute(self, verbose=False):
     '''
     Brute-force combinatoric method to check topological embedding
     '''
     # Clear table and nodemap:
     self.T = { superN:{ subN:None for subN in self.subD.nodes } 
               for superN in self.superD.nodes }  
     self.AB_nodemap = None
     N = len(self.subD.nodes) # number of subdesign nodes
     if len(self.superD.nodes) < N:
         # shortcut - fewer nodes in superdesign
         self.AB_nodemap = None
         return False
     # compute number of matchings:
     num_matchings = math.factorial(len(self.subD.nodes))*math.factorial(len(self.superD.nodes))
     count = 0
     for sub_perm in permutations(self.subD.nodes):
         for super_perm in permutations(self.superD.nodes, N):
             #sys.stdout.write("%d  \r" % count )
             #sys.stdout.flush()    # carriage returns don't work in Eclipse :-(
             if verbose:
                 print str(count) + " / " + str(num_matchings)
                 count += 1
             self.AB_nodemap = dict (zip(sub_perm, super_perm))
             if self.check_vertex2vertex():
                 if self.check_edge2path():
                     if self.check_vertex_disjointness():
                         return True
     # no embedding found.
     self.AB_nodemap = None
     return False

    def odf_sh(self):
        r""" Calculates the real analytical ODF in terms of Spherical
        Harmonics.
        """
        # Number of Spherical Harmonics involved in the estimation
        J = (self.radial_order + 1) * (self.radial_order + 2) // 2

        # Compute the Spherical Harmonics Coefficients
        c_sh = np.zeros(J)
        counter = 0

        for l in range(0, self.radial_order + 1, 2):
            for n in range(l, int((self.radial_order + l) / 2) + 1):
                for m in range(-l, l + 1):

                    j = int(l + m + (2 * np.array(range(0, l, 2)) + 1).sum())

                    Cnl = (
                        ((-1) ** (n - l / 2)) /
                        (2.0 * (4.0 * np.pi ** 2 * self.zeta) ** (3.0 / 2.0)) *
                        ((2.0 * (4.0 * np.pi ** 2 * self.zeta) ** (3.0 / 2.0) *
                          factorial(n - l)) /
                         (gamma(n + 3.0 / 2.0))) ** (1.0 / 2.0)
                    )
                    Gnl = (gamma(l / 2 + 3.0 / 2.0) * gamma(3.0 / 2.0 + n)) / \
                        (gamma(l + 3.0 / 2.0) * factorial(n - l)) * \
                        (1.0 / 2.0) ** (-l / 2 - 3.0 / 2.0)
                    Fnl = hyp2f1(-n + l, l / 2 + 3.0 / 2.0, l + 3.0 / 2.0, 2.0)

                    c_sh[j] += self._shore_coef[counter] * Cnl * Gnl * Fnl
                    counter += 1

        return c_sh

def basis(coords,b1,b2,nmax):

    hermites = np.loadtxt('hermite_coeffs.txt')
    xrot = coords[:,0]/b1
    yrot = coords[:,1]/b2
    npix = coords.shape[0]
    all_basis = np.zeros((npix,nmax+1,nmax+1))    

    gauss = np.exp(-0.5*(np.array(xrot)**2+np.array(yrot)**2))    
    n1 = 0                                  
    n2 = 0
    for n1 in range(0,nmax):
        n2 = 0
        while (n1+n2) <= nmax:
            norm = m.sqrt(m.pow(2,n1+n2)*m.pi*b1*b2*m.factorial(n1)*m.factorial(n2))
            k=0
            h1=0.
            h2=0.
            while (k <= n1):        
                h1 = h1+hermites[n1, k]*(np.array(xrot))**(n1-k)
                k=k+1
            k=0
            while (k <= n2):        
                h2 = h2+hermites[n2, k]*(np.array(yrot))**(n2-k)
                k=k+1
                
            all_basis[:, n1, n2] = gauss/norm*h1*h2;
            n2 = n2+1

    return all_basis

def n_choose_k(n: int, k: int) -> int:
    """ Return the binomial coefficient for n choose k.

    The binomial coefficient is defined as:

        n choose k = n! / (k!(n-k)!)

    For a number of objects n, which k choices allowed, this function reports
    the number of ways to make the selection if order is disregarded.

    We define n choose k == 0 for k > n.

    Args:
        n (int): A positive integer indicating the possible number of items to
            select.
        k (int): A positive integer indicating the number of items selected.

    Returns:
        int: The result of n choose k.

    Raises:
        ValueError: if n or k are not non-negative integers.
    """
    # Check that n and k are allowed
    if not countable.is_nonnegative_integer(n):
        raise ValueError("n is a not a non-negative integer")
    if not countable.is_nonnegative_integer(k):
        raise ValueError("k is a not a non-negative integer")
    # We define n choose k for k > n as 0, that is, there are no way
    # to choose more items than exist in a set.
    if k > n:
        return 0
    # Compute and return
    return math.factorial(n) // (math.factorial(k) * math.factorial(n - k))

def g_statistic(X, p=None, idx=None):
    """
    return g statistic and p value
    arguments: 
        X - the periodicity profile (e.g. DFT magnitudes, autocorrelation etc)
            X needs to contain only those period values being considered, 
            i.e. only periods in the range [llim, ulim]
    """
    # X should be real
    X = abs(numpy.array(X))
    if p is None: 
        power = X.max(0)
        idx = X.argmax(0)
    else:
        assert idx is not None
        power = X[idx]
    g_obs = power/X.sum()
    M = numpy.floor(1/g_obs)
    pmax = len(X)
    
    result = numpy.zeros((int(M+1),), float)
    pmax_fact = factorial(pmax)
    for index in xrange(1, min(pmax, int(M))+1):
        v = (-1)**(index-1)*pmax_fact/factorial(pmax-index)/factorial(index)
        v *= (1-index*g_obs)**(pmax-1)
        result[index] = v
    
    p_val = result.sum()
    return g_obs, p_val

    def get_score(self, subtract_cards):
        same_form = 0
        form_larger = 0

        if self.straight:
            for i in range(0, 12 - len(self.cards)):
                and_cards = True
                for j in range(len(self.cards)):
                    and_cards = and_cards and subtract_cards[(i + j)* 4 + self.cards[j].suit]

                if and_cards:
                    same_form += 1
                    if self.cards[0].rank > i:
                        form_larger += 1
        else:
            for i in range(0, 12):
                and_cards = True
                for j in range(len(self.cards)):
                    and_cards = and_cards and subtract_cards[i* 4 + self.cards[j].suit]

                if and_cards:
                    same_form += 1
                    if self.cards[0].rank > i:
                        form_larger += 1

            count_rank_2 = 0
            for j in range(12*4, len(subtract_cards)):
                if subtract_cards[j]:
                    count_rank_2 += 1

            if count_rank_2 >= len(self.cards):
                same_form += math.factorial(count_rank_2)/(math.factorial(count_rank_2 - len(self.cards)) * math.factorial(len(self.cards)))

        return same_form if same_form == 0 else form_larger * 1.0/same_form

def GetWaveFunction(n,l,m,x,y,z):
    r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
    theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))
    phi = lambda x,y,z: numpy.arctan2(y,x)
    #phi = lambda x,y,z: numpy.arctan(y/x)
    a0 = 1.
    R = lambda r,n,l: numpy.sqrt(((2.0/n/a0)**3)*(math.factorial(n-l-1))/(math.factorial(n+l))/2/n)*(2*r/n/a0)**l * numpy.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0)
    #R = lambda r,n,l: (2*r/n/a0)**l * numpy.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0)
    WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
    absWF = lambda r,theta,phi,n,l,m: abs((WF(r,theta,phi,n,l,m))**2)
    wf = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
    w = absWF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
    w[numpy.isnan(w)]=0
    #wf[numpy.isnan(wf)]=0

    phase = numpy.angle(wf)
    realwf = numpy.imag(wf)*10
    #realwf = numpy.fabs(numpy.real(wf)*10)# need to by 10 because of the pg.isosurface() function
    #phase = numpy.array(w)
    #xsize = w.size/w[0].size
    #ysize = w[0].size/w[0][0].size
    #zsize = w[0][0].size
    #for i in range(0,xsize):
    #    for j in range(0,ysize):
    #        for k in range(0,zsize):
    #            #angular = numpy.arctan2(wf[i][j][k].imag,wf[i][j][k].real)
    #            #if angular < 0:
    #            #    phase[i][j][k] = angular+2*numpy.pi
    #            #else:
    #            #    phase[i][j][k] = angular
    #            phase[i][j][k] = numpy.arctan2(wf[i][j][k].imag,wf[i][j][k].real)
    phase[numpy.isnan(phase)]=0
    realwf[numpy.isnan(realwf)]=0
    return w,phase, realwf

def fast_nth_perm(digits, n):

    #temporary storage variables
    curr_num = n - 1
    result = ""

    #go through the digits from the biggest down (i.e. 9 to 0)
    for i in reversed(range(len(digits))):
        #do integer division to get the quotient of curr_num / i!
        quotient = curr_num / math.factorial(int(i))
        #also get the remainder with mod
        remainder = curr_num % math.factorial(int(i))

        #the quotient is the index of the current digit in the desired
        #permutation. Add it to the result.
        result += (digits[quotient])
        # and remove it from the list of digits (each digit can only
        # appear once)
        digits = "%s%s" % (digits[:quotient], digits[quotient+1:])

        #update curr num to be the remainder of the previous
        #curr_num's division.
        curr_num = remainder
        
    return result

def tabela(n):
    somatorio = lambda termo: sum(termo(float(i)) for i in xrange(n))
    mysin = lambda x: somatorio(lambda i: (x ** (2 * i + 1)) * (-1) ** i / factorial(2 * i + 1))
    mycos = lambda x: somatorio(lambda i: (x ** (2 * i)) * (-1) ** i / factorial(2 * i))
    mytan = lambda x: mysin(x) / mycos(x)
    mypi = somatorio(lambda k: 4 * (-1) ** k / (2 * k + 1))
    mye = somatorio(lambda k: 1.0 / factorial(k))
    mypie = mypi / mye

    # calcular valores
    valores = [
        ['n={0}'.format(n), 'Ve', 'Va', 'Eabs', 'Erel'],
        ['pi', pi, mypi],
        ['e', e, mye],
        ['pi/e', pie, mypie],
        ['sen(pi/e)', sin(pie), mysin(mypie)],
        ['cos(pi/e)', cos(pie), mycos(mypie)],
        ['tg(pi/e)', tan(pie), mytan(mypie)],
    ]

    # calcular erros
    for v in valores[1:]:
        ve, va = v[1], v[2]
        eabs = ve - va
        erel = eabs / ve
        v.append(eabs)
        v.append(erel)

    # imprimir tabela e uma linha em branco
    print_table(valores)
    print

def naiveBayesMultinomial(value,frequency,adjustment = 1):
    if not value: return 0
    lim = 80
    if frequency > 92:
        return math.log(math.pow(value*adjustment,lim)/math.factorial(lim))
    else:
        return math.log(math.pow(value*adjustment,frequency)/math.factorial(frequency))

def nSphereVolume(n):
  '''
  Calculate the volume of a unit n-sphere. 

  TYPICAL USAGE
  =============
  
    >>> nSphereVolume(0)
    1
    >>> nSphereVolume(1)
    2
    >>> nSphereVolume(2)/pi
    1.0
    >>> nSphereVolume(3)/(pi**2)
    0.5
  
  REFERENCES
  ==========
  
    Equation 5.19.4, NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.0.6 of 2013-05-06.
  
  @param n: dimension of the n-sphere
  @type n: int
  @return: volume of the unit n-sphere
  @rtype: 
  '''
  if n%2 == 0:
    d = n/2
    return(pi**d/factorial(d))
  else:
    d = (n-1)/2
    return((2*((4*pi)**d)*factorial(d))/(factorial(2*d+1)))