Python function.expand_mul函数代码示例

 def _cholesky(self, hermitian=True):
     """Helper function of cholesky.
     Without the error checks.
     To be used privately.
     Implements the Cholesky-Banachiewicz algorithm.
     Returns L such that L*L.H == self if hermitian flag is True,
     or L*L.T == self if hermitian is False.
     """
     L = zeros(self.rows, self.rows)
     if hermitian:
         for i in range(self.rows):
             for j in range(i):
                 L[i, j] = (1 / L[j, j])*expand_mul(self[i, j] -
                     sum(L[i, k]*L[j, k].conjugate() for k in range(j)))
             Lii2 = expand_mul(self[i, i] -
                 sum(L[i, k]*L[i, k].conjugate() for k in range(i)))
             if Lii2.is_positive is False:
                 raise ValueError("Matrix must be positive-definite")
             L[i, i] = sqrt(Lii2)
     else:
         for i in range(self.rows):
             for j in range(i):
                 L[i, j] = (1 / L[j, j])*(self[i, j] -
                     sum(L[i, k]*L[j, k] for k in range(j)))
             L[i, i] = sqrt(self[i, i] -
                 sum(L[i, k]**2 for k in range(i)))
     return self._new(L)

 def _LDLdecomposition(self, hermitian=True):
     """Helper function of LDLdecomposition.
     Without the error checks.
     To be used privately.
     Returns L and D such that L*D*L.H == self if hermitian flag is True,
     or L*D*L.T == self if hermitian is False.
     """
     # https://en.wikipedia.org/wiki/Cholesky_decomposition#LDL_decomposition_2
     D = zeros(self.rows, self.rows)
     L = eye(self.rows)
     if hermitian:
         for i in range(self.rows):
             for j in range(i):
                 L[i, j] = (1 / D[j, j])*expand_mul(self[i, j] - sum(
                     L[i, k]*L[j, k].conjugate()*D[k, k] for k in range(j)))
             D[i, i] = expand_mul(self[i, i] -
                 sum(L[i, k]*L[i, k].conjugate()*D[k, k] for k in range(i)))
             if D[i, i].is_positive is False:
                 raise ValueError("Matrix must be positive-definite")
     else:
         for i in range(self.rows):
             for j in range(i):
                 L[i, j] = (1 / D[j, j])*(self[i, j] - sum(
                     L[i, k]*L[j, k]*D[k, k] for k in range(j)))
             D[i, i] = self[i, i] - sum(L[i, k]**2*D[k, k] for k in range(i))
     return self._new(L), self._new(D)

 def valid(x):
     # this is used to see if gen=x satisfies the
     # relational by substituting it into the
     # expanded form and testing against 0, e.g.
     # if expr = x*(x + 1) < 2 then e = x*(x + 1) - 2
     # and expanded_e = x**2 + x - 2; the test is
     # whether a given value of x satisfies
     # x**2 + x - 2 < 0
     #
     # expanded_e, expr and gen used from enclosing scope
     v = expanded_e.subs(gen, expand_mul(x))
     try:
         r = expr.func(v, 0)
     except TypeError:
         r = S.false
     if r in (S.true, S.false):
         return r
     if v.is_real is False:
         return S.false
     else:
         v = v.n(2)
         if v.is_comparable:
             return expr.func(v, 0)
         # not comparable or couldn't be evaluated
         raise NotImplementedError(
             'relationship did not evaluate: %s' % r)

def test_issue_11230():
    # a specific test that always failed
    a, b, f, k, l, i = symbols('a b f k l i')
    p = [a*b*f*k*l, a*i*k**2*l, f*i*k**2*l]
    R, C = cse(p)
    assert not any(i.is_Mul for a in C for i in a.args)

    # random tests for the issue
    from random import choice
    from sympy.core.function import expand_mul
    s = symbols('a:m')
    # 35 Mul tests, none of which should ever fail
    ex = [Mul(*[choice(s) for i in range(5)]) for i in range(7)]
    for p in subsets(ex, 3):
        p = list(p)
        R, C = cse(p)
        assert not any(i.is_Mul for a in C for i in a.args)
        for ri in reversed(R):
            for i in range(len(C)):
                C[i] = C[i].subs(*ri)
        assert p == C
    # 35 Add tests, none of which should ever fail
    ex = [Add(*[choice(s[:7]) for i in range(5)]) for i in range(7)]
    for p in subsets(ex, 3):
        p = list(p)
        was = R, C = cse(p)
        assert not any(i.is_Add for a in C for i in a.args)
        for ri in reversed(R):
            for i in range(len(C)):
                C[i] = C[i].subs(*ri)
        # use expand_mul to handle cases like this:
        # p = [a + 2*b + 2*e, 2*b + c + 2*e, b + 2*c + 2*g]
        # x0 = 2*(b + e) is identified giving a rebuilt p that
        # is now `[a + 2*(b + e), c + 2*(b + e), b + 2*c + 2*g]`
        assert p == [expand_mul(i) for i in C]

def sqrtdenest(expr, max_iter=3):
    """Denests sqrts in an expression that contain other square roots
    if possible, otherwise returns the expr unchanged. This is based on the
    algorithms of [1].

    Examples
    ========

    >>> from sympy.simplify.sqrtdenest import sqrtdenest
    >>> from sympy import sqrt
    >>> sqrtdenest(sqrt(5 + 2 * sqrt(6)))
    sqrt(2) + sqrt(3)

    See Also
    ========
    sympy.solvers.solvers.unrad

    References
    ==========
    [1] http://researcher.watson.ibm.com/researcher/files/us-fagin/symb85.pdf

    [2] D. J. Jeffrey and A. D. Rich, 'Symplifying Square Roots of Square Roots
    by Denesting' (available at http://www.cybertester.com/data/denest.pdf)

    """
    expr = expand_mul(sympify(expr))
    for i in range(max_iter):
        z = _sqrtdenest0(expr)
        if expr == z:
            return expr
        expr = z
    return expr

    def _eval_as_leading_term(self, x):
        from sympy import expand_mul, factor_terms

        old = self

        expr = expand_mul(self)
        if not expr.is_Add:
            return expr.as_leading_term(x)

        infinite = [t for t in expr.args if t.is_infinite]

        expr = expr.func(*[t.as_leading_term(x) for t in expr.args]).removeO()
        if not expr:
            # simple leading term analysis gave us 0 but we have to send
            # back a term, so compute the leading term (via series)
            return old.compute_leading_term(x)
        elif expr is S.NaN:
            return old.func._from_args(infinite)
        elif not expr.is_Add:
            return expr
        else:
            plain = expr.func(*[s for s, _ in expr.extract_leading_order(x)])
            rv = factor_terms(plain, fraction=False)
            rv_simplify = rv.simplify()
            # if it simplifies to an x-free expression, return that;
            # tests don't fail if we don't but it seems nicer to do this
            if x not in rv_simplify.free_symbols:
                if rv_simplify.is_zero and plain.is_zero is not True:
                    return (expr - plain)._eval_as_leading_term(x)
                return rv_simplify
            return rv

def sqrtdenest(expr, max_iter=3):
    """Denests sqrts in an expression that contain other square roots
    if possible, otherwise returns the expr unchanged. This is based on the
    algorithms of [1].

    Examples
    ========

    >>> from sympy.simplify.sqrtdenest import sqrtdenest
    >>> from sympy import sqrt
    >>> sqrtdenest(sqrt(5 + 2 * sqrt(6)))
    sqrt(2) + sqrt(3)

    See Also
    ========
    sympy.solvers.solvers.unrad

    References
    ==========
    [1] http://www.almaden.ibm.com/cs/people/fagin/symb85.pdf
    """
    expr = expand_mul(sympify(expr))
    for i in range(max_iter):
        z = _sqrtdenest0(expr)
        if expr == z:
            return expr
        expr = z
    return expr

 def _bell_poly(n, prev):
     s = 1
     a = 1
     for k in xrange(2, n + 1):
         a = a * (n - k + 1) // (k - 1)
         s += a * prev[k - 1]
     return expand_mul(_sym * s)

    def _bell_incomplete_poly(n, k, symbols):
        r"""
        The second kind of Bell polynomials (incomplete Bell polynomials).

        Calculated by recurrence formula:

        .. math:: B_{n,k}(x_1, x_2, \dotsc, x_{n-k+1}) =
                \sum_{m=1}^{n-k+1}
                \x_m \binom{n-1}{m-1} B_{n-m,k-1}(x_1, x_2, \dotsc, x_{n-m-k})

        where
            B_{0,0} = 1;
            B_{n,0} = 0; for n>=1
            B_{0,k} = 0; for k>=1

        """
        if (n == 0) and (k == 0):
            return S.One
        elif (n == 0) or (k == 0):
            return S.Zero
        s = S.Zero
        a = S.One
        for m in xrange(1, n - k + 2):
            s += a * bell._bell_incomplete_poly(
                n - m, k - 1, symbols) * symbols[m - 1]
            a = a * (n - m) / m
        return expand_mul(s)

 def _LDLdecomposition(self):
     """Helper function of LDLdecomposition.
     Without the error checks.
     To be used privately.
     """
     # https://en.wikipedia.org/wiki/Cholesky_decomposition#LDL_decomposition_2
     D = zeros(self.rows, self.rows)
     L = eye(self.rows)
     for i in range(self.rows):
         for j in range(i):
             L[i, j] = (1 / D[j, j])*expand_mul(self[i, j] - sum(
                 L[i, k]*L[j, k].conjugate()*D[k, k] for k in range(j)))
         D[i, i] = expand_mul(self[i, i] -
             sum(L[i, k]*L[i, k].conjugate()*D[k, k] for k in range(i)))
         if D[i, i].is_positive is False:
             raise ValueError("Matrix must be positive-definite")
     return self._new(L), self._new(D)

 def _cholesky(self):
     """Helper function of cholesky.
     Without the error checks.
     To be used privately.
     Implements the Cholesky-Banachiewicz algorithm.
     """
     L = zeros(self.rows, self.rows)
     for i in range(self.rows):
         for j in range(i):
             L[i, j] = (1 / L[j, j])*expand_mul(self[i, j] -
                 sum(L[i, k]*L[j, k].conjugate() for k in range(j)))
         Lii2 = expand_mul(self[i, i] -
             sum(L[i, k]*L[i, k].conjugate() for k in range(i)))
         if Lii2.is_positive is False:
             raise ValueError("Matrix must be positive-definite")
         L[i, i] = sqrt(Lii2)
     return self._new(L)

def convolution_fwht(a, b):
    """
    Performs dyadic (*bitwise-XOR*) convolution using Fast Walsh Hadamard
    Transform.

    The convolution is automatically padded to the right with zeros, as the
    *radix-2 FWHT* requires the number of sample points to be a power of 2.

    Parameters
    ==========

    a, b : iterables
        The sequences for which convolution is performed.

    Examples
    ========

    >>> from sympy import symbols, S, I
    >>> from sympy.discrete.convolutions import convolution_fwht

    >>> u, v, x, y = symbols('u v x y')
    >>> convolution_fwht([u, v], [x, y])
    [u*x + v*y, u*y + v*x]

    >>> convolution_fwht([2, 3], [4, 5])
    [23, 22]
    >>> convolution_fwht([2, 5 + 4*I, 7], [6*I, 7, 3 + 4*I])
    [56 + 68*I, -10 + 30*I, 6 + 50*I, 48 + 32*I]

    >>> convolution_fwht([S(33)/7, S(55)/6, S(7)/4], [S(2)/3, 5])
    [2057/42, 1870/63, 7/6, 35/4]

    References
    ==========

    .. [1] https://www.radioeng.cz/fulltexts/2002/02_03_40_42.pdf
    .. [2] https://en.wikipedia.org/wiki/Hadamard_transform

    """

    if not a or not b:
        return []

    a, b = a[:], b[:]
    n = max(len(a), len(b))

    if n&(n - 1): # not a power of 2
        n = 2**n.bit_length()

    # padding with zeros
    a += [S.Zero]*(n - len(a))
    b += [S.Zero]*(n - len(b))

    a, b = fwht(a), fwht(b)
    a = [expand_mul(x*y) for x, y in zip(a, b)]
    a = ifwht(a)

    return a

 def _combine_inverse(lhs, rhs):
     """
     Returns lhs - rhs, but treats oo like a symbol so oo - oo
     returns 0, instead of a nan.
     """
     from sympy.core.function import expand_mul
     from sympy.core.symbol import Dummy
     inf = (S.Infinity, S.NegativeInfinity)
     if lhs.has(*inf) or rhs.has(*inf):
         oo = Dummy('oo')
         reps = {
             S.Infinity: oo,
             S.NegativeInfinity: -oo}
         ireps = dict([(v, k) for k, v in reps.items()])
         eq = expand_mul(lhs.xreplace(reps) - rhs.xreplace(reps))
         if eq.has(oo):
             eq = eq.replace(
                 lambda x: x.is_Pow and x.base == oo,
                 lambda x: x.base)
         return eq.xreplace(ireps)
     else:
         return expand_mul(lhs - rhs)

 def eval(cls, p, q):
     from sympy.simplify.simplify import nsimplify
     if q.is_Number:
         float = not q.is_Rational
         pnew = expand_mul(p)
         if pnew.is_Number:
             float = float or not pnew.is_Rational
             if not float:
                 return pnew % q
             return Float(nsimplify(pnew) % nsimplify(q))
         elif pnew.is_Add and pnew.args[0].is_Number:
             r, p = pnew.as_two_terms()
             p += Mod(r, q)
     return Mod(p, q, evaluate=False)

def convolution_fft(a, b, dps=None):
    """
    Performs linear convolution using Fast Fourier Transform.

    Parameters
    ==========

    a, b : iterables
        The sequences for which convolution is performed.
    dps : Integer
        Specifies the number of decimal digits for precision.

    Examples
    ========

    >>> from sympy import S, I
    >>> from sympy.discrete.convolutions import convolution_fft

    >>> convolution_fft([2, 3], [4, 5])
    [8, 22, 15]
    >>> convolution_fft([2, 5], [6, 7, 3])
    [12, 44, 41, 15]
    >>> convolution_fft([1 + 2*I, 4 + 3*I], [S(5)/4, 6])
    [5/4 + 5*I/2, 11 + 63*I/4, 24 + 18*I]

    References
    ==========

    .. [1] https://en.wikipedia.org/wiki/Convolution_theorem
    .. [2] https://en.wikipedia.org/wiki/Discrete_Fourier_transform_(general%29

    """

    a, b = a[:], b[:]
    n = m = len(a) + len(b) - 1 # convolution size

    if n > 0 and n&(n - 1): # not a power of 2
        n = 2**n.bit_length()

    # padding with zeros
    a += [S.Zero]*(n - len(a))
    b += [S.Zero]*(n - len(b))

    a, b = fft(a, dps), fft(b, dps)
    a = [expand_mul(x*y) for x, y in zip(a, b)]
    a = ifft(a, dps)[:m]

    return a