我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用sympy.S。
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_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_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_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_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_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]
def test_sqrtdenest_rec(): assert sqrtdenest(sqrt(-4*sqrt(14) - 2*r6 + 4*sqrt(21) + 33)) == \ -r2 + r3 + 2*r7 assert sqrtdenest(sqrt(-28*r7 - 14*r5 + 4*sqrt(35) + 82)) == \ -7 + r5 + 2*r7 assert sqrtdenest(sqrt(6*r2/11 + 2*sqrt(22)/11 + 6*sqrt(11)/11 + 2)) == \ sqrt(11)*(r2 + 3 + sqrt(11))/11 assert sqrtdenest(sqrt(468*r3 + 3024*r2 + 2912*r6 + 19735)) == \ 9*r3 + 26 + 56*r6 z = sqrt(-490*r3 - 98*sqrt(115) - 98*sqrt(345) - 2107) assert sqrtdenest(z) == sqrt(-1)*(7*r5 + 7*r15 + 7*sqrt(23)) z = sqrt(-4*sqrt(14) - 2*r6 + 4*sqrt(21) + 34) assert sqrtdenest(z) == z assert sqrtdenest(sqrt(-8*r2 - 2*r5 + 18)) == -r10 + 1 + r2 + r5 assert sqrtdenest(sqrt(8*r2 + 2*r5 - 18)) == \ sqrt(-1)*(-r10 + 1 + r2 + r5) assert sqrtdenest(sqrt(8*r2/3 + 14*r5/3 + S(154)/9)) == \ -r10/3 + r2 + r5 + 3 assert sqrtdenest(sqrt(sqrt(2*r6 + 5) + sqrt(2*r7 + 8))) == \ sqrt(1 + r2 + r3 + r7) assert sqrtdenest(sqrt(4*r15 + 8*r5 + 12*r3 + 24)) == 1 + r3 + r5 + r15 w = 1 + r2 + r3 + r5 + r7 assert sqrtdenest(sqrt((w**2).expand())) == w z = sqrt((w**2).expand() + 1) assert sqrtdenest(z) == z z = sqrt(2*r10 + 6*r2 + 4*r5 + 12 + 10*r15 + 30*r3) assert sqrtdenest(z) == z
def test_lattice_make_args(): assert join.make_args(0) == set([0]) assert join.make_args(1) == set([1]) assert join.make_args(join(2, 3, 4)) == set([S(2), S(3), S(4)])
def test_apply_finite_diff(): x, h = symbols('x h') f = Function('f') assert (apply_finite_diff(1, [x-h, x+h], [f(x-h), f(x+h)], x) - (f(x+h)-f(x-h))/(2*h)).simplify() == 0 assert (apply_finite_diff(1, [5, 6, 7], [f(5), f(6), f(7)], 5) - (-S(3)/2*f(5) + 2*f(6) - S(1)/2*f(7))).simplify() == 0
def test_gosper_normal(): assert gosper_normal(4*n + 5, 2*(4*n + 1)*(2*n + 3), n) == \ (Poly(S(1)/4, n), Poly(n + S(3)/2), Poly(n + S(1)/4))
def test_gosper_term(): assert gosper_term((4*k + 1)*factorial( k)/factorial(2*k + 1), k) == (-k - S(1)/2)/(k + S(1)/4)
def test_gosper_sum_parametric(): assert gosper_sum(binomial(S(1)/2, m - j + 1)*binomial(S(1)/2, m + j), (j, 1, n)) == \ n*(1 + m - n)*(-1 + 2*m + 2*n)*binomial(S(1)/2, 1 + m - n)* \ binomial(S(1)/2, m + n)/(m*(1 + 2*m))
def test_gosper_sum_AeqB_part1(): f1a = n**4 f1b = n**3*2**n f1c = 1/(n**2 + sqrt(5)*n - 1) f1d = n**4*4**n/binomial(2*n, n) f1e = factorial(3*n)/(factorial(n)*factorial(n + 1)*factorial(n + 2)*27**n) f1f = binomial(2*n, n)**2/((n + 1)*4**(2*n)) f1g = (4*n - 1)*binomial(2*n, n)**2/((2*n - 1)**2*4**(2*n)) f1h = n*factorial(n - S(1)/2)**2/factorial(n + 1)**2 g1a = m*(m + 1)*(2*m + 1)*(3*m**2 + 3*m - 1)/30 g1b = 26 + 2**(m + 1)*(m**3 - 3*m**2 + 9*m - 13) g1c = (m + 1)*(m*(m**2 - 7*m + 3)*sqrt(5) - ( 3*m**3 - 7*m**2 + 19*m - 6))/(2*m**3*sqrt(5) + m**4 + 5*m**2 - 1)/6 g1d = -S(2)/231 + 2*4**m*(m + 1)*(63*m**4 + 112*m**3 + 18*m**2 - 22*m + 3)/(693*binomial(2*m, m)) g1e = -S(9)/2 + (81*m**2 + 261*m + 200)*factorial( 3*m + 2)/(40*27**m*factorial(m)*factorial(m + 1)*factorial(m + 2)) g1f = (2*m + 1)**2*binomial(2*m, m)**2/(4**(2*m)*(m + 1)) g1g = -binomial(2*m, m)**2/4**(2*m) g1h = 4*pi -(2*m + 1)**2*(3*m + 4)*factorial(m - S(1)/2)**2/factorial(m + 1)**2 g = gosper_sum(f1a, (n, 0, m)) assert g is not None and simplify(g - g1a) == 0 g = gosper_sum(f1b, (n, 0, m)) assert g is not None and simplify(g - g1b) == 0 g = gosper_sum(f1c, (n, 0, m)) assert g is not None and simplify(g - g1c) == 0 g = gosper_sum(f1d, (n, 0, m)) assert g is not None and simplify(g - g1d) == 0 g = gosper_sum(f1e, (n, 0, m)) assert g is not None and simplify(g - g1e) == 0 g = gosper_sum(f1f, (n, 0, m)) assert g is not None and simplify(g - g1f) == 0 g = gosper_sum(f1g, (n, 0, m)) assert g is not None and simplify(g - g1g) == 0 g = gosper_sum(f1h, (n, 0, m)) # need to call rewrite(gamma) here because we have terms involving # factorial(1/2) assert g is not None and simplify(g - g1h).rewrite(gamma) == 0
def symbol_array(base, n=None): """ Generates a string arrary with str_array and replaces each string in array with Symbol of same name. """ symbol_str_lst = str_array(base, n) result = [] for symbol_str in symbol_str_lst: result.append(S(symbol_str)) return tuple(result)
def is_pell_transformation_ok(eq): """ Test whether X*Y, X, or Y terms are present in the equation after transforming the equation using the transformation returned by transformation_to_pell(). If they are not present we are good. Moreover, coefficient of X**2 should be a divisor of coefficient of Y**2 and the constant term. """ A, B = transformation_to_DN(eq) u = (A*Matrix([X, Y]) + B)[0] v = (A*Matrix([X, Y]) + B)[1] simplified = _mexpand(Subs(eq, (x, y), (u, v)).doit()) coeff = dict([reversed(t.as_independent(*[X, Y])) for t in simplified.args]) for term in [X*Y, X, Y]: if term in coeff.keys(): return False for term in [X**2, Y**2, Integer(1)]: if term not in coeff.keys(): coeff[term] = Integer(0) if coeff[X**2] != 0: return isinstance(S(coeff[Y**2])/coeff[X**2], Integer) and isinstance(S(coeff[Integer(1)])/coeff[X**2], Integer) return True
def E_nl(n, Z=1): """ Returns the energy of the state (n, l) in Hartree atomic units. The energy doesn't depend on "l". Examples ======== >>> from sympy import var >>> from sympy.physics.hydrogen import E_nl >>> var("n Z") (n, Z) >>> E_nl(n, Z) -Z**2/(2*n**2) >>> E_nl(1) -1/2 >>> E_nl(2) -1/8 >>> E_nl(3) -1/18 >>> E_nl(3, 47) -2209/18 """ n, Z = S(n), S(Z) if n.is_integer and (n < 1): raise ValueError("'n' must be positive integer") return -Z**2/(2*n**2)
def test_represent_rotation(): assert represent(Rotation(0, pi/2, 0)) == \ Matrix( [[WignerD( S( 1)/2, S( 1)/2, S( 1)/2, 0, pi/2, 0), WignerD( S(1)/2, S(1)/2, -S(1)/2, 0, pi/2, 0)], [WignerD(S(1)/2, -S(1)/2, S(1)/2, 0, pi/2, 0), WignerD(S(1)/2, -S(1)/2, -S(1)/2, 0, pi/2, 0)]]) assert represent(Rotation(0, pi/2, 0), doit=True) == \ Matrix([[sqrt(2)/2, -sqrt(2)/2], [sqrt(2)/2, sqrt(2)/2]])
def test_innerproduct(): assert InnerProduct(JzBra(1, 1), JzKet(1, 1)).doit() == 1 assert InnerProduct( JzBra(S(1)/2, S(1)/2), JzKet(S(1)/2, -S(1)/2)).doit() == 0 assert InnerProduct(JzBra(j, m), JzKet(j, m)).doit() == 1 assert InnerProduct(JzBra(1, 0), JyKet(1, 1)).doit() == I/sqrt(2) assert InnerProduct( JxBra(S(1)/2, S(1)/2), JzKet(S(1)/2, S(1)/2)).doit() == -sqrt(2)/2 assert InnerProduct(JyBra(1, 1), JzKet(1, 1)).doit() == S(1)/2 assert InnerProduct(JxBra(1, -1), JyKet(1, 1)).doit() == 0
def test_jzket(): j, m = symbols('j m') # j not integer or half integer raises(ValueError, lambda: JzKet(S(2)/3, -S(1)/3)) raises(ValueError, lambda: JzKet(S(2)/3, m)) # j < 0 raises(ValueError, lambda: JzKet(-1, 1)) raises(ValueError, lambda: JzKet(-1, m)) # m not integer or half integer raises(ValueError, lambda: JzKet(j, -S(1)/3)) # abs(m) > j raises(ValueError, lambda: JzKet(1, 2)) raises(ValueError, lambda: JzKet(1, -2)) # j-m not integer raises(ValueError, lambda: JzKet(1, S(1)/2))
def test_doit(): assert Wigner3j(1/2, -1/2, 1/2, 1/2, 0, 0).doit() == -sqrt(2)/2 assert Wigner6j(1, 2, 3, 2, 1, 2).doit() == sqrt(21)/105 assert Wigner9j( 2, 1, 1, S(3)/2, S(1)/2, 1, S(1)/2, S(1)/2, 0).doit() == sqrt(2)/12 assert CG(1/2, 1/2, 1/2, -1/2, 1, 0).doit() == sqrt(2)/2
def test_entropy(): up = JzKet(S(1)/2, S(1)/2) down = JzKet(S(1)/2, -S(1)/2) d = Density((up, 0.5), (down, 0.5)) # test for density object ent = entropy(d) assert entropy(d) == 0.5*log(2) assert d.entropy() == 0.5*log(2) np = import_module('numpy', min_module_version='1.4.0') if np: #do this test only if 'numpy' is available on test machine np_mat = represent(d, format='numpy') ent = entropy(np_mat) assert isinstance(np_mat, np.matrixlib.defmatrix.matrix) assert ent.real == 0.69314718055994529 assert ent.imag == 0 scipy = import_module('scipy', __import__kwargs={'fromlist': ['sparse']}) if scipy and np: #do this test only if numpy and scipy are available mat = represent(d, format="scipy.sparse") assert isinstance(mat, scipy_sparse_matrix) assert ent.real == 0.69314718055994529 assert ent.imag == 0
def test_clebsch_gordan_docs(): assert clebsch_gordan(S(3)/2, S(1)/2, 2, S(3)/2, S(1)/2, 2) == 1 assert clebsch_gordan(S(3)/2, S(1)/2, 1, S(3)/2, -S(1)/2, 1) == sqrt(3)/2 assert clebsch_gordan(S(3)/2, S(1)/2, 1, -S(1)/2, S(1)/2, 0) == -sqrt(2)/2
def test_wigner(): def tn(a, b): return (a - b).n(64) < S('1e-64') assert tn(wigner_9j(1, 1, 1, 1, 1, 1, 1, 1, 0, prec=64), S(1)/18) assert wigner_9j(3, 3, 2, 3, 3, 2, 3, 3, 2) == 3221*sqrt( 70)/(246960*sqrt(105)) - 365/(3528*sqrt(70)*sqrt(105)) assert wigner_6j(5, 5, 5, 5, 5, 5) == Rational(1, 52) assert tn(wigner_6j(8, 8, 8, 8, 8, 8, prec=64), -S(12219)/965770)
def test_gaunt(): def tn(a, b): return (a - b).n(64) < S('1e-64') assert gaunt(1, 0, 1, 1, 0, -1) == -1/(2*sqrt(pi)) assert tn(gaunt( 10, 10, 12, 9, 3, -12, prec=64), (-S(98)/62031) * sqrt(6279)/sqrt(pi))
def test_racah(): assert racah(3,3,3,3,3,3) == Rational(-1,14) assert racah(2,2,2,2,2,2) == Rational(-3,70) assert racah(7,8,7,1,7,7, prec=4).is_Float assert racah(5.5,7.5,9.5,6.5,8,9) == -719*sqrt(598)/1158924 assert abs(racah(5.5,7.5,9.5,6.5,8,9, prec=4) - (-0.01517)) < S('1e-4')
def test_hydrogen_energies(): assert E_nl(n, Z) == -Z**2/(2*n**2) assert E_nl(n) == -1/(2*n**2) assert E_nl(1, 47) == -S(47)**2/(2*1**2) assert E_nl(2, 47) == -S(47)**2/(2*2**2) assert E_nl(1) == -S(1)/(2*1**2) assert E_nl(2) == -S(1)/(2*2**2) assert E_nl(3) == -S(1)/(2*3**2) assert E_nl(4) == -S(1)/(2*4**2) assert E_nl(100) == -S(1)/(2*100**2) raises(ValueError, lambda: E_nl(0))
