Python sympy.factorial() Examples

def probability(self, event):
    # Specializations for optimization.
    if isinstance(event, FiniteProductEvent):
      assert len(self._spaces) == len(event.events)
      return sympy.prod([
          space.probability(event_slice)
          for space, event_slice in zip(self._spaces, event.events)])

    if isinstance(event, CountLevelSetEvent) and self.all_spaces_equal():
      space = self._spaces[0]
      counts = event.counts
      probabilities = {
          value: space.probability(DiscreteEvent({value}))
          for value in six.iterkeys(counts)
      }

      num_events = sum(six.itervalues(counts))
      assert num_events == len(self._spaces)
      # Multinomial coefficient:
      coeff = (
          sympy.factorial(num_events) / sympy.prod(
              [sympy.factorial(i) for i in six.itervalues(counts)]))
      return coeff * sympy.prod([
          pow(probabilities[value], counts[value])
          for value in six.iterkeys(counts)
      ])

    raise ValueError('Unhandled event type {}'.format(type(event))) 

def test_factorial():
    n = sympy.Symbol('n')
    assert theano_code(sympy.factorial(n)) 

def compute_dobrodeev(n, I0, I2, I22, I4, pm_type, i, j, k, symbolic=False):
    """Compute some helper quantities used in

    L.N. Dobrodeev,
    Cubature rules with equal coefficients for integrating functions with
    respect to symmetric domains,
    USSR Computational Mathematics and Mathematical Physics,
    Volume 18, Issue 4, 1978, Pages 27-34,
    <https://doi.org/10.1016/0041-5553(78)90064-2>.
    """
    t = 1 if pm_type == "I" else -1

    fact = sympy.factorial if symbolic else math.factorial
    sqrt = sympy.sqrt if symbolic else numpy.sqrt

    L = comb(n, i) * 2 ** i
    M = fact(n) // (fact(j) * fact(k) * fact(n - j - k)) * 2 ** (j + k)
    N = L + M
    F = I22 / I0 - I2 ** 2 / I0 ** 2 + (I4 / I0 - I22 / I0) / n
    R = (
        -(j + k - i) / i * I2 ** 2 / I0 ** 2
        + (j + k - 1) / n * I4 / I0
        - (n - 1) / n * I22 / I0
    )
    H = (
        1
        / i
        * (
            (j + k - i) * I2 ** 2 / I0 ** 2
            + (j + k) / n * ((i - 1) * I4 / I0 - (n - 1) * I22 / I0)
        )
    )
    Q = L / M * R + H - t * 2 * I2 / I0 * (j + k - i) / i * sqrt(L / M * F)

    G = 1 / N
    a = sqrt(n / i * (I2 / I0 + t * sqrt(M / L * F)))
    b = sqrt(n / (j + k) * (I2 / I0 - t * sqrt(L / M * F) + t * sqrt(k / j * Q)))
    c = sqrt(n / (j + k) * (I2 / I0 - t * sqrt(L / M * F) - t * sqrt(j / k * Q)))
    return G, a, b, c 

def comb(a, b):
    if sys.version < "3.8":
        try:
            binom = math.factorial(a) // math.factorial(b) // math.factorial(a - b)
        except ValueError:
            binom = 0
        return binom
    return math.comb(a, b) 

def gamma_n_2(n, symbolic):
    # gamma(n / 2)
    frac = sympy.Rational if symbolic else lambda a, b: a / b
    sqrt = sympy.sqrt if symbolic else math.sqrt
    pi = sympy.pi if symbolic else math.pi

    if n % 2 == 0:
        return math.factorial(n // 2 - 1)

    n2 = n // 2
    return frac(math.factorial(2 * n2), 4 ** n2 * math.factorial(n2)) * sqrt(pi) 

def __init__(self, **traits):
        TaylorPoly.__init__(self, **traits)
        #Declare the analytical function
  
        Z=sympy.Function("Z")
        Ax,Ay,Kx,Ky=sympy.symbols(("Ax","Ay","Kx","Ky"))
        x, y =sympy.symbols('xy')
        Z=(Ax*x**2+Ay*y**2)/(1+sympy.sqrt(1-(1+Kx)*Ax**2*x**2-(1+Ky)*Ay**2*y**2));
        
        #Calculate taylor polynomial coheficients
        cohef=[[Z, ],]
        order=self.n
        for i in range(0, order+1, 2):
            if i!=0:
                cohef.append([sympy.diff(cohef[i/2-1][0], y, 2), ])
            for j in range(2, order-i+1, 2):
                cohef[i/2].append(sympy.diff(cohef[i/2][j/2 -1], x, 2))
        

        A_x=self.Ax
        A_y=self.Ay
        K_x=self.Kx
        K_y=self.Ky
        
        c=zeros((self.n+1, self.n+1))
        for i in range(0, order/2+1):
            for j in range(0,order/2- i+1):
                cohef[j][i]=cohef[j][i].subs(x, 0).subs(y, 0).subs(Ax, A_x).subs(Ay, A_y).subs(Kx, K_x).subs(Ky, K_y)/(sympy.factorial(2*i)*sympy.factorial(2*j))
                c[2*j, 2*i]=cohef[j][i].evalf()
        
        # Add the high order corrections
        if len(self.ho_cohef.shape)==2:
            cx, cy = c.shape 
            dx, dy =self.ho_cohef.shape
            mx=array((cx, dx)).max()
            my=array((cy, dy)).max()
            self.cohef=zeros((mx, my))
            self.cohef[0:cx, 0:cy]=c
            self.cohef[0:dy, 0:dy]=self.cohef[0:dy, 0:dy]+self.ho_cohef
        else:
            self.cohef=c 

def convert_postfix(postfix):
    if hasattr(postfix, 'exp'):
        exp_nested = postfix.exp()
    else:
        exp_nested = postfix.exp_nofunc()

    exp = convert_exp(exp_nested)
    for op in postfix.postfix_op():
        if op.BANG():
            if isinstance(exp, list):
                raise Exception("Cannot apply postfix to derivative")
            exp = sympy.factorial(exp, evaluate=False)
        elif op.eval_at():
            ev = op.eval_at()
            at_b = None
            at_a = None
            if ev.eval_at_sup():
                at_b = do_subs(exp, ev.eval_at_sup()) 
            if ev.eval_at_sub():
                at_a = do_subs(exp, ev.eval_at_sub())
            if at_b != None and at_a != None:
                exp = sympy.Add(at_b, -1 * at_a, evaluate=False)
            elif at_b != None:
                exp = at_b
            elif at_a != None:
                exp = at_a
            
    return exp 

def _get_Ylm(self, l, m):
        """
        Compute an expression for spherical harmonic of order (l,m)
        in terms of Cartesian unit vectors, :math:`\hat{z}`
        and :math:`\hat{x} + i \hat{y}`

        Parameters
        ----------
        l : int
            the degree of the harmonic
        m : int
            the order of the harmonic; |m| < l

        Returns
        -------
        expr :
            a sympy expression that corresponds to the
            requested Ylm

        References
        ----------
        https://en.wikipedia.org/wiki/Spherical_harmonics
        """
        import sympy as sp

        # the relevant cartesian and spherical symbols
        x, y, z, r = sp.symbols('x y z r', real=True, positive=True)
        xhat, yhat, zhat = sp.symbols('xhat yhat zhat', real=True, positive=True)
        xpyhat = sp.Symbol('xpyhat', complex=True)
        phi, theta = sp.symbols('phi theta')
        defs = [(sp.sin(phi), y/sp.sqrt(x**2+y**2)),
                (sp.cos(phi), x/sp.sqrt(x**2+y**2)),
                (sp.cos(theta), z/sp.sqrt(x**2 + y**2 + z**2))
                ]

        # the cos(theta) dependence encoded by the associated Legendre poly
        expr = sp.assoc_legendre(l, m, sp.cos(theta))

        # the exp(i*m*phi) dependence
        expr *= sp.expand_trig(sp.cos(m*phi)) + sp.I*sp.expand_trig(sp.sin(m*phi))

        # simplifying optimizations
        expr = sp.together(expr.subs(defs)).subs(x**2 + y**2 + z**2, r**2)
        expr = expr.expand().subs([(x/r, xhat), (y/r, yhat), (z/r, zhat)])
        expr = expr.factor().factor(extension=[sp.I]).subs(xhat+sp.I*yhat, xpyhat)
        expr = expr.subs(xhat**2 + yhat**2, 1-zhat**2).factor()

        # and finally add the normalization
        amp = sp.sqrt((2*l+1) / (4*numpy.pi) * sp.factorial(l-m) / sp.factorial(l+m))
        expr *= amp

        return expr