Search code examples
pythonnumpyvectorizationnumpy-ndarraybinomial-theorem

Vectorized approach for calculating numerator for generalized binomial theorem?


I am trying to program the generalized binomial theorem, where n can be any rational number, using a vectorized approach. Formula attached as image below.

The numerator for each term is n, n×(n-1) , n×(n-1)×(n-2) and so on. I have assigned 0.5 to n to and am trying to generate 5 terms.

So far, I have an array of the products of the numerator: [ 0.5 -0.5 -1.5 -2.5 -3.5] using

def num_products(number_of_terms):
    r = np.arange(1,number_of_terms+1)
    num_prod = np.array(n-r+1)
    return num_prod

But want to create an array of the numerators for each terms, like this (w each item in array shown separated by commas):

[ 0.5, 0.5×-0.5, 0.5×-0.5×-1.5, 0.5×-0.5×-1.5×-2.5, 0.5×-0.5×-1.5×-2.5×-3.5]

Does anyone know how to do this using arrays (vectorized approach)? I am trying to make it very quick to compute the terms so I can have a larger number of terms and increase accuracy of the result.

Formula for generalized binomial theorem


Solution

  • Each term x*(x-1)*(x-2)*...*(x - n + 1) is known as a falling factorial. The wikipedia article also describes the rising factorial x*(x+1)*...*(x + n - 1). Some computational libraries include implementations of these. For example, mpmath has mpmath.ff and mpmath.rf.

    SciPy implements the rising factorial as scipy.special.poch. The falling factorial can be implemented in terms of the rising factorial, as in

    from scipy.special import poch
    
    
    def ff(x, m):
        return poch(x - m + 1, m)
    

    Because poch is implemented as a NumPy "ufunc", it handles broadcasting, and therefore so does ff. This means you can pass in an array of values for m to compute all the corresponding falling factorials with one call.

    For example, to get the first six numerator terms (including the initial term 1) of the generalized binomial with n = 0.5, call ff(0.5, np.arange(6)):

    In [38]: ff(0.5, np.arange(6))
    Out[38]: array([ 1.     ,  0.5    , -0.25   ,  0.375  , -0.9375 ,  3.28125])
    

    That is the same as [1, 0.5, 0.5×-0.5, 0.5×-0.5×-1.5, 0.5×-0.5×-1.5×-2.5, 0.5×-0.5×-1.5×-2.5×-3.5]:

    In [40]: [1, 0.5, 0.5*-0.5, 0.5*-0.5*-1.5, 0.5*-0.5*-1.5*-2.5, 0.5*-0.5*-1.5*-2.5*-3.5]
    Out[40]: [1, 0.5, -0.25, 0.375, -0.9375, 3.28125]
    

    So if you don't mind the dependence on SciPy, you can use ff defined above to do what you want.