Python math 模块，factorial() 实例源码

def test_figurate(self):
        """Test that figurate outputs correct values for, n = 0..nmax, d = 1..dmax
        by comparing against evaluation using explicit formulae for linear, triangular 
        and tetrahedral numbers, and also evaluation using math.factorial."""
        nmax = 20
        dmax = 20
        for n in range(0,nmax):
            # Linear numbers
            self.assertEqual( n,figurate(n,1) )
            # Triangular numbers
            self.assertEqual( int( round( n*(n+1)/2 ) ), figurate(n,2) )
            # Tetrahedral numbers
            self.assertEqual( int( round( n*(n+1)*(n+2)/6 ) ), figurate(n,3) )
            # General
            for d in range(1,dmax):
                self.assertEqual( self.binomial_factorial(n+d-1,d), figurate(n,d) )

def _get_M(self,delta_t):
        n=len(self.a)
        A=np.zeros(n)
        for i,ai in enumerate(self.a):
            Ae=[self.a[i]*(-1)**j*math.factorial(i)/math.factorial(j)/math.factorial(i-j)/((delta_t)**i) for j in range(i+1)]# Elementery A to assemblate in A
            for j,aej in enumerate(Ae):
                A[j]+=aej

        n=len(self.b)
        B=np.zeros(n)
        for i,ai in enumerate(self.b):
            Be=[self.b[i]*(-1)**j*math.factorial(i)/math.factorial(j)/math.factorial(i-j)/((delta_t)**i) for j in range(i+1)]# Elementery A to assemblate in A
            for j,bej in enumerate(Be):
                B[j]+=bej

        Mo=[-x/B[0] for x in B[1:][::-1]]
        Mi=[x/B[0] for x in A[::-1]]
        return (Mi,Mo)

def __init__(self, index):
        self.points = numpy.linspace(-1.0, 1.0, index+1)
        self.degree = index + 1 if index % 2 == 0 else index

        # Formula (26) from
        # <http://mathworld.wolfram.com/Newton-CotesFormulas.html>.
        # Note that Sympy carries out all operations in rationals, i.e.,
        # _exactly_. Only at the end, the rational is converted into a float.
        n = index
        self.weights = numpy.empty(n+1)
        t = sympy.Symbol('t')
        for r in range(n+1):
            # Compare with get_weights().
            f = sympy.prod([(t - i) for i in range(n+1) if i != r])
            alpha = 2 * \
                (-1)**(n-r) * sympy.integrate(f, (t, 0, n)) \
                / (math.factorial(r) * math.factorial(n-r)) \
                / index
            self.weights[r] = alpha
        return

def __init__(self, index):
        self.points = numpy.linspace(-1.0, 1.0, index+2)[1:-1]
        self.degree = index if (index+1) % 2 == 0 else index - 1
        #
        n = index+1
        self.weights = numpy.empty(n-1)
        t = sympy.Symbol('t')
        for r in range(1, n):
            # Compare with get_weights().
            f = sympy.prod([(t - i) for i in range(1, n) if i != r])
            alpha = 2 * \
                (-1)**(n-r+1) * sympy.integrate(f, (t, 0, n)) \
                / (math.factorial(r-1) * math.factorial(n-1-r)) \
                / n
            self.weights[r-1] = alpha
        return

def poisson_distribution(gamma):
    """Computes the probability of events for a Poisson distribution.

    Args:
        gamma: The average number of events to occur in an interval.

    Returns:
        The probability distribution of k events occuring.
        This is a function taking one parameter (k) and returning the
        probability of k events occuring.
    """
    constant_factor = math.e ** (-gamma)

    def probability(k):
        """The probability of k events occuring for a Poisson distribution.

        Args:
            k: The number of the events occuring in an interval.

        Returns:
            The probability of k events occuring in the given distribution.
        """
        return (constant_factor * gamma ** k) / (math.factorial(k))

    return probability

def getPermutation(self, n, k):
        """
        :type n: int
        :type k: int
        :rtype: str
        """
        k = k-1
        L = [l for l in range(1,n+1)]
        S = ""
        for l in xrange(n):
            base = factorial(n-l-1)
            ind = k // base
            k = k % base
            S += str(L[ind])
            L.remove(L[ind])



        return S

def non_increasing(n_digit=100):
    # a decreasing number can be obtained by reversing an increasing number
    # to consider possible tailing zeros of a decreasing number, I dissect
    # leading part of the coded binary strings into several cases:
    # 0, 10, 110, 1110, 11110, ...
    r = 8
    ori = n_digit + 8
    m = 1
    total = 0
    while ori >= 9:
        cur = factorial(ori) / factorial(r) / factorial(ori-r)
        total += (cur*m)
        m += 1
        ori -= 1

    return total

def savitzky_golay(y, window_size, order, deriv=0, rate=1):
    import numpy as np
    from math import factorial
    try:
        window_size = np.abs(np.int(window_size))
        order = np.abs(np.int(order))
    except ValueError, msg:
        raise ValueError("window_size and order have to be of type int")
    if window_size % 2 != 1 or window_size < 1:
        raise TypeError("window_size size must be a positive odd number")
    if window_size < order + 2:
        raise TypeError("window_size is too small for the polynomials order")
    order_range = range(order+1)
    half_window = (window_size -1) // 2
    b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
    m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
    firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
    lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
    y = np.concatenate((firstvals, y, lastvals))
    return np.convolve( m[::-1], y, mode='valid')

def savitzky_golay(y, window_size, order, deriv=0, rate=1):
    import numpy as np
    from math import factorial
    try:
        window_size = np.abs(np.int(window_size))
        order = np.abs(np.int(order))
    except ValueError, msg:
        raise ValueError("window_size and order have to be of type int")
    if window_size % 2 != 1 or window_size < 1:
        raise TypeError("window_size size must be a positive odd number")
    if window_size < order + 2:
        raise TypeError("window_size is too small for the polynomials order")
    order_range = range(order+1)
    half_window = (window_size -1) // 2
    b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
    m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
    firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
    lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
    y = np.concatenate((firstvals, y, lastvals))
    return np.convolve( m[::-1], y, mode='valid')

def savitzky_golay(y, window_size, order, deriv=0, rate=1):
    import numpy as np
    from math import factorial
    try:
        window_size = np.abs(np.int(window_size))
        order = np.abs(np.int(order))
    except ValueError, msg:
        raise ValueError("window_size and order have to be of type int")
    if window_size % 2 != 1 or window_size < 1:
        raise TypeError("window_size size must be a positive odd number")
    if window_size < order + 2:
        raise TypeError("window_size is too small for the polynomials order")
    order_range = range(order+1)
    half_window = (window_size -1) // 2
    b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
    m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
    firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
    lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
    y = np.concatenate((firstvals, y, lastvals))
    return np.convolve( m[::-1], y, mode='valid')

def test_quartic(self):
        import math

        # Use a fixed seed so the test is deterministic and round
        # the nodes to 8 bits of precision to avoid round-off.
        nodes = utils.get_random_nodes(
            shape=(15, 2), seed=11112222, num_bits=8)

        lambda_vals = (0.125, 0.375, 0.5)
        index = 0
        expected = np.asfortranarray([[0.0, 0.0]])
        for k in range(4 + 1):
            for j in range(4 + 1 - k):
                i = 4 - j - k
                denom = (math.factorial(i) * math.factorial(j) *
                         math.factorial(k))
                coeff = 24 / denom
                expected += (
                    coeff * lambda_vals[0]**i * lambda_vals[1]**j *
                    lambda_vals[2]**k * nodes[[index], :])
                index += 1

        result = self._call_function_under_test(nodes, 4, *lambda_vals)
        self.assertEqual(result, expected)

def perm(n):
  """Permutations.
     n = number of elements

     Notation: P
                n

     All elements included: yes
     Can elements repeat:   no
     Order matters:         yes

     See: Number of ways there are to rearrange n elements.

     Practical example: amount of numbers with 3 distinct digits.
     Let the elements be: 1, 2, 3:

       123, 132, 213, 231, 312, 321
  """
  return factorial(n)

def dpois(lmbda):
  """Poisson Distribution
     lmbda = average number of successes per unit interval

     Used to determine the probability of an amount of
     successes occuring in a fixed interval (time, area…)

     This doesn't return a value, but rather the specified Poisson function
  """
  def p(k):
    if 0 <= k:
      return (exp(-lmbda) * lmbda**k) / factorial(k)
    else:
      return 0

  # Allow accessing the used 'lmbda' value from the function
  p.__dict__['lmbda'] = lmbda
  p.__dict__['expected'] = lmbda
  p.__dict__['variance'] = lmbda
  return p

def getPermutation(self, n, k):
        """
        :type n: int
        :type k: int
        :rtype: str
        """
        import math
        f = math.factorial(n)
        l = [i + 1 for i in range(n)]
        ans = ""
        while l:
            f /= len(l)
            for i in range(len(l)):
                if i * f < k <= (i + 1) * f:
                    ans += str(l[i])
                    l = l[0:i] + l[i + 1:]
                    k -= (i * f)
        return ans

def savitzky_golay(y, window_size, order, deriv=0, rate=1):


    try:
        window_size = np.abs(np.int(window_size))
        order = np.abs(np.int(order))
    except ValueError, msg:
        raise ValueError("window_size and order have to be of type int")
    if window_size % 2 != 1 or window_size < 1:
        raise TypeError("window_size size must be a positive odd number")
    if window_size < order + 2:
        raise TypeError("window_size is too small for the polynomials order")
    order_range = range(order+1)
    half_window = (window_size -1) // 2
    # precompute coefficients
    b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
    m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
    # pad the signal at the extremes with
    # values taken from the signal itself
    firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
    lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
    y = np.concatenate((firstvals, y, lastvals))
    return np.convolve( m[::-1], y, mode='valid')

def solve_rpn(expression):
    operations = {
        '+': lambda x, y: x + y,
        '-': lambda x, y: x - y,
        '*': lambda x, y: x * y,
        '/': lambda x, y: x / y,
        '//': lambda x, y: x // y,
        '%': lambda x, y: x % y,
        '^': lambda x, y: x ** y
    }

    stack = collections.deque()

    for token in expression.split():
        if token == '!':
            stack.append(math.factorial(stack.pop()))
        elif token in operations:
            arguments = [stack.pop(), stack.pop()][::-1]
            stack.append(operations[token](*arguments))
        else:
            stack.append(float(token))

    return stack[0]

def findPermutation(s, rank):
    n = len(s)
    step = factorial(n)
    perm = ''
    s = sorted(s)
    # if rank == 0:
    #     return s
    index = 0
    for i in range(0, n):
        step /= (n - i)
        index = int((rank) // step)
        # if index > 0:
        perm += s[index]
        s = s[:index] + s[index+1:]
        rank -= (index) * step
        # else:  # index == 0
        #     perm += s[0]
        #     s = s[1:]
        # print(perm, '\t***', rank)
    return perm

def listPosition(word):
    s = word
    n = len(s)
    # print(factorial(n))
    rem = 1
    raw_order = sorted(s)
    # print(raw_order)
    for idx, letter in enumerate(s):
        l = len(raw_order)
        F = factorial(l)//repeat_product(raw_order)
        position = raw_order.index(letter)
        if position == 0:
            raw_order = raw_order[1:]
            continue
        rem += F//l*position
        raw_order = raw_order[:position] + raw_order[position+1:]
    return rem

def on_message(message):
    if message.author.bot:
        return
    nums = locate_numbers(message.content)
    if len(nums) > 0:
        guild = message.guild
        oh_no = "You've uttered "
        if guild is not None:
            name = message.author.nick or message.author.name
            oh_no = name + " has uttered "
        oh_no += "some factorials!" if len(nums) > 1 else "a factorial!"
        maxed = len(nums) > max_mentions
        if maxed:
            nums = nums[:max_mentions]
        for figure in nums:
            oh_no += "".join(
                ["\n``", figure[0], str(figure[1]), get_factorial(figure),
                 "``"])
        if maxed:
            oh_no += "\nand others..."
        await message.channel.send(oh_no)

def generateSampleNumberForBins(num_per_coder, n, k):
    """
        Generate the number of articles necessary to assign each coder * number
        of articles, given the bin size.
        Formula is: (n - 1)! / (k - 1)! (n - 1)!
    """

    a = factorial(n - 1) / ( factorial( k - 1 ) * factorial(n - k) )

    ## the number per coder doesn't divide evenly into the number of appearances
    ## subtract the remainder
    remainder = num_per_coder % a
    num_per_coder -= remainder

    ## get number of bins
    num_bins = len(list(combinations(range(0,n), k)))

    return int(num_bins * num_per_coder / a)

#####
### DUPE HANDLING
#####

def getPermutation(self, n, k):
        """
        :type n: int
        :type k: int
        :rtype: str
        """
        numbers = range(1, n+1)
        permutation = ''
        k -= 1
        while n > 0:
            n -= 1
            # get the index of current digit
            index, k = divmod(k, math.factorial(n))
            permutation += str(numbers[index])
            # remove handled number
            numbers.remove(numbers[index])

        return permutation

def binom(n: int, k: int) -> int:
    '''Computes the binomial coefficient.
    This shows how many combinations of
    *n* things taken in groups of size *k*.

    :param n: size of the universe
    :param k: size of each subset

    :returns: the number of combinations

    >>> binom(52, 5)
    2598960
    >>> binom(52, 0)
    1
    >>> binom(52, 52)
    1
    '''
    return factorial(n) // (factorial(k) * factorial(n-k))

def stirling2(n, k):
    """

    The Stirling numbers of the second kind,
    written S(n,k) or :math:`\lbrace\textstyle{n\atop k}\rbrace`
    count the number of ways to partition a set of n labelled objects
    into k nonempty unlabelled subsets.

    ..  math::

        \lbrace\textstyle{n\atop n}\rbrace = 1  \\

        \lbrace\textstyle{n\atop 1}\rbrace = 1  \\

        \lbrace\textstyle{n\atop k}\rbrace = k \lbrace\textstyle{n-1 \atop k}\rbrace + \lbrace\textstyle{n-1 \atop k-1}\rbrace

    Or

    ..  math::

        \left\{ {n \atop k}\right\} = \frac{1}{k!}\sum_{j=0}^{k} (-1)^{k-j} \binom{k}{j} j^n
    """

    return 1/factorial(k)*sum( (-1 if (k-j)%2 else 1)*binom(k,j)*j**n for j in range(0,k+1) )

def getPermutation(self, n, k):
        """
        :type n: int
        :type k: int
        :rtype: str
        """
        if n <= 0 or k <= 0:
            return ""
        from math import factorial, ceil
        solution = []
        digits = [str(i) for i in range(0, n + 1)]
        while n > 0:
            n2 = factorial(n - 1)
            index = int(ceil(1.0 * k / n2))
            solution.append(digits.pop(index))
            n -= 1
            k -= n2 * (index - 1)
        return ''.join(solution)

def findPermutation(s, rank):
    n = len(s)
    step = factorial(n)
    perm = ''
    s = sorted(s)
    # if rank == 0:
    #     return s
    index = 0
    for i in range(0, n):
        step /= (n - i)
        index = int((rank) // step)
        # if index > 0:
        perm += s[index]
        s = s[:index] + s[index+1:]
        rank -= (index) * step
        # else:  # index == 0
        #     perm += s[0]
        #     s = s[1:]
        # print(perm, '\t***', rank)
    return perm

def listPosition(word):
    s = word
    n = len(s)
    # print(factorial(n))
    rem = 1
    raw_order = sorted(s)
    # print(raw_order)
    for idx, letter in enumerate(s):
        l = len(raw_order)
        F = factorial(l)//repeat_product(raw_order)
        position = raw_order.index(letter)
        if position == 0:
            raw_order = raw_order[1:]
            continue
        rem += F//l*position
        raw_order = raw_order[:position] + raw_order[position+1:]
    return rem

def incomplete_gamma_lower(a, x, precision=9):
    # Numerical approximation of the lower incomplete gamma function to the given precision
    # https://en.wikipedia.org/wiki/Incomplete_gamma_function
    # http://mathworld.wolfram.com/IncompleteGammaFunction.html

    max_diff = 10 ** -precision

    def summand(k):
        return ((-x) ** k) / factorial(k) / (a + k)

    xs = x ** a
    sum_ = summand(0)
    prev_value = xs * sum_  # for k=0
    sum_ += summand(1)
    cur_value = xs * sum_  # for k=1
    k = 1
    while abs(cur_value - prev_value) > max_diff:
        k += 1
        prev_value = cur_value
        sum_ += summand(k)
        cur_value = xs * sum_

    return cur_value

def table_es(nstate,nparticle,spinz,dtype=np.int64):
    '''
    This function generates the table of binary representations of a particle-conserved and spin-conserved basis.
    '''
    n,nup,ndw=nstate/2,(nparticle+int(2*spinz))/2,(nparticle-int(2*spinz))/2
    result=np.zeros(factorial(n)/factorial(nup)/factorial(n-nup)*factorial(n)/factorial(ndw)/factorial(n-ndw),dtype=dtype)
    buff_up=list(combinations(xrange(1,2*n,2),nup))
    buff_dw=list(combinations(xrange(0,2*n,2),ndw))
    count=0
    for vup in buff_up:
        buff=0
        for num in vup:
            buff+=(1<<num)
        for vdw in buff_dw:
            basis=buff
            for num in vdw:
                basis+=(1<<num)
            result[count]=basis
            count+=1
    result.sort()
    return result

def set_parameters(self, X):
        # number of possible rules, i.e. rule space italic(A) prior
        self.patternSpace = np.ones(self.maxlen+1)
        # This patternSpace is an approximation
        # because the original code allows
        # the following situation, take tic-tac-toe
        # 1_O == 1 and 1_O_neg == 1, which is impossible
        numAttributes = len(X.columns)
        for i in range(1, self.maxlen+1):
            tmp = 1
            for j in range(numAttributes-i+1,numAttributes+1): 
                tmp *= j
            self.patternSpace[i] = tmp / math.factorial(i)

        if self.alpha_l == None:
            self.alpha_l = [1 for i in range(self.maxlen+1)]
        if self.beta_l == None:
            self.beta_l = [(self.patternSpace[i]*100+1) for i in range(self.maxlen+1)]

def set_parameters(self, X):
        # number of possible rules, i.e. rule space italic(A) prior
        self.patternSpace = np.ones(self.maxlen+1)
        # This patternSpace is an approximation
        # because the original code allows
        # the following situation, take tic-tac-toe
        # 1_O == 1 and 1_O_neg == 1, which is impossible
        numAttributes = len(X.columns)
        for i in range(1, self.maxlen+1):
            tmp = 1
            for j in range(numAttributes-i+1,numAttributes+1): 
                tmp *= j
            self.patternSpace[i] = tmp / math.factorial(i)

        if self.alpha_l == None:
            self.alpha_l = [1 for i in range(self.maxlen+1)]
        if self.beta_l == None:
            self.beta_l = [(self.patternSpace[i]*100+1) for i in range(self.maxlen+1)]

def test_perm(self):
        ''' Permutations. '''
        fs_ord = set()
        fs_unord = set()
        for fs in util.factorize(512, 3):
            fs_ord.add(fs)
            fs_unord.add(frozenset(fs))

        cnt = 0
        for fs in fs_unord:
            if len(fs) == 3:
                # Permutations.
                cnt += math.factorial(3)
            elif len(fs) == 2:
                # Permutations of a, a, b.
                cnt += 3
            else:
                # Pattern a, a, a.
                cnt += 1
        self.assertEqual(len(fs_ord), cnt)

def combinations(n, r):
    return fact(n) / fact(r) / fact(n-r)

def combination(n, k):
    """ Code for combination formula i.e. "n choose k" """
    from math import factorial
    return factorial(n) / (factorial(n - k) * factorial(k))

def __init__(self, margin, beta, beta_min, scale):
        self.margin = int(margin)
        self.beta = float(beta)
        self.beta_min = float(beta_min)
        self.scale = float(scale)
        self.c_map = []
        self.k_map = []
        c_m_n = lambda m, n: math.factorial(n) / math.factorial(m) / math.factorial(n-m)
        for i in range(margin+1):
            self.c_map.append(c_m_n(i, margin))
            self.k_map.append(math.cos(i * math.pi / margin))

def estimate_event_probability(self, r, n, p):
        # Binomial Event Discriminator
        from_timestamp = self.min_timestamp + datetime.timedelta(days=365)
        to_timestamp = self.max_timestamp
        timestamps, values = self.load_monitor_data(from_timestamp, to_timestamp, None)
        values = np.array(values, dtype=float)
        prob = math.factorial(n) / (math.factorial(r) * math.factorial(n-r)) * math.pow(p, r) * (math.pow(1-p, n-r))

def fOperator(stack, z, mode):
    if mode == 1:   # num
        stack.append(math.factorial(int(z)))
    elif mode == 2: # str
        stack.append([''.join(p) for p in itertools.permutations(z)])
    elif mode == 3: # list
        stack.append([list(p) for p in itertools.permutations(z)])
    else:
        monadNotImplemented(mode, '')

# g

def KOperator(stack, x, y, mode):
    if mode == 1:   # num, num
        n = int(x)
        k = int(y)
        if k < 0 or k > n:
            stack.append(0)
        else:
            stack.append(math.factorial(n)/(math.factorial(k)*math.factorial(n-k)))
    #elif mode == 2: # num, str
    elif mode == 3 or mode == 7: # num, list
        n = int(x if mode == 3 else y)
        l = y if mode == 3 else x

        def subsets(l, n):
            if n > len(l) or n < 0:
                return []
            elif n == len(l):
                return [l]
            elif n == 0:
                return [[]]
            elif n == 1:
                return [[i] for i in l]
            else:
                result = []
                for i in range(len(l)-n+1):
                    result += [[l[i]] + s for s in subsets(l[i+1:], n-1)]
                return result

        stack.append(subsets(l, n))
    #elif mode == 4: # str, num
    elif mode == 5: # str, str
        stack.append(''.join(c for c in x if c in y))
    #elif mode == 6: # str, list
    #elif mode == 8: # list, str
    elif mode == 9: # list, list
        stack.append([i for i in x if i in y])
    else:
        dyadNotImplemented(mode, '')

# ?

def _hypervolume_couboid(self, cuboid):
        """Computes the hypervolume of a single fuzzified cuboid."""

        all_dims = [dim for domain in self._core._domains.values() for dim in domain]
        n = len(all_dims)

        # calculating the factor in front of the sum
        weight_product = 1.0
        for (dom, dom_weight) in self._weights._domain_weights.items():
            for (dim, dim_weight) in self._weights._dimension_weights[dom].items():
                weight_product *= dom_weight * sqrt(dim_weight)
        factor = self._mu / (self._c**n * weight_product)

        # outer sum
        outer_sum = 0.0        
        for i in range(0, n+1):
            # inner sum
            inner_sum = 0.0
            subsets = list(itertools.combinations(all_dims, i))
            for subset in subsets:
                # first product
                first_product = 1.0
                for dim in set(all_dims) - set(subset):
                    dom = filter(lambda (x,y): dim in y, self._core._domains.items())[0][0]
                    w_dom = self._weights._domain_weights[dom]
                    w_dim = self._weights._dimension_weights[dom][dim]
                    b = cuboid._p_max[dim] - cuboid._p_min[dim]
                    first_product *= w_dom * sqrt(w_dim) * b * self._c

                # second product
                second_product = 1.0
                reduced_domain_structure = self._reduce_domains(self._core._domains, subset)
                for (dom, dims) in reduced_domain_structure.items():
                    n_domain = len(dims)
                    second_product *= factorial(n_domain) * (pi ** (n_domain/2.0))/(gamma((n_domain/2.0) + 1))

                inner_sum += first_product * second_product

            outer_sum += inner_sum
        return factor * outer_sum

def count_n_choose_k(n, k):
    """Returns the number of k combinations from the set n (n choose k)."""
    return (math.factorial(n) /
            math.factorial(k) /
            math.factorial(n - k))

def binomial_factorial(self,n,k):
        """Calculate binomial coefficient using math.factorial, for testing against binomial
        coefficients generated by other means."""
        if n >= k: 
            return int( round( math.factorial(n) / ( math.factorial(k) * math.factorial(n - k) ) ) )
        else:
            return 0

def test_binomial(self):
        """Test that binomial outputs correct values for, n = 0..nmax, k = 0, n
        by comparing against evaluation of math.factorial."""
        nmax = 20
        for n in range(0,nmax):
            for k in range(0,n):
                self.assertEqual( binomial(n,k), self.binomial_factorial(n,k) )

def k_comb(n, k):
    assert n > k, "Check values in k_comb!"
    return int(math.factorial(n) / (math.factorial(k) * math.factorial(n - k)))

def test_scheme(scheme, tol):
    degree = check_degree(
            lambda poly: quadpy.e1r.integrate(poly, scheme),
            lambda k: math.factorial(k[0]),
            1,
            scheme.degree + 1,
            tol=tol
            )
    assert degree == scheme.degree, \
        'Observed: {}   expected: {}'.format(degree, scheme.degree)
    return

def integrate_monomial_over_enr(k):
    if numpy.any(k % 2 == 1):
        return 0
    n = len(k)
    return 2 * math.factorial(sum(k) + n - 1) * numpy.prod([
        math.gamma((kk+1) / 2.0) for kk in k
        ]) / math.gamma((sum(k) + n) / 2)

def plot(scheme, show_axes=True):
    ax = plt.gca()
    plt.axis('equal')

    if not show_axes:
        ax.set_axis_off()

    n = 2
    I0 = 2*math.factorial(n-1)*math.pi**(0.5*n) / math.gamma(0.5*n)

    helpers.plot_disks(
        plt, scheme.points, scheme.weights, I0
        )
    return

def _size_cond(self, size):
        n = int(size)
        # r = 2
        # f = math.factorial
        # return int(f(n) / f(r) / f(n - r))
        return int((n * (n - 1)) / 2)

def size_cond(size):
    n = size
    r = 2
    f = math.factorial
    return int(f(n) / f(r) / f(n-r))

def nCr(n,r):
    f = math.factorial
    return f(n) / f(r) / f(n-r)

def facSum(n):
    sn = str(n)
    sumfac = 0
    for val in sn:
        valInt = int(val)
        sumfac += factorial(valInt)
    return sumfac

def sumOfDigFac(n):
    stringn = str(n)
    sumoffac = 0
    for i in xrange(0, len(stringn)):
        sumoffac += math.factorial(int(stringn[i]))
    if sumoffac != n:
        return False
    return True