Python sympy 模块，var() 实例源码

def constrained_by(self):
        """
        returns a list of parameters that constrain this parameter
        """
        if self._is_constraint is None:
            return []
        params = []
        for var in self.is_constraint._vars:
            param = var.get_parameter()
            if param.uniqueid != self.uniqueid:
                params.append(param)
        return params

    #~ @property
    #~ def in_constraints(self):
        #~ """
        #~ returns a list the labels of the constraints in which this parameter constrains another
        #~
        #~ you can then access all of the parameters of a given constraint via bundle.get_constraint(constraint)
        #~ """
        #~ return [param.uniquetwig for param in self.in_constraints_expressions]

def related_to(self):
        """
        returns a list of all parameters that are either constrained by or constrain this parameter
        """
        params = []
        constraints = self.in_constraints
        if self.is_constraint is not None:
            constraints.append(self.is_constraint)

        for constraint in constraints:
            for var in constraint._vars:
                param = var.get_parameter()
                if param not in params and param.uniqueid != self.uniqueid:
                    params.append(param)

        return params


#    def sin(self):
#        return np.sin(self.get_value(unit=u.rad))

def test_sylvester():
    x = var('x')

    assert sylvester(x**3 -7, 0, x) == sylvester(x**3 -7, 0, x, 1) == Matrix([[0]])
    assert sylvester(0, x**3 -7, x) == sylvester(0, x**3 -7, x, 1) == Matrix([[0]])
    assert sylvester(x**3 -7, 0, x, 2) == Matrix([[0]])
    assert sylvester(0, x**3 -7, x, 2) == Matrix([[0]])

    assert sylvester(x**3 -7, 7, x).det() == sylvester(x**3 -7, 7, x, 1).det() == 343
    assert sylvester(7, x**3 -7, x).det() == sylvester(7, x**3 -7, x, 1).det() == 343
    assert sylvester(x**3 -7, 7, x, 2).det() == -343
    assert sylvester(7, x**3 -7, x, 2).det() == 343

    assert sylvester(3, 7, x).det() == sylvester(3, 7, x, 1).det() == sylvester(3, 7, x, 2).det() == 1

    assert sylvester(3, 0, x).det() == sylvester(3, 0, x, 1).det() == sylvester(3, 0, x, 2).det() == 1

    assert sylvester(x - 3, x - 8, x) == sylvester(x - 3, x - 8, x, 1) == sylvester(x - 3, x - 8, x, 2) == Matrix([[1, -3], [1, -8]])

    assert sylvester(x**3 - 7*x + 7, 3*x**2 - 7, x) == sylvester(x**3 - 7*x + 7, 3*x**2 - 7, x, 1) == Matrix([[1, 0, -7,  7,  0], [0, 1,  0, -7,  7], [3, 0, -7,  0,  0], [0, 3,  0, -7,  0], [0, 0,  3,  0, -7]])

    assert sylvester(x**3 - 7*x + 7, 3*x**2 - 7, x, 2) == Matrix([
[1, 0, -7,  7,  0,  0], [0, 3,  0, -7,  0,  0], [0, 1,  0, -7,  7,  0], [0, 0,  3,  0, -7,  0], [0, 0,  1,  0, -7,  7], [0, 0,  0,  3,  0, -7]])

def test_print_as_full():
    xvar = sympy.var('xvar')
    yvar = sympy.var('yvar')
    res1 = '''
# m
# m_num=2
m[0, 0] += 1
m[0, 1] += 2*xvar + (yvar*yvar*yvar)
'''
    res1 = res1.strip()

    m = sympy.Matrix([[1, 2*xvar + yvar**3]])
    assert print_as_full(m, 'm') == res1

    res2 = '''
# subs
# yvar = xvar
# m
# m_num=2
m[0, 0] += 1
m[0, 1] += (yvar*yvar*yvar) + 2*yvar
'''
    res2 = res2.strip()
    assert print_as_full(m, 'm', subs={xvar: yvar}) == res2

def _make_z1():
    var("z1")

# make z2 with call-depth = 2

def __make_z2():
    var("z2")

def test_var():
    var("a")
    assert a == Symbol("a")

    var("b bb cc zz _x")
    assert b == Symbol("b")
    assert bb == Symbol("bb")
    assert cc == Symbol("cc")
    assert zz == Symbol("zz")
    assert _x == Symbol("_x")

    v = var(['d', 'e', 'fg'])
    assert d == Symbol('d')
    assert e == Symbol('e')
    assert fg == Symbol('fg')

    # check return value
    assert v == [d, e, fg]

    # see if var() really injects into global namespace
    raises(NameError, lambda: z1)
    _make_z1()
    assert z1 == Symbol("z1")

    raises(NameError, lambda: z2)
    _make_z2()
    assert z2 == Symbol("z2")

def test_var_return():
    raises(ValueError, lambda: var(''))
    v2 = var('q')
    v3 = var('q p')

    assert v2 == Symbol('q')
    assert v3 == (Symbol('q'), Symbol('p'))

def test_var_accepts_comma():
    v1 = var('x y z')
    v2 = var('x,y,z')
    v3 = var('x,y z')

    assert v1 == v2
    assert v1 == v3

def test_var_cls():
    f = var('f', cls=Function)

    assert isinstance(f, FunctionClass)

    g, h = var('g,h', cls=Function)

    assert isinstance(g, FunctionClass)
    assert isinstance(h, FunctionClass)

def check_expression(expr, var_symbols):
    """Does eval(expr) both in Sage and SymPy and does other checks."""

    # evaluate the expression in the context of Sage:
    if var_symbols:
        sage.var(var_symbols)
    a = globals().copy()
    # safety checks...
    assert not "sin" in a
    a.update(sage.__dict__)
    assert "sin" in a
    e_sage = eval(expr, a)
    assert not isinstance(e_sage, sympy.Basic)

    # evaluate the expression in the context of SymPy:
    if var_symbols:
        sympy.var(var_symbols)
    b = globals().copy()
    assert not "sin" in b
    b.update(sympy.__dict__)
    assert "sin" in b
    b.update(sympy.__dict__)
    e_sympy = eval(expr, b)
    assert isinstance(e_sympy, sympy.Basic)

    # Do the actual checks:
    assert sympy.S(e_sage) == e_sympy
    assert e_sage == sage.SR(e_sympy)

def test_issue_4023():
    sage.var("a x")
    log = sage.log
    i = sympy.integrate(log(x)/a, (x, a, a + 1))
    i2 = sympy.simplify(i)
    s = sage.SR(i2)
    assert s == (a*log(1 + a) - a*log(a) + log(1 + a) - 1)/a

# This string contains Sage doctests, that execute all the functions above.
# When you add a new function, please add it here as well.

def energy_corrections(perturbation, n, a=10, mass=0.5):
    """
    Calculating first two order corrections due to perturbation theory and
    returns tuple where zero element is unperturbated energy, and two second
    is corrections

    ``n``
        the "nodal" quantum number. Corresponds to the number of nodes in the
        wavefunction.  n >= 0
    ``a``
        width of the well. a > 0
    ``mass``
        mass.

    """
    x, _a = var("x _a")

    Vnm = lambda n, m, a: Integral(X_n(n, a, x) * X_n(m, a, x)
        * perturbation.subs({_a: a}), (x, 0, a)).n()

    # As we know from theory for V0*r/a we will just V(n, n-1) and V(n, n+1)
    #   wouldn't equals zero

    return (E_n(n, a, mass).evalf(),

            Vnm(n, n, a).evalf(),

            (Vnm(n, n - 1, a)**2/(E_n(n, a, mass) - E_n(n - 1, a, mass))
           + Vnm(n, n + 1, a)**2/(E_n(n, a, mass) - E_n(n + 1, a, mass))).evalf())

def main():
    print()
    print("Applying perturbation theory to calculate the ground state energy")
    print("of the infinite 1D box of width ``a`` with a perturbation")
    print("which is linear in ``x``, up to second order in perturbation.")
    print()

    x, _a = var("x _a")
    perturbation = .1 * x / _a

    E1 = energy_corrections(perturbation, 1)
    print("Energy for first term (n=1):")
    print("E_1^{(0)} = ", E1[0])
    print("E_1^{(1)} = ", E1[1])
    print("E_1^{(2)} = ", E1[2])
    print()

    E2 = energy_corrections(perturbation, 2)
    print("Energy for second term (n=2):")
    print("E_2^{(0)} = ", E2[0])
    print("E_2^{(1)} = ", E2[1])
    print("E_2^{(2)} = ", E2[2])
    print()

    E3 = energy_corrections(perturbation, 3)
    print("Energy for third term (n=3):")
    print("E_3^{(0)} = ", E3[0])
    print("E_3^{(1)} = ", E3[1])
    print("E_3^{(2)} = ", E3[2])
    print()

def constrains(self):
        """
        returns a list of parameters that are constrained by this parameter
        """
        params = []
        for constraint in self.in_constraints:
            for var in constraint._vars:
                param = var.get_parameter()
                if param.component == constraint.component and param.qualifier == constraint.qualifier:
                    if param not in params and param.uniqueid != self.uniqueid:
                        params.append(param)
        return params

def vars(self):
        """
        return all the variables in a PS
        """

        return ParameterSet([var.get_parameter() for var in self._vars])

def _get_var(self, param=None, **kwargs):
        if not isinstance(param, Parameter):
            if isinstance(param, str) and 'twig' not in kwargs.keys():
                kwargs['twig'] = param

            param = self.get_parameter(**kwargs)

        varids = [var.unique_label for var in self._vars]
        if param.uniqueid not in varids:
            raise KeyError("{} was not found in expression".format(param.uniquetwig))
        return self._vars[varids.index(param.uniqueid)]

def set_value(self, value, **kwargs):
        """

        kwargs are passed on to filter
        """
        _orig_value = deepcopy(self.get_value())

        if self._bundle is None:
            raise ValueError("ConstraintParameters must be attached from the bundle, and cannot be standalone")
        value = str(value) # <-- in case unicode
        # if the user wants to see the expression, we'll replace all
        # var.safe_label with var.curly_label
        self._value, self._vars = self._parse_expr(value)
        #~ print "***", self.uniquetwig, self.uniqueid
        self._add_history(redo_func='set_value', redo_kwargs={'value': value, 'uniqueid': self.uniqueid}, undo_func='set_value', undo_kwargs={'value': _orig_value, 'uniqueid': self.uniqueid})

def _update_bookkeeping(self):
        # do bookkeeping on parameters
        self._remove_bookkeeping()
        for var in self._vars:
            param = var.get_parameter()
            if param.qualifier == self.qualifier and param.component == self.component:
                # then this is the currently constrained parameter
                param._is_constraint = self.uniqueid
                if param.uniqueid in param._in_constraints:
                    param._in_constraints.remove(self.uniqueid)
            else:
                # then this is a constraining parameter
                if self.uniqueid not in param._in_constraints:
                    param._in_constraints.append(self.uniqueid)

def _remove_bookkeeping(self):
        for var in self._vars:
            param = var.get_parameter()
            if param._is_constraint == self.uniqueid:
                param._is_constraint = None
            if self.uniqueid in param._in_constraints:
                param._in_constraints.remove(self.uniqueid)

def get_value(self):
        """
        """
        # for access to the sympy-safe expr, just use self._expr
        expr = self._value
        for var in self._vars:
            # update to current unique twig
            var.update_user_label()  # update curly label
            #~ print "***", expr, var.safe_label, var.curly_label
            expr = expr.replace(str(var.safe_label), str(var.curly_label))

        return expr

def PSC(F, G, x):
    """PSC(F, G, x) returns a set with the non-zero principal subresultant
    coefficients (psc) of the two polynomials F and G with respect to the
    variable x.

    If the degree of the polynomial F is strictly less than the degree of
    the polynomial G, then F and G are interchanged.

    Extended psc beyond the n-th are not considered, where n is the minimum
    of the degrees of F and G with respect to x.

    >>> PSC(0, 0, var('x'))
    set()
    >>> PSC(poly('2*x'), poly('3*y'), var('x'))
    {3*y}
    >>> PSC(poly('2*x'), poly('3*y + 5*x**2'), var('x')) == {12*var('y'), 2}
    True
    >>> PSC(poly('x**3'), poly('x**3 + x'), var('x'))
    {1}
    """
    subs = subresultants(F, G, x)
    s = set()
    i = len(subs) - 1
    if i < 0:
        return s
    currDeg = degree(subs[i], x)
    while i > 0:
        nextDeg = degree(subs[i-1], x)
        s.add( LC(subs[i], x)**(nextDeg-currDeg) )
        currDeg = nextDeg
        i -= 1
    return s


# Execute doctest when run from the command line

def _make_z1():
    var("z1")

# make z2 with call-depth = 2

def __make_z2():
    var("z2")

def test_var():
    var("a")
    assert a == Symbol("a")

    var("b bb cc zz _x")
    assert b == Symbol("b")
    assert bb == Symbol("bb")
    assert cc == Symbol("cc")
    assert zz == Symbol("zz")
    assert _x == Symbol("_x")

    v = var(['d', 'e', 'fg'])
    assert d == Symbol('d')
    assert e == Symbol('e')
    assert fg == Symbol('fg')

    # check return value
    assert v == [d, e, fg]

    # see if var() really injects into global namespace
    raises(NameError, lambda: z1)
    _make_z1()
    assert z1 == Symbol("z1")

    raises(NameError, lambda: z2)
    _make_z2()
    assert z2 == Symbol("z2")

def test_var_return():
    raises(ValueError, lambda: var(''))
    v2 = var('q')
    v3 = var('q p')

    assert v2 == Symbol('q')
    assert v3 == (Symbol('q'), Symbol('p'))

def test_var_accepts_comma():
    v1 = var('x y z')
    v2 = var('x,y,z')
    v3 = var('x,y z')

    assert v1 == v2
    assert v1 == v3

def test_var_cls():
    f = var('f', cls=Function)

    assert isinstance(f, FunctionClass)

    g, h = var('g,h', cls=Function)

    assert isinstance(g, FunctionClass)
    assert isinstance(h, FunctionClass)

def check_expression(expr, var_symbols, only_from_sympy=False):
    """
    Does eval(expr) both in Sage and SymPy and does other checks.
    """

    # evaluate the expression in the context of Sage:
    if var_symbols:
        sage.var(var_symbols)
    a = globals().copy()
    # safety checks...
    a.update(sage.__dict__)
    assert "sin" in a
    is_different = False
    try:
        e_sage = eval(expr, a)
        assert not isinstance(e_sage, sympy.Basic)
    except (NameError, TypeError):
        is_different = True
        pass

    # evaluate the expression in the context of SymPy:
    if var_symbols:
        sympy_vars = sympy.var(var_symbols)
    b = globals().copy()
    b.update(sympy.__dict__)
    assert "sin" in b
    b.update(sympy.__dict__)
    e_sympy = eval(expr, b)
    assert isinstance(e_sympy, sympy.Basic)

    # Sympy func may have specific _sage_ method
    if is_different:
        _sage_method = getattr(e_sympy.func, "_sage_")
        e_sage = _sage_method(sympy.S(e_sympy))

    # Do the actual checks:
    if not only_from_sympy:
        assert sympy.S(e_sage) == e_sympy
    assert e_sage == sage.SR(e_sympy)

def test_issue_4023():
    sage.var("a x")
    log = sage.log
    i = sympy.integrate(log(x)/a, (x, a, a + 1))
    i2 = sympy.simplify(i)
    s = sage.SR(i2)
    assert s == (a*log(1 + a) - a*log(a) + log(1 + a) - 1)/a

def test_undefined_function():
    f = sympy.Function('f')
    sf = sage.function('f')
    x = sympy.symbols('x')
    sx = sage.var('x')
    assert bool(sf(sx) == f(x)._sage_())
    #assert bool(f == sympy.sympify(sf))

# This string contains Sage doctests, that execute all the functions above.
# When you add a new function, please add it here as well.

def test_subresultants_bezout():
    x = var('x')

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    assert subresultants_bezout(p, q, x) == subresultants(p, q, x)
    assert subresultants_bezout(p, q, x)[-1] == sylvester(p, q, x).det()
    assert subresultants_bezout(p, q, x) != euclid_amv(p, q, x)
    amv_factors = [1, 1, -1, 1, -1, 1]
    assert subresultants_bezout(p, q, x) == [i*j for i, j in zip(amv_factors, modified_subresultants_amv(p, q, x))]

    p = x**3 - 7*x + 7
    q = 3*x**2 - 7
    assert subresultants_bezout(p, q, x) == euclid_amv(p, q, x)

def test_modified_subresultants_bezout():
    x = var('x')

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    amv_factors = [1, 1, -1, 1, -1, 1]
    assert modified_subresultants_bezout(p, q, x) == [i*j for i, j in zip(amv_factors, subresultants_amv(p, q, x))]
    assert modified_subresultants_bezout(p, q, x)[-1] != sylvester(p + x**8, q, x).det()
    assert modified_subresultants_bezout(p, q, x) != sturm_amv(p, q, x)

    p = x**3 - 7*x + 7
    q = 3*x**2 - 7
    assert modified_subresultants_bezout(p, q, x) == sturm_amv(p, q, x)
    assert modified_subresultants_bezout(-p, q, x) != sturm_amv(-p, q, x)

def test_sturm_pg():
    x = var('x')

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    assert sturm_pg(p, q, x)[-1] != sylvester(p, q, x, 2).det()
    sam_factors = [1, 1, -1, -1, 1, 1]
    assert sturm_pg(p, q, x) == [i*j for i,j in zip(sam_factors, euclid_pg(p, q, x))]

    p = -9*x**5 - 5*x**3 - 9
    q = -45*x**4 - 15*x**2
    assert sturm_pg(p, q, x, 1)[-1] == sylvester(p, q, x, 1).det()
    assert sturm_pg(p, q, x)[-1] != sylvester(p, q, x, 2).det()
    assert sturm_pg(-p, q, x)[-1] == sylvester(-p, q, x, 2).det()
    assert sturm_pg(-p, q, x) == modified_subresultants_pg(-p, q, x)

def test_sturm_q():
    x = var('x')

    p = x**3 - 7*x + 7
    q = 3*x**2 - 7
    assert sturm_q(p, q, x) == sturm(p)
    assert sturm_q(-p, -q, x) != sturm(-p)

def test_sturm_amv():
    x = var('x')

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    assert sturm_amv(p, q, x)[-1] != sylvester(p, q, x, 2).det()
    sam_factors = [1, 1, -1, -1, 1, 1]
    assert sturm_amv(p, q, x) == [i*j for i,j in zip(sam_factors, euclid_amv(p, q, x))]

    p = -9*x**5 - 5*x**3 - 9
    q = -45*x**4 - 15*x**2
    assert sturm_amv(p, q, x, 1)[-1] == sylvester(p, q, x, 1).det()
    assert sturm_amv(p, q, x)[-1] != sylvester(p, q, x, 2).det()
    assert sturm_amv(-p, q, x)[-1] == sylvester(-p, q, x, 2).det()
    assert sturm_pg(-p, q, x) == modified_subresultants_pg(-p, q, x)

def test_euclid_q():
    x = var('x')

    p = x**3 - 7*x + 7
    q = 3*x**2 - 7
    assert euclid_q(p, q, x)[-1] == -sturm(p)[-1]

def test_euclid_amv():
    x = var('x')

    p = x**3 - 7*x + 7
    q = 3*x**2 - 7
    assert euclid_amv(p, q, x)[-1] == sylvester(p, q, x).det()
    assert euclid_amv(p, q, x) == subresultants_amv(p, q, x)

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    assert euclid_amv(p, q, x)[-1] != sylvester(p, q, x, 2).det()
    sam_factors = [1, 1, -1, -1, 1, 1]
    assert euclid_amv(p, q, x) == [i*j for i,j in zip(sam_factors, sturm_amv(p, q, x))]

def test_modified_subresultants_pg():
    x = var('x')

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    amv_factors = [1, 1, -1, 1, -1, 1]
    assert modified_subresultants_pg(p, q, x) == [i*j for i, j in zip(amv_factors, subresultants_pg(p, q, x))]
    assert modified_subresultants_pg(p, q, x)[-1] != sylvester(p + x**8, q, x).det()
    assert modified_subresultants_pg(p, q, x) != sturm_pg(p, q, x)

    p = x**3 - 7*x + 7
    q = 3*x**2 - 7
    assert modified_subresultants_pg(p, q, x) == sturm_pg(p, q, x)
    assert modified_subresultants_pg(-p, q, x) != sturm_pg(-p, q, x)

def test_subresultants_pg():
    x = var('x')

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    assert subresultants_pg(p, q, x) == subresultants(p, q, x)
    assert subresultants_pg(p, q, x)[-1] == sylvester(p, q, x).det()
    assert subresultants_pg(p, q, x) != euclid_pg(p, q, x)
    amv_factors = [1, 1, -1, 1, -1, 1]
    assert subresultants_pg(p, q, x) == [i*j for i, j in zip(amv_factors, modified_subresultants_amv(p, q, x))]

    p = x**3 - 7*x + 7
    q = 3*x**2 - 7
    assert subresultants_pg(p, q, x) == euclid_pg(p, q, x)

def test_subresultants_amv_q():
    x = var('x')

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    assert subresultants_amv_q(p, q, x) == subresultants(p, q, x)
    assert subresultants_amv_q(p, q, x)[-1] == sylvester(p, q, x).det()
    assert subresultants_amv_q(p, q, x) != euclid_amv(p, q, x)
    amv_factors = [1, 1, -1, 1, -1, 1]
    assert subresultants_amv_q(p, q, x) == [i*j for i, j in zip(amv_factors, modified_subresultants_amv(p, q, x))]

    p = x**3 - 7*x + 7
    q = 3*x**2 - 7
    assert subresultants_amv(p, q, x) == euclid_amv(p, q, x)

def test_quo_z():
    x = var('x')

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    assert quo_z(p, -q, x) != pquo(p, -q, x)

def test_subresultants_amv():
    x = var('x')

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    assert subresultants_amv(p, q, x) == subresultants(p, q, x)
    assert subresultants_amv(p, q, x)[-1] == sylvester(p, q, x).det()
    assert subresultants_amv(p, q, x) != euclid_amv(p, q, x)
    amv_factors = [1, 1, -1, 1, -1, 1]
    assert subresultants_amv(p, q, x) == [i*j for i, j in zip(amv_factors, modified_subresultants_amv(p, q, x))]

    p = x**3 - 7*x + 7
    q = 3*x**2 - 7
    assert subresultants_amv(p, q, x) == euclid_amv(p, q, x)

def test_modified_subresultants_amv():
    x = var('x')

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    amv_factors = [1, 1, -1, 1, -1, 1]
    assert modified_subresultants_amv(p, q, x) == [i*j for i, j in zip(amv_factors, subresultants_amv(p, q, x))]
    assert modified_subresultants_amv(p, q, x)[-1] != sylvester(p + x**8, q, x).det()
    assert modified_subresultants_amv(p, q, x) != sturm_amv(p, q, x)

    p = x**3 - 7*x + 7
    q = 3*x**2 - 7
    assert modified_subresultants_amv(p, q, x) == sturm_amv(p, q, x)
    assert modified_subresultants_amv(-p, q, x) != sturm_amv(-p, q, x)

def test_subresultants_rem():
    x = var('x')

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    assert subresultants_rem(p, q, x) == subresultants(p, q, x)
    assert subresultants_rem(p, q, x)[-1] == sylvester(p, q, x).det()
    assert subresultants_rem(p, q, x) != euclid_amv(p, q, x)
    amv_factors = [1, 1, -1, 1, -1, 1]
    assert subresultants_rem(p, q, x) == [i*j for i, j in zip(amv_factors, modified_subresultants_amv(p, q, x))]

    p = x**3 - 7*x + 7
    q = 3*x**2 - 7
    assert subresultants_rem(p, q, x) == euclid_amv(p, q, x)

def test_subresultants_vv():
    x = var('x')

    p = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5
    q = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21
    assert subresultants_vv(p, q, x) == subresultants(p, q, x)
    assert subresultants_vv(p, q, x)[-1] == sylvester(p, q, x).det()
    assert subresultants_vv(p, q, x) != euclid_amv(p, q, x)
    amv_factors = [1, 1, -1, 1, -1, 1]
    assert subresultants_vv(p, q, x) == [i*j for i, j in zip(amv_factors, modified_subresultants_amv(p, q, x))]

    p = x**3 - 7*x + 7
    q = 3*x**2 - 7
    assert subresultants_vv(p, q, x) == euclid_amv(p, q, x)

def test_print_as_sparse():
    xvar = sympy.var('xvar')
    yvar = sympy.var('yvar')
    m = sympy.Matrix([[1, 2*xvar + yvar**3]])

    res1 = """
# test_m
# test_m_num=2
c += 1
test_mr[c] = row+0
test_mc[c] = col+0
test_mv[c] += 1
c += 1
test_mr[c] = row+0
test_mc[c] = col+1
test_mv[c] += 2*xvar + (yvar*yvar*yvar)
"""
    res1 = res1.strip()
    print_as_sparse(m, 'test_m', print_file=False) == res1
    print_as_sparse(m, 'test_m', print_file=False, full_symmetric=True) == res1

    res2 = """
# test_m
# test_m_num=2
c += 1
test_mr[c] = row+0
test_mc[c] = col+0
test_mv[c] += 1
c += 1
test_mr[c] = row+0
test_mc[c] = col+1
test_mv[c] += (yvar*yvar*yvar) + 2*yvar
"""
    res2 = res2.strip()
    print_as_sparse(m, 'test_m', print_file=False, subs={xvar: yvar}) == res2

def test_print_as_array():
    xvar = sympy.var('xvar')
    yvar = sympy.var('yvar')
    m = sympy.Matrix([[1, 2*xvar + yvar**3, + 2*xvar*(yvar+1)]])
    res1 = '''
# m
# m_num=3
m[pos+0] += 1
m[pos+1] += 2*xvar + (yvar*yvar*yvar)
m[pos+2] += 2*xvar*(yvar + 1)
    '''
    res1 = res1.strip()
    assert print_as_array(m, 'm') == res1

    res2 = '''
# cdefs
cdef double x0
# subs
x0 = 2*xvar
# m
# m_num=3
m[pos+0] += 1
m[pos+1] += x0 + (yvar*yvar*yvar)
m[pos+2] += x0*(yvar + 1)
'''
    res2 = res2.strip()
    assert print_as_array(m, 'm', use_cse=True) == res2