Python Poly.subs方法代码示例

# 需要导入模块: from sympy.polys import Poly [as 别名]
# 或者: from sympy.polys.Poly import subs [as 别名]
def test_issue_8438():
    p = Poly([1, y, -2, -3], x).as_expr()
    roots = roots_cubic(Poly(p, x), x)
    z = -S(3)/2 - 7*I/2  # this will fail in code given in commit msg
    post = [r.subs(y, z) for r in roots]
    assert set(post) == \
    set(roots_cubic(Poly(p.subs(y, z), x)))
    # /!\ if p is not made an expression, this is *very* slow
    assert all(p.subs({y: z, x: i}).n(2, chop=True) == 0 for i in post)

# 需要导入模块: from sympy.polys import Poly [as 别名]
# 或者: from sympy.polys.Poly import subs [as 别名]
def test_roots_quartic():
    assert roots_quartic(Poly(x**4, x)) == [0, 0, 0, 0]
    assert roots_quartic(Poly(x**4 + x**3, x)) in [
        [-1, 0, 0, 0],
        [0, -1, 0, 0],
        [0, 0, -1, 0],
        [0, 0, 0, -1]
    ]
    assert roots_quartic(Poly(x**4 - x**3, x)) in [
        [1, 0, 0, 0],
        [0, 1, 0, 0],
        [0, 0, 1, 0],
        [0, 0, 0, 1]
    ]

    lhs = roots_quartic(Poly(x**4 + x, x))
    rhs = [S.Half + I*sqrt(3)/2, S.Half - I*sqrt(3)/2, S.Zero, -S.One]

    assert sorted(lhs, key=hash) == sorted(rhs, key=hash)

    # test of all branches of roots quartic
    for i, (a, b, c, d) in enumerate([(1, 2, 3, 0),
                                      (3, -7, -9, 9),
                                      (1, 2, 3, 4),
                                      (1, 2, 3, 4),
                                      (-7, -3, 3, -6),
                                      (-3, 5, -6, -4),
                                      (6, -5, -10, -3)]):
        if i == 2:
            c = -a*(a**2/S(8) - b/S(2))
        elif i == 3:
            d = a*(a*(3*a**2/S(256) - b/S(16)) + c/S(4))
        eq = x**4 + a*x**3 + b*x**2 + c*x + d
        ans = roots_quartic(Poly(eq, x))
        assert all(eq.subs(x, ai).n(chop=True) == 0 for ai in ans)

    # not all symbolic quartics are unresolvable
    eq = Poly(q*x + q/4 + x**4 + x**3 + 2*x**2 - Rational(1, 3), x)
    sol = roots_quartic(eq)
    assert all(verify_numerically(eq.subs(x, i), 0) for i in sol)
    z = symbols('z', negative=True)
    eq = x**4 + 2*x**3 + 3*x**2 + x*(z + 11) + 5
    zans = roots_quartic(Poly(eq, x))
    assert all([verify_numerically(eq.subs(((x, i), (z, -1))), 0) for i in zans])
    # but some are (see also issue 4989)
    # it's ok if the solution is not Piecewise, but the tests below should pass
    eq = Poly(y*x**4 + x**3 - x + z, x)
    ans = roots_quartic(eq)
    assert all(type(i) == Piecewise for i in ans)
    reps = (
        dict(y=-Rational(1, 3), z=-Rational(1, 4)),  # 4 real
        dict(y=-Rational(1, 3), z=-Rational(1, 2)),  # 2 real
        dict(y=-Rational(1, 3), z=-2))  # 0 real
    for rep in reps:
        sol = roots_quartic(Poly(eq.subs(rep), x))
        assert all([verify_numerically(w.subs(rep) - s, 0) for w, s in zip(ans, sol)])

# 需要导入模块: from sympy.polys import Poly [as 别名]
# 或者: from sympy.polys.Poly import subs [as 别名]
def max_onepiece(x, f: Poly, g: Poly, l, u):
    roots = sorted(set((f - g).real_roots()))
    new_polynomial_pieces = []
    new_bounds = [l]
    for r in roots:
        if l < r < u:
            m = (r + new_bounds[-1]) / 2
            if f.subs(x, m) >= g.subs(x, m):
                new_polynomial_pieces.append(f)
            else:
                new_polynomial_pieces.append(g)
            new_bounds.append(r)
    new_bounds.append(u)
    return PiecewisePolynomial(new_polynomial_pieces, new_bounds)

# 需要导入模块: from sympy.polys import Poly [as 别名]
# 或者: from sympy.polys.Poly import subs [as 别名]
def test_roots_quartic():
    assert roots_quartic(Poly(x ** 4, x)) == [0, 0, 0, 0]
    assert roots_quartic(Poly(x ** 4 + x ** 3, x)) in [[-1, 0, 0, 0], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, -1]]
    assert roots_quartic(Poly(x ** 4 - x ** 3, x)) in [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]

    lhs = roots_quartic(Poly(x ** 4 + x, x))
    rhs = [S.Half + I * sqrt(3) / 2, S.Half - I * sqrt(3) / 2, S.Zero, -S.One]

    assert sorted(lhs, key=hash) == sorted(rhs, key=hash)

    # test of all branches of roots quartic
    for i, (a, b, c, d) in enumerate(
        [(1, 2, 3, 0), (3, -7, -9, 9), (1, 2, 3, 4), (1, 2, 3, 4), (-7, -3, 3, -6), (-3, 5, -6, -4), (6, -5, -10, -3)]
    ):
        if i == 2:
            c = -a * (a ** 2 / S(8) - b / S(2))
        elif i == 3:
            d = a * (a * (3 * a ** 2 / S(256) - b / S(16)) + c / S(4))
        eq = x ** 4 + a * x ** 3 + b * x ** 2 + c * x + d
        ans = roots_quartic(Poly(eq, x))
        assert all(eq.subs(x, ai).n(chop=True) == 0 for ai in ans)

    # not all symbolic quartics are unresolvable
    eq = Poly(q * x + q / 4 + x ** 4 + x ** 3 + 2 * x ** 2 - Rational(1, 3), x)
    sol = roots_quartic(eq)
    assert all(test_numerically(eq.subs(x, i), 0) for i in sol)
    # but some are (see also iss 1890)
    raises(PolynomialError, lambda: roots_quartic(Poly(y * x ** 4 + x + z, x)))

# 需要导入模块: from sympy.polys import Poly [as 别名]
# 或者: from sympy.polys.Poly import subs [as 别名]
def gosper_term(f, n):
    r"""
    Compute Gosper's hypergeometric term for ``f``.

    Suppose ``f`` is a hypergeometric term such that:

    .. math::
        s_n = \sum_{k=0}^{n-1} f_k

    and `f_k` doesn't depend on `n`. Returns a hypergeometric
    term `g_n` such that `g_{n+1} - g_n = f_n`.

    **Examples**

    >>> from sympy.concrete.gosper import gosper_term
    >>> from sympy.functions import factorial
    >>> from sympy.abc import n

    >>> gosper_term((4*n + 1)*factorial(n)/factorial(2*n + 1), n)
    (-n - 1/2)/(n + 1/4)

    """
    r = hypersimp(f, n)

    if r is None:
        return None  # 'f' is *not* a hypergeometric term

    p, q = r.as_numer_denom()

    A, B, C = gosper_normal(p, q, n)
    B = B.shift(-1)

    N = S(A.degree())
    M = S(B.degree())
    K = S(C.degree())

    if (N != M) or (A.LC() != B.LC()):
        D = set([K - max(N, M)])
    elif not N:
        D = set([K - N + 1, S(0)])
    else:
        D = set([K - N + 1, (B.nth(N - 1) - A.nth(N - 1)) / A.LC()])

    for d in set(D):
        if not d.is_Integer or d < 0:
            D.remove(d)

    if not D:
        return None  # 'f(n)' is *not* Gosper-summable

    d = max(D)

    coeffs = symbols("c:%s" % (d + 1), cls=Dummy)
    domain = A.get_domain().inject(*coeffs)

    x = Poly(coeffs, n, domain=domain)
    H = A * x.shift(1) - B * x - C

    solution = solve(H.coeffs(), coeffs)

    if solution is None:
        return None  # 'f(n)' is *not* Gosper-summable

    x = x.as_expr().subs(solution)

    for coeff in coeffs:
        if coeff not in solution:
            x = x.subs(coeff, 0)

    if x is S.Zero:
        return None  # 'f(n)' is *not* Gosper-summable
    else:
        return B.as_expr() * x / C.as_expr()

# 需要导入模块: from sympy.polys import Poly [as 别名]
# 或者: from sympy.polys.Poly import subs [as 别名]
    def _ratsimpmodprime(a, b, allsol, N=0, D=0):
        r"""
        Computes a rational simplification of ``a/b`` which minimizes
        the sum of the total degrees of the numerator and the denominator.

        The algorithm proceeds by looking at ``a * d - b * c`` modulo
        the ideal generated by ``G`` for some ``c`` and ``d`` with degree
        less than ``a`` and ``b`` respectively.
        The coefficients of ``c`` and ``d`` are indeterminates and thus
        the coefficients of the normalform of ``a * d - b * c`` are
        linear polynomials in these indeterminates.
        If these linear polynomials, considered as system of
        equations, have a nontrivial solution, then `\frac{a}{b}
        \equiv \frac{c}{d}` modulo the ideal generated by ``G``. So,
        by construction, the degree of ``c`` and ``d`` is less than
        the degree of ``a`` and ``b``, so a simpler representation
        has been found.
        After a simpler representation has been found, the algorithm
        tries to reduce the degree of the numerator and denominator
        and returns the result afterwards.

        As an extension, if quick=False, we look at all possible degrees such
        that the total degree is less than *or equal to* the best current
        solution. We retain a list of all solutions of minimal degree, and try
        to find the best one at the end.
        """
        c, d = a, b
        steps = 0

        maxdeg = a.total_degree() + b.total_degree()
        if quick:
            bound = maxdeg - 1
        else:
            bound = maxdeg
        while N + D <= bound:
            if (N, D) in tested:
                break
            tested.add((N, D))

            M1 = staircase(N)
            M2 = staircase(D)
            debug('%s / %s: %s, %s' % (N, D, M1, M2))

            Cs = symbols("c:%d" % len(M1), cls=Dummy)
            Ds = symbols("d:%d" % len(M2), cls=Dummy)
            ng = Cs + Ds

            c_hat = Poly(
                sum([Cs[i] * M1[i] for i in range(len(M1))]), opt.gens + ng)
            d_hat = Poly(
                sum([Ds[i] * M2[i] for i in range(len(M2))]), opt.gens + ng)

            r = reduced(a * d_hat - b * c_hat, G, opt.gens + ng,
                        order=opt.order, polys=True)[1]

            S = Poly(r, gens=opt.gens).coeffs()
            sol = solve(S, Cs + Ds, particular=True, quick=True)

            if sol and not all([s == 0 for s in sol.values()]):
                c = c_hat.subs(sol)
                d = d_hat.subs(sol)

                # The "free" variables occurring before as parameters
                # might still be in the substituted c, d, so set them
                # to the value chosen before:
                c = c.subs(dict(list(zip(Cs + Ds, [1] * (len(Cs) + len(Ds))))))
                d = d.subs(dict(list(zip(Cs + Ds, [1] * (len(Cs) + len(Ds))))))

                c = Poly(c, opt.gens)
                d = Poly(d, opt.gens)
                if d == 0:
                    raise ValueError('Ideal not prime?')

                allsol.append((c_hat, d_hat, S, Cs + Ds))
                if N + D != maxdeg:
                    allsol = [allsol[-1]]

                break

            steps += 1
            N += 1
            D += 1

        if steps > 0:
            c, d, allsol = _ratsimpmodprime(c, d, allsol, N, D - steps)
            c, d, allsol = _ratsimpmodprime(c, d, allsol, N - steps, D)

        return c, d, allsol