def test_eval_trace():
up = JzKet(S(1)/2, S(1)/2)
down = JzKet(S(1)/2, -S(1)/2)
d = Density((up, 0.5), (down, 0.5))
t = Tr(d)
assert t.doit() == 1
#test dummy time dependent states
class TestTimeDepKet(TimeDepKet):
def _eval_trace(self, bra, **options):
return 1
x, t = symbols('x t')
k1 = TestTimeDepKet(0, 0.5)
k2 = TestTimeDepKet(0, 1)
d = Density([k1, 0.5], [k2, 0.5])
assert d.doit() == (0.5 * OuterProduct(k1, k1.dual) +
0.5 * OuterProduct(k2, k2.dual))
t = Tr(d)
assert t.doit() == 1
python类S的实例源码
def test_wavefunction():
a = 1/Z
R = {
(1, 0): 2*sqrt(1/a**3) * exp(-r/a),
(2, 0): sqrt(1/(2*a**3)) * exp(-r/(2*a)) * (1 - r/(2*a)),
(2, 1): S(1)/2 * sqrt(1/(6*a**3)) * exp(-r/(2*a)) * r/a,
(3, 0): S(2)/3 * sqrt(1/(3*a**3)) * exp(-r/(3*a)) *
(1 - 2*r/(3*a) + S(2)/27 * (r/a)**2),
(3, 1): S(4)/27 * sqrt(2/(3*a**3)) * exp(-r/(3*a)) *
(1 - r/(6*a)) * r/a,
(3, 2): S(2)/81 * sqrt(2/(15*a**3)) * exp(-r/(3*a)) * (r/a)**2,
(4, 0): S(1)/4 * sqrt(1/a**3) * exp(-r/(4*a)) *
(1 - 3*r/(4*a) + S(1)/8 * (r/a)**2 - S(1)/192 * (r/a)**3),
(4, 1): S(1)/16 * sqrt(5/(3*a**3)) * exp(-r/(4*a)) *
(1 - r/(4*a) + S(1)/80 * (r/a)**2) * (r/a),
(4, 2): S(1)/64 * sqrt(1/(5*a**3)) * exp(-r/(4*a)) *
(1 - r/(12*a)) * (r/a)**2,
(4, 3): S(1)/768 * sqrt(1/(35*a**3)) * exp(-r/(4*a)) * (r/a)**3,
}
for n, l in R:
assert simplify(R_nl(n, l, r, Z) - R[(n, l)]) == 0
def test_transpose():
Sq = MatrixSymbol('Sq', n, n)
assert transpose(A) == Transpose(A)
assert Transpose(A).shape == (m, n)
assert Transpose(A*B).shape == (l, n)
assert transpose(Transpose(A)) == A
assert isinstance(Transpose(Transpose(A)), Transpose)
assert adjoint(Transpose(A)) == Adjoint(Transpose(A))
assert conjugate(Transpose(A)) == Adjoint(A)
assert Transpose(eye(3)).doit() == eye(3)
assert Transpose(S(5)).doit() == S(5)
assert Transpose(Matrix([[1, 2], [3, 4]])).doit() == Matrix([[1, 3], [2, 4]])
assert transpose(trace(Sq)) == trace(Sq)
assert trace(Transpose(Sq)) == trace(Sq)
assert Transpose(Sq)[0, 1] == Sq[1, 0]
assert Transpose(A*B).doit() == Transpose(B) * Transpose(A)
def test_spde():
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(t**2 + 1, t)]})
raises(NonElementaryIntegralException, lambda: spde(Poly(t, t), Poly((t - 1)*(t**2 + 1), t), Poly(1, t), 0, DE))
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(t, t)]})
assert spde(Poly(t**2 + x*t*2 + x**2, t), Poly(t**2/x**2 + (2/x - 1)*t, t),
Poly(t**2/x**2 + (2/x - 1)*t, t), 0, DE) == \
(Poly(0, t), Poly(0, t), 0, Poly(0, t), Poly(1, t))
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(t0/x**2, t0), Poly(1/x, t)]})
assert spde(Poly(t**2, t), Poly(-t**2/x**2 - 1/x, t),
Poly((2*x - 1)*t**4 + (t0 + x)/x*t**3 - (t0 + 4*x**2)/(2*x)*t**2 + x*t, t), 3, DE) == \
(Poly(0, t), Poly(0, t), 0, Poly(0, t),
Poly(t0*t**2/2 + x**2*t**2 - x**2*t, t))
DE = DifferentialExtension(extension={'D': [Poly(1, x)]})
assert spde(Poly(x**2 + x + 1, x), Poly(-2*x - 1, x), Poly(x**5/2 +
3*x**4/4 + x**3 - x**2 + 1, x), 4, DE) == \
(Poly(0, x), Poly(x/2 - S(1)/4, x), 2, Poly(x**2 + x + 1, x), Poly(5*x/4, x))
assert spde(Poly(x**2 + x + 1, x), Poly(-2*x - 1, x), Poly(x**5/2 +
3*x**4/4 + x**3 - x**2 + 1, x), n, DE) == \
(Poly(0, x), Poly(x/2 - S(1)/4, x), -2 + n, Poly(x**2 + x + 1, x), Poly(5*x/4, x))
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(1, t)]})
raises(NonElementaryIntegralException, lambda: spde(Poly((t - 1)*(t**2 + 1)**2, t), Poly((t - 1)*(t**2 + 1), t), Poly(1, t), 0, DE))
def test_solve_poly_rde_no_cancel():
# deg(b) large
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(1 + t**2, t)]})
assert solve_poly_rde(Poly(t**2 + 1, t), Poly(t**3 + (x + 1)*t**2 + t + x + 2, t),
oo, DE) == Poly(t + x, t)
# deg(b) small
DE = DifferentialExtension(extension={'D': [Poly(1, x)]})
assert solve_poly_rde(Poly(0, x), Poly(x/2 - S(1)/4, x), oo, DE) == \
Poly(x**2/4 - x/4, x)
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(t**2 + 1, t)]})
assert solve_poly_rde(Poly(2, t), Poly(t**2 + 2*t + 3, t), 1, DE) == \
Poly(t + 1, t, x)
# deg(b) == deg(D) - 1
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(t**2 + 1, t)]})
assert no_cancel_equal(Poly(1 - t, t),
Poly(t**3 + t**2 - 2*x*t - 2*x, t), oo, DE) == \
(Poly(t**2, t), 1, Poly((-2 - 2*x)*t - 2*x, t))
def test_is_deriv_k():
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(1/x, t1), Poly(1/(x + 1), t2)],
'L_K': [1, 2], 'E_K': [], 'L_args': [x, x + 1], 'E_args': []})
assert is_deriv_k(Poly(2*x**2 + 2*x, t2), Poly(1, t2), DE) == \
([(t1, 1), (t2, 1)], t1 + t2, 2)
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(1/x, t1), Poly(t2, t2)],
'L_K': [1], 'E_K': [2], 'L_args': [x], 'E_args': [x]})
assert is_deriv_k(Poly(x**2*t2**3, t2), Poly(1, t2), DE) == \
([(x, 3), (t1, 2)], 2*t1 + 3*x, 1)
# TODO: Add more tests, including ones with exponentials
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(2/x, t1)],
'L_K': [1], 'E_K': [], 'L_args': [x**2], 'E_args': []})
assert is_deriv_k(Poly(x, t1), Poly(1, t1), DE) == \
([(t1, S(1)/2)], t1/2, 1)
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(2/(1 + x), t0)],
'L_K': [1], 'E_K': [], 'L_args': [x**2 + 2*x + 1], 'E_args': []})
assert is_deriv_k(Poly(1 + x, t0), Poly(1, t0), DE) == \
([(t0, S(1)/2)], t0/2, 1)
def test_dup_clear_denoms():
assert dup_clear_denoms([], QQ, ZZ) == (ZZ(1), [])
assert dup_clear_denoms([QQ(1)], QQ, ZZ) == (ZZ(1), [QQ(1)])
assert dup_clear_denoms([QQ(7)], QQ, ZZ) == (ZZ(1), [QQ(7)])
assert dup_clear_denoms([QQ(7, 3)], QQ) == (ZZ(3), [QQ(7)])
assert dup_clear_denoms([QQ(7, 3)], QQ, ZZ) == (ZZ(3), [QQ(7)])
assert dup_clear_denoms(
[QQ(3), QQ(1), QQ(0)], QQ, ZZ) == (ZZ(1), [QQ(3), QQ(1), QQ(0)])
assert dup_clear_denoms(
[QQ(1), QQ(1, 2), QQ(0)], QQ, ZZ) == (ZZ(2), [QQ(2), QQ(1), QQ(0)])
assert dup_clear_denoms([QQ(3), QQ(
1), QQ(0)], QQ, ZZ, convert=True) == (ZZ(1), [ZZ(3), ZZ(1), ZZ(0)])
assert dup_clear_denoms([QQ(1), QQ(
1, 2), QQ(0)], QQ, ZZ, convert=True) == (ZZ(2), [ZZ(2), ZZ(1), ZZ(0)])
assert dup_clear_denoms(
[EX(S(3)/2), EX(S(9)/4)], EX) == (EX(4), [EX(6), EX(9)])
assert dup_clear_denoms([EX(7)], EX) == (EX(1), [EX(7)])
assert dup_clear_denoms([EX(sin(x)/x), EX(0)], EX) == (EX(x), [EX(sin(x)), EX(0)])
def test_syzygy():
R = QQ.old_poly_ring(x, y, z)
M = R.free_module(1).submodule([x*y], [y*z], [x*z])
S = R.free_module(3).submodule([0, x, -y], [z, -x, 0])
assert M.syzygy_module() == S
M2 = M / ([x*y*z],)
S2 = R.free_module(3).submodule([z, 0, 0], [0, x, 0], [0, 0, y])
assert M2.syzygy_module() == S2
F = R.free_module(3)
assert F.submodule(*F.basis()).syzygy_module() == F.submodule()
R2 = QQ.old_poly_ring(x, y, z) / [x*y*z]
M3 = R2.free_module(1).submodule([x*y], [y*z], [x*z])
S3 = R2.free_module(3).submodule([z, 0, 0], [0, x, 0], [0, 0, y])
assert M3.syzygy_module() == S3
def test_in_terms_of_generators():
R = QQ.old_poly_ring(x, order="ilex")
M = R.free_module(2).submodule([2*x, 0], [1, 2])
assert M.in_terms_of_generators(
[x, x]) == [R.convert(S(1)/4), R.convert(x/2)]
raises(ValueError, lambda: M.in_terms_of_generators([1, 0]))
M = R.free_module(2) / ([x, 0], [1, 1])
SM = M.submodule([1, x])
assert SM.in_terms_of_generators([2, 0]) == [R.convert(-2/(x - 1))]
R = QQ.old_poly_ring(x, y) / [x**2 - y**2]
M = R.free_module(2)
SM = M.submodule([x, 0], [0, y])
assert SM.in_terms_of_generators(
[x**2, x**2]) == [R.convert(x), R.convert(y)]
def test_special_printers():
class IntervalPrinter(LambdaPrinter):
"""Use ``lambda`` printer but print numbers as ``mpi`` intervals. """
def _print_Integer(self, expr):
return "mpi('%s')" % super(IntervalPrinter, self)._print_Integer(expr)
def _print_Rational(self, expr):
return "mpi('%s')" % super(IntervalPrinter, self)._print_Rational(expr)
def intervalrepr(expr):
return IntervalPrinter().doprint(expr)
expr = sympy.sqrt(sympy.sqrt(2) + sympy.sqrt(3)) + sympy.S(1)/2
func0 = lambdify((), expr, modules="mpmath", printer=intervalrepr)
func1 = lambdify((), expr, modules="mpmath", printer=IntervalPrinter)
func2 = lambdify((), expr, modules="mpmath", printer=IntervalPrinter())
mpi = type(mpmath.mpi(1, 2))
assert isinstance(func0(), mpi)
assert isinstance(func1(), mpi)
assert isinstance(func2(), mpi)
def test_limit_seq():
e = binomial(2*n, n) / Sum(binomial(2*k, k), (k, 1, n))
assert limit_seq(e) == S(3) / 4
assert limit_seq(e, m) == e
e = (5*n**3 + 3*n**2 + 4) / (3*n**3 + 4*n - 5)
assert limit_seq(e, n) == S(5) / 3
e = (harmonic(n) * Sum(harmonic(k), (k, 1, n))) / (n * harmonic(2*n)**2)
assert limit_seq(e, n) == 1
e = Sum(k**2 * Sum(2**m/m, (m, 1, k)), (k, 1, n)) / (2**n*n)
assert limit_seq(e, n) == 4
e = (Sum(binomial(3*k, k) * binomial(5*k, k), (k, 1, n)) /
(binomial(3*n, n) * binomial(5*n, n)))
assert limit_seq(e, n) == S(84375) / 83351
e = Sum(harmonic(k)**2/k, (k, 1, 2*n)) / harmonic(n)**3
assert limit_seq(e, n) == S(1) / 3
raises(ValueError, lambda: limit_seq(e * m))
def test_limit_seq_fail():
# improve Summation algorithm or add ad-hoc criteria
e = (harmonic(n)**3 * Sum(1/harmonic(k), (k, 1, n)) /
(n * Sum(harmonic(k)/k, (k, 1, n))))
assert limit_seq(e, n) == 2
# No unique dominant term
e = (Sum(2**k * binomial(2*k, k) / k**2, (k, 1, n)) /
(Sum(2**k/k*2, (k, 1, n)) * Sum(binomial(2*k, k), (k, 1, n))))
assert limit_seq(e, n) == S(3) / 7
# Simplifications of summations needs to be improved.
e = n**3*Sum(2**k/k**2, (k, 1, n))**2 / (2**n * Sum(2**k/k, (k, 1, n)))
assert limit_seq(e, n) == 2
e = (harmonic(n) * Sum(2**k/k, (k, 1, n)) /
(n * Sum(2**k*harmonic(k)/k**2, (k, 1, n))))
assert limit_seq(e, n) == 1
e = (Sum(2**k*factorial(k) / k**2, (k, 1, 2*n)) /
(Sum(4**k/k**2, (k, 1, n)) * Sum(factorial(k), (k, 1, 2*n))))
assert limit_seq(e, n) == S(3) / 16
def __init__(self, list_of_poly, parent):
# the parent ring for this operator
# must be an RecurrenceOperatorAlgebra object
self.parent = parent
# sequence of polynomials in n for each power of Sn
# represents the operator
# convert the expressions into ring elements using from_sympy
if isinstance(list_of_poly, list):
for i, j in enumerate(list_of_poly):
if isinstance(j, int):
list_of_poly[i] = self.parent.base.from_sympy(S(j))
elif not isinstance(j, self.parent.base.dtype):
list_of_poly[i] = self.parent.base.from_sympy(j)
self.listofpoly = list_of_poly
self.order = len(self.listofpoly) - 1
def test_differentiate_finite():
x, y = symbols('x y')
f = Function('f')
res0 = differentiate_finite(f(x, y) + exp(42), x, y)
xm, xp, ym, yp = [v + sign*S(1)/2 for v, sign in product([x, y], [-1, 1])]
ref0 = f(xm, ym) + f(xp, yp) - f(xm, yp) - f(xp, ym)
assert (res0 - ref0).simplify() == 0
g = Function('g')
res1 = differentiate_finite(f(x)*g(x) + 42, x)
ref1 = (-f(x - S(1)/2) + f(x + S(1)/2))*g(x) + \
(-g(x - S(1)/2) + g(x + S(1)/2))*f(x)
assert (res1 - ref1).simplify() == 0
res2 = differentiate_finite(f(x) + x**3 + 42, x, points=[x-1, x+1],
evaluate=False)
ref2 = (f(x + 1) + (x + 1)**3 - f(x - 1) - (x - 1)**3)/2
assert (res2 - ref2).simplify() == 0
def test_boson_states():
a = BosonOp("a")
# Fock states
n = 3
assert (BosonFockBra(0) * BosonFockKet(1)).doit() == 0
assert (BosonFockBra(1) * BosonFockKet(1)).doit() == 1
assert qapply(BosonFockBra(n) * Dagger(a)**n * BosonFockKet(0)) \
== sqrt(prod(range(1, n+1)))
# Coherent states
alpha1, alpha2 = 1.2, 4.3
assert (BosonCoherentBra(alpha1) * BosonCoherentKet(alpha1)).doit() == 1
assert (BosonCoherentBra(alpha2) * BosonCoherentKet(alpha2)).doit() == 1
assert abs((BosonCoherentBra(alpha1) * BosonCoherentKet(alpha2)).doit() -
exp(-S(1) / 2 * (alpha1 - alpha2) ** 2)) < 1e-12
assert qapply(a * BosonCoherentKet(alpha1)) == \
alpha1 * BosonCoherentKet(alpha1)
def test_eval_trace():
up = JzKet(S(1)/2, S(1)/2)
down = JzKet(S(1)/2, -S(1)/2)
d = Density((up, 0.5), (down, 0.5))
t = Tr(d)
assert t.doit() == 1
#test dummy time dependent states
class TestTimeDepKet(TimeDepKet):
def _eval_trace(self, bra, **options):
return 1
x, t = symbols('x t')
k1 = TestTimeDepKet(0, 0.5)
k2 = TestTimeDepKet(0, 1)
d = Density([k1, 0.5], [k2, 0.5])
assert d.doit() == (0.5 * OuterProduct(k1, k1.dual) +
0.5 * OuterProduct(k2, k2.dual))
t = Tr(d)
assert t.doit() == 1
def test_clebsch_gordan3():
j_1 = S(3)/2
j_2 = S(3)/2
m = S(3)
j = S(3)
m_1 = S(3)/2
m_2 = S(3)/2
assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1
j_1 = S(3)/2
j_2 = S(3)/2
m = S(2)
j = S(2)
m_1 = S(3)/2
m_2 = S(1)/2
assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1/sqrt(2)
j_1 = S(3)/2
j_2 = S(3)/2
m = S(2)
j = S(3)
m_1 = S(3)/2
m_2 = S(1)/2
assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1/sqrt(2)
def test_clebsch_gordan5():
j_1 = S(5)/2
j_2 = S(1)
m = S(7)/2
j = S(7)/2
m_1 = S(5)/2
m_2 = 1
assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1
j_1 = S(5)/2
j_2 = S(1)
m = S(5)/2
j = S(5)/2
m_1 = S(5)/2
m_2 = 0
assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == sqrt(5)/sqrt(7)
j_1 = S(5)/2
j_2 = S(1)
m = S(3)/2
j = S(3)/2
m_1 = S(1)/2
m_2 = 1
assert clebsch_gordan(j_1, j_2, j, m_1, m_2, m) == 1/sqrt(15)
def test_wavefunction():
a = 1/Z
R = {
(1, 0): 2*sqrt(1/a**3) * exp(-r/a),
(2, 0): sqrt(1/(2*a**3)) * exp(-r/(2*a)) * (1 - r/(2*a)),
(2, 1): S(1)/2 * sqrt(1/(6*a**3)) * exp(-r/(2*a)) * r/a,
(3, 0): S(2)/3 * sqrt(1/(3*a**3)) * exp(-r/(3*a)) *
(1 - 2*r/(3*a) + S(2)/27 * (r/a)**2),
(3, 1): S(4)/27 * sqrt(2/(3*a**3)) * exp(-r/(3*a)) *
(1 - r/(6*a)) * r/a,
(3, 2): S(2)/81 * sqrt(2/(15*a**3)) * exp(-r/(3*a)) * (r/a)**2,
(4, 0): S(1)/4 * sqrt(1/a**3) * exp(-r/(4*a)) *
(1 - 3*r/(4*a) + S(1)/8 * (r/a)**2 - S(1)/192 * (r/a)**3),
(4, 1): S(1)/16 * sqrt(5/(3*a**3)) * exp(-r/(4*a)) *
(1 - r/(4*a) + S(1)/80 * (r/a)**2) * (r/a),
(4, 2): S(1)/64 * sqrt(1/(5*a**3)) * exp(-r/(4*a)) *
(1 - r/(12*a)) * (r/a)**2,
(4, 3): S(1)/768 * sqrt(1/(35*a**3)) * exp(-r/(4*a)) * (r/a)**3,
}
for n, l in R:
assert simplify(R_nl(n, l, r, Z) - R[(n, l)]) == 0
def test_transpose():
Sq = MatrixSymbol('Sq', n, n)
assert transpose(A) == Transpose(A)
assert Transpose(A).shape == (m, n)
assert Transpose(A*B).shape == (l, n)
assert transpose(Transpose(A)) == A
assert isinstance(Transpose(Transpose(A)), Transpose)
assert adjoint(Transpose(A)) == Adjoint(Transpose(A))
assert conjugate(Transpose(A)) == Adjoint(A)
assert Transpose(eye(3)).doit() == eye(3)
assert Transpose(S(5)).doit() == S(5)
assert Transpose(Matrix([[1, 2], [3, 4]])).doit() == Matrix([[1, 3], [2, 4]])
assert transpose(trace(Sq)) == trace(Sq)
assert trace(Transpose(Sq)) == trace(Sq)
assert Transpose(Sq)[0, 1] == Sq[1, 0]
assert Transpose(A*B).doit() == Transpose(B) * Transpose(A)
def test_spde():
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(t**2 + 1, t)]})
raises(NonElementaryIntegralException, lambda: spde(Poly(t, t), Poly((t - 1)*(t**2 + 1), t), Poly(1, t), 0, DE))
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(t, t)]})
assert spde(Poly(t**2 + x*t*2 + x**2, t), Poly(t**2/x**2 + (2/x - 1)*t, t),
Poly(t**2/x**2 + (2/x - 1)*t, t), 0, DE) == \
(Poly(0, t), Poly(0, t), 0, Poly(0, t), Poly(1, t))
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(t0/x**2, t0), Poly(1/x, t)]})
assert spde(Poly(t**2, t), Poly(-t**2/x**2 - 1/x, t),
Poly((2*x - 1)*t**4 + (t0 + x)/x*t**3 - (t0 + 4*x**2)/(2*x)*t**2 + x*t, t), 3, DE) == \
(Poly(0, t), Poly(0, t), 0, Poly(0, t),
Poly(t0*t**2/2 + x**2*t**2 - x**2*t, t))
DE = DifferentialExtension(extension={'D': [Poly(1, x)]})
assert spde(Poly(x**2 + x + 1, x), Poly(-2*x - 1, x), Poly(x**5/2 +
3*x**4/4 + x**3 - x**2 + 1, x), 4, DE) == \
(Poly(0, x), Poly(x/2 - S(1)/4, x), 2, Poly(x**2 + x + 1, x), Poly(5*x/4, x))
assert spde(Poly(x**2 + x + 1, x), Poly(-2*x - 1, x), Poly(x**5/2 +
3*x**4/4 + x**3 - x**2 + 1, x), n, DE) == \
(Poly(0, x), Poly(x/2 - S(1)/4, x), -2 + n, Poly(x**2 + x + 1, x), Poly(5*x/4, x))
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(1, t)]})
raises(NonElementaryIntegralException, lambda: spde(Poly((t - 1)*(t**2 + 1)**2, t), Poly((t - 1)*(t**2 + 1), t), Poly(1, t), 0, DE))
DE = DifferentialExtension(extension={'D': [Poly(1, x)]})
assert spde(Poly(x**2 - x, x), Poly(1, x), Poly(9*x**4 - 10*x**3 + 2*x**2, x), 4, DE) == (Poly(0, x), Poly(0, x), 0, Poly(0, x), Poly(3*x**3 - 2*x**2, x))
assert spde(Poly(x**2 - x, x), Poly(x**2 - 5*x + 3, x), Poly(x**7 - x**6 - 2*x**4 + 3*x**3 - x**2, x), 5, DE) == \
(Poly(1, x), Poly(x + 1, x), 1, Poly(x**4 - x**3, x), Poly(x**3 - x**2, x))
def test_solve_poly_rde_no_cancel():
# deg(b) large
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(1 + t**2, t)]})
assert solve_poly_rde(Poly(t**2 + 1, t), Poly(t**3 + (x + 1)*t**2 + t + x + 2, t),
oo, DE) == Poly(t + x, t)
# deg(b) small
DE = DifferentialExtension(extension={'D': [Poly(1, x)]})
assert solve_poly_rde(Poly(0, x), Poly(x/2 - S(1)/4, x), oo, DE) == \
Poly(x**2/4 - x/4, x)
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(t**2 + 1, t)]})
assert solve_poly_rde(Poly(2, t), Poly(t**2 + 2*t + 3, t), 1, DE) == \
Poly(t + 1, t, x)
# deg(b) == deg(D) - 1
DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(t**2 + 1, t)]})
assert no_cancel_equal(Poly(1 - t, t),
Poly(t**3 + t**2 - 2*x*t - 2*x, t), oo, DE) == \
(Poly(t**2, t), 1, Poly((-2 - 2*x)*t - 2*x, t))
def test_syzygy():
R = QQ.old_poly_ring(x, y, z)
M = R.free_module(1).submodule([x*y], [y*z], [x*z])
S = R.free_module(3).submodule([0, x, -y], [z, -x, 0])
assert M.syzygy_module() == S
M2 = M / ([x*y*z],)
S2 = R.free_module(3).submodule([z, 0, 0], [0, x, 0], [0, 0, y])
assert M2.syzygy_module() == S2
F = R.free_module(3)
assert F.submodule(*F.basis()).syzygy_module() == F.submodule()
R2 = QQ.old_poly_ring(x, y, z) / [x*y*z]
M3 = R2.free_module(1).submodule([x*y], [y*z], [x*z])
S3 = R2.free_module(3).submodule([z, 0, 0], [0, x, 0], [0, 0, y])
assert M3.syzygy_module() == S3
def test_in_terms_of_generators():
R = QQ.old_poly_ring(x, order="ilex")
M = R.free_module(2).submodule([2*x, 0], [1, 2])
assert M.in_terms_of_generators(
[x, x]) == [R.convert(S(1)/4), R.convert(x/2)]
raises(ValueError, lambda: M.in_terms_of_generators([1, 0]))
M = R.free_module(2) / ([x, 0], [1, 1])
SM = M.submodule([1, x])
assert SM.in_terms_of_generators([2, 0]) == [R.convert(-2/(x - 1))]
R = QQ.old_poly_ring(x, y) / [x**2 - y**2]
M = R.free_module(2)
SM = M.submodule([x, 0], [0, y])
assert SM.in_terms_of_generators(
[x**2, x**2]) == [R.convert(x), R.convert(y)]
def test_numexpr_printer():
if not numexpr:
skip("numexpr not installed.")
# if translation/printing is done incorrectly then evaluating
# a lambdified numexpr expression will throw an exception
from sympy.printing.lambdarepr import NumExprPrinter
from sympy import S
blacklist = ('where', 'complex', 'contains')
arg_tuple = (x, y, z) # some functions take more than one argument
for sym in NumExprPrinter._numexpr_functions.keys():
if sym in blacklist:
continue
ssym = S(sym)
if hasattr(ssym, '_nargs'):
nargs = ssym._nargs[0]
else:
nargs = 1
args = arg_tuple[:nargs]
f = lambdify(args, ssym(*args), modules='numexpr')
assert f(*(1, )*nargs) is not None
def test_issue9474():
mods = [None, 'math']
if numpy:
mods.append('numpy')
if mpmath:
mods.append('mpmath')
for mod in mods:
f = lambdify(x, sympy.S(1)/x, modules=mod)
assert f(2) == 0.5
f = lambdify(x, floor(sympy.S(1)/x), modules=mod)
assert f(2) == 0
if mpmath:
f = lambdify(x, sympy.S(1)/sympy.Abs(x), modules=['mpmath'])
assert isinstance(f(2), mpmath.mpf)
for absfunc, modules in product([Abs, abs], mods):
f = lambdify(x, absfunc(x), modules=modules)
assert f(-1) == 1
assert f(1) == 1
assert f(3+4j) == 5
def test_special_printers():
class IntervalPrinter(LambdaPrinter):
"""Use ``lambda`` printer but print numbers as ``mpi`` intervals. """
def _print_Integer(self, expr):
return "mpi('%s')" % super(IntervalPrinter, self)._print_Integer(expr)
def _print_Rational(self, expr):
return "mpi('%s')" % super(IntervalPrinter, self)._print_Rational(expr)
def intervalrepr(expr):
return IntervalPrinter().doprint(expr)
expr = sympy.sqrt(sympy.sqrt(2) + sympy.sqrt(3)) + sympy.S(1)/2
func0 = lambdify((), expr, modules="mpmath", printer=intervalrepr)
func1 = lambdify((), expr, modules="mpmath", printer=IntervalPrinter)
func2 = lambdify((), expr, modules="mpmath", printer=IntervalPrinter())
mpi = type(mpmath.mpi(1, 2))
assert isinstance(func0(), mpi)
assert isinstance(func1(), mpi)
assert isinstance(func2(), mpi)
def test_functions():
assert residue(1/sin(x), x, 0) == 1
assert residue(2/sin(x), x, 0) == 2
assert residue(1/sin(x)**2, x, 0) == 0
assert residue(1/sin(x)**5, x, 0) == S(3)/8
def test_expressions_failing():
assert residue(1/(x**4 + 1), x, exp(I*pi/4)) == -(S(1)/4 + I/4)/sqrt(2)
n = Symbol('n', integer=True, positive=True)
assert residue(exp(z)/(z - pi*I/4*a)**n, z, I*pi*a) == \
exp(I*pi*a/4)/factorial(n - 1)
assert residue(1/(x**2 + a**2)**2, x, a*I) == -I/4/a**3
def test_sqrtdenest():
d = {sqrt(5 + 2 * r6): r2 + r3,
sqrt(5. + 2 * r6): sqrt(5. + 2 * r6),
sqrt(5. + 4*sqrt(5 + 2 * r6)): sqrt(5.0 + 4*r2 + 4*r3),
sqrt(r2): sqrt(r2),
sqrt(5 + r7): sqrt(5 + r7),
sqrt(3 + sqrt(5 + 2*r7)):
3*r2*(5 + 2*r7)**(S(1)/4)/(2*sqrt(6 + 3*r7)) +
r2*sqrt(6 + 3*r7)/(2*(5 + 2*r7)**(S(1)/4)),
sqrt(3 + 2*r3): 3**(S(3)/4)*(r6/2 + 3*r2/2)/3}
for i in d:
assert sqrtdenest(i) == d[i]
