diff --git a/.mailmap b/.mailmap index d11f98085f3d..2684bc38b47f 100644 --- a/.mailmap +++ b/.mailmap @@ -547,3 +547,4 @@ Prince Gupta LAPTOP-AS1M2R8B\codem MaqOwais Vaibhav Bhat VBhat97 Bhaskar Joshi BhaskarJoshi-01 +Miguel Torres Costa diff --git a/.travis.yml b/.travis.yml index 87465944b534..6ea942839c2a 100644 --- a/.travis.yml +++ b/.travis.yml @@ -33,7 +33,7 @@ jobs: # all of the dependencies are supported on 3.8. env: - TEST_ASCII="true" - - TEST_OPT_DEPENDENCY="matchpy numpy scipy gmpy2 matplotlib aesara llvmlite autowrap cython wurlitzer python-symengine=0.5.1 tensorflow numexpr ipython antlr-python-runtime>=4.7,<4.8 antlr>=4.7,<4.8 cloudpickle pyglet pycosat lfortran python-clang lxml" + - TEST_OPT_DEPENDENCY="matchpy numpy scipy gmpy2 matplotlib aesara llvmlite autowrap cython wurlitzer python-symengine=0.5.1 tensorflow numexpr ipython antlr-python-runtime>=4.7,<4.8 antlr>=4.7,<4.8 cloudpickle pyglet pycosat python-sat lfortran python-clang lxml" - TEST_SAGE="true" - SYMPY_STRICT_COMPILER_CHECKS=1 addons: diff --git a/README.md b/README.md index 1ecc50422774..859a90a18dd2 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ [![Zenodo Badge](https://zenodo.org/badge/18918/sympy/sympy.svg)](https://zenodo.org/badge/latestdoi/18918/sympy/sympy) [![codecov Badge](https://codecov.io/gh/sympy/sympy/branch/master/graph/badge.svg)](https://codecov.io/gh/sympy/sympy) -[![SymPy Banner](banner.svg)](https://sympy.org/) +[![SymPy Banner](https://github.com/sympy/sympy/raw/master/banner.svg)](https://sympy.org/) See the AUTHORS file for the list of authors. diff --git a/bin/test_external_imports.py b/bin/test_external_imports.py index 9d77de4ed119..62e8df6bcf22 100755 --- a/bin/test_external_imports.py +++ b/bin/test_external_imports.py @@ -25,7 +25,7 @@ # These libraries are optional, but are always imported at SymPy import time # when they are installed. External libraries should only be added to this # list if they are required for core SymPy functionality. -hard_optional_dependencies = ['gmpy', 'gmpy2', 'pycosat'] +hard_optional_dependencies = ['gmpy', 'gmpy2', 'pycosat', 'python-sat'] import sys import os diff --git a/bin/test_submodule_imports.py b/bin/test_submodule_imports.py index 86df5bcea95c..1fb634ab5f64 100755 --- a/bin/test_submodule_imports.py +++ b/bin/test_submodule_imports.py @@ -30,7 +30,6 @@ 'deprecated', 'discrete', 'external', - 'equation', 'functions', 'geometry', 'integrals', diff --git a/doc/src/modules/physics/control/lti.rst b/doc/src/modules/physics/control/lti.rst index cc9ae87bcd64..4d683a4ffa8c 100644 --- a/doc/src/modules/physics/control/lti.rst +++ b/doc/src/modules/physics/control/lti.rst @@ -5,5 +5,16 @@ Control API lti === -.. automodule:: sympy.physics.control.lti +.. module:: sympy.physics.control.lti + +.. autoclass:: TransferFunction + :members: + +.. autoclass:: Series + :members: + +.. autoclass:: Parallel + :members: + +.. autoclass:: Feedback :members: diff --git a/doc/src/modules/solvers/diophantine.rst b/doc/src/modules/solvers/diophantine.rst index 5afe1f8d0f5a..06edbfc73f6e 100644 --- a/doc/src/modules/solvers/diophantine.rst +++ b/doc/src/modules/solvers/diophantine.rst @@ -324,7 +324,7 @@ References .. [3] Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0,[online], Available: http://www.alpertron.com.ar/METHODS.HTM .. [4] Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online], - Available: http://www.jpr2718.org/ax2p.pdf + Available: https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf User Functions -------------- diff --git a/doc/src/modules/solvers/ode.rst b/doc/src/modules/solvers/ode.rst index 4322f08069b0..1bc87c03bad8 100644 --- a/doc/src/modules/solvers/ode.rst +++ b/doc/src/modules/solvers/ode.rst @@ -111,7 +111,8 @@ Bernoulli Liouville ^^^^^^^^^ -.. autofunction:: sympy.solvers.ode.ode::ode_Liouville +.. autoclass:: sympy.solvers.ode.single::Liouville + :members: Riccati_special_minus2 ^^^^^^^^^^^^^^^^^^^^^^ @@ -141,7 +142,8 @@ nth_order_reducible separable ^^^^^^^^^ -.. autofunction:: sympy.solvers.ode.ode::ode_separable +.. autoclass:: sympy.solvers.ode.single::Separable + :members: almost_linear ^^^^^^^^^^^^^ @@ -154,7 +156,7 @@ linear_coefficients separable_reduced ^^^^^^^^^^^^^^^^^ -.. autofunction:: sympy.solvers.ode.ode::ode_separable_reduced +.. autofunction:: sympy.solvers.ode.single::SeparableReduced lie_group ^^^^^^^^^ diff --git a/setup.cfg b/setup.cfg index a440ccbcf088..f37df4129444 100644 --- a/setup.cfg +++ b/setup.cfg @@ -58,6 +58,8 @@ ignore_missing_imports = True ignore_missing_imports = True [mypy-pymc3.*] ignore_missing_imports = True +[mypy-python-sat.*] +ignore_missing_imports = True [mypy-pytest.*] ignore_missing_imports = True [mypy-_pytest.*] diff --git a/setup.py b/setup.py index 54521df1c3c4..615052e87837 100755 --- a/setup.py +++ b/setup.py @@ -93,7 +93,6 @@ 'sympy.crypto', 'sympy.diffgeom', 'sympy.discrete', - 'sympy.equation', 'sympy.external', 'sympy.functions', 'sympy.functions.combinatorial', @@ -359,7 +358,6 @@ def run(self): 'sympy.crypto.tests', 'sympy.diffgeom.tests', 'sympy.discrete.tests', - 'sympy.equation.tests', 'sympy.external.tests', 'sympy.functions.combinatorial.tests', 'sympy.functions.elementary.tests', diff --git a/sympy/__init__.py b/sympy/__init__.py index ce2d8caf4aca..dadf42ea9616 100644 --- a/sympy/__init__.py +++ b/sympy/__init__.py @@ -152,8 +152,6 @@ def __sympy_debug(): inverse_mobius_transform, convolution, covering_product, intersecting_product) -from .equation import (Equation, Eqn,) - from .simplify import (simplify, hypersimp, hypersimilar, logcombine, separatevars, posify, besselsimp, kroneckersimp, signsimp, bottom_up, nsimplify, FU, fu, sqrtdenest, cse, use, epath, EPath, hyperexpand, @@ -378,9 +376,6 @@ def __sympy_debug(): 'inverse_mobius_transform', 'convolution', 'covering_product', 'intersecting_product', - # sympy.equation - 'Equation', 'Eqn', - # sympy.simplify 'simplify', 'hypersimp', 'hypersimilar', 'logcombine', 'separatevars', 'posify', 'besselsimp', 'kroneckersimp', 'signsimp', 'bottom_up', diff --git a/sympy/calculus/tests/test_util.py b/sympy/calculus/tests/test_util.py index 41514e7210c0..97799ce521d3 100644 --- a/sympy/calculus/tests/test_util.py +++ b/sympy/calculus/tests/test_util.py @@ -5,13 +5,14 @@ periodicity, lcim, AccumBounds, is_convex, stationary_points, minimum, maximum) from sympy.core import Add, Mul, Pow +from sympy.core.expr import unchanged from sympy.sets.sets import (Interval, FiniteSet, EmptySet, Complement, Union) -from sympy.testing.pytest import raises, _both_exp_pow +from sympy.testing.pytest import raises, _both_exp_pow, XFAIL from sympy.abc import x a = Symbol('a', real=True) - +B = AccumBounds def test_function_range(): x, y, a, b = symbols('x y a b') @@ -310,253 +311,325 @@ def test_issue_19869(): def test_AccumBounds(): - assert AccumBounds(1, 2).args == (1, 2) - assert AccumBounds(1, 2).delta is S.One - assert AccumBounds(1, 2).mid == Rational(3, 2) - assert AccumBounds(1, 3).is_real == True - - assert AccumBounds(1, 1) is S.One - - assert AccumBounds(1, 2) + 1 == AccumBounds(2, 3) - assert 1 + AccumBounds(1, 2) == AccumBounds(2, 3) - assert AccumBounds(1, 2) + AccumBounds(2, 3) == AccumBounds(3, 5) - - assert -AccumBounds(1, 2) == AccumBounds(-2, -1) - - assert AccumBounds(1, 2) - 1 == AccumBounds(0, 1) - assert 1 - AccumBounds(1, 2) == AccumBounds(-1, 0) - assert AccumBounds(2, 3) - AccumBounds(1, 2) == AccumBounds(0, 2) - - assert x + AccumBounds(1, 2) == Add(AccumBounds(1, 2), x) - assert a + AccumBounds(1, 2) == AccumBounds(1 + a, 2 + a) - assert AccumBounds(1, 2) - x == Add(AccumBounds(1, 2), -x) - - assert AccumBounds(-oo, 1) + oo == AccumBounds(-oo, oo) - assert AccumBounds(1, oo) + oo is oo - assert AccumBounds(1, oo) - oo == AccumBounds(-oo, oo) - assert (-oo - AccumBounds(-1, oo)) is -oo - assert AccumBounds(-oo, 1) - oo is -oo - - assert AccumBounds(1, oo) - oo == AccumBounds(-oo, oo) - assert AccumBounds(-oo, 1) - (-oo) == AccumBounds(-oo, oo) - assert (oo - AccumBounds(1, oo)) == AccumBounds(-oo, oo) - assert (-oo - AccumBounds(1, oo)) is -oo - - assert AccumBounds(1, 2)/2 == AccumBounds(S.Half, 1) - assert 2/AccumBounds(2, 3) == AccumBounds(Rational(2, 3), 1) - assert 1/AccumBounds(-1, 1) == AccumBounds(-oo, oo) - - assert abs(AccumBounds(1, 2)) == AccumBounds(1, 2) - assert abs(AccumBounds(-2, -1)) == AccumBounds(1, 2) - assert abs(AccumBounds(-2, 1)) == AccumBounds(0, 2) - assert abs(AccumBounds(-1, 2)) == AccumBounds(0, 2) + assert B(1, 2).args == (1, 2) + assert B(1, 2).delta is S.One + assert B(1, 2).mid == Rational(3, 2) + assert B(1, 3).is_real == True + + assert B(1, 1) is S.One + + assert B(1, 2) + 1 == B(2, 3) + assert 1 + B(1, 2) == B(2, 3) + assert B(1, 2) + B(2, 3) == B(3, 5) + + assert -B(1, 2) == B(-2, -1) + + assert B(1, 2) - 1 == B(0, 1) + assert 1 - B(1, 2) == B(-1, 0) + assert B(2, 3) - B(1, 2) == B(0, 2) + + assert x + B(1, 2) == Add(B(1, 2), x) + assert a + B(1, 2) == B(1 + a, 2 + a) + assert B(1, 2) - x == Add(B(1, 2), -x) + + assert B(-oo, 1) + oo == B(-oo, oo) + assert B(1, oo) + oo is oo + assert B(1, oo) - oo == B(-oo, oo) + assert (-oo - B(-1, oo)) is -oo + assert B(-oo, 1) - oo is -oo + + assert B(1, oo) - oo == B(-oo, oo) + assert B(-oo, 1) - (-oo) == B(-oo, oo) + assert (oo - B(1, oo)) == B(-oo, oo) + assert (-oo - B(1, oo)) is -oo + + assert B(1, 2)/2 == B(S.Half, 1) + assert 2/B(2, 3) == B(Rational(2, 3), 1) + assert 1/B(-1, 1) == B(-oo, oo) + + assert abs(B(1, 2)) == B(1, 2) + assert abs(B(-2, -1)) == B(1, 2) + assert abs(B(-2, 1)) == B(0, 2) + assert abs(B(-1, 2)) == B(0, 2) c = Symbol('c') - raises(ValueError, lambda: AccumBounds(0, c)) - raises(ValueError, lambda: AccumBounds(1, -1)) + raises(ValueError, lambda: B(0, c)) + raises(ValueError, lambda: B(1, -1)) + r = Symbol('r', real=True) + raises(ValueError, lambda: B(r, r - 1)) def test_AccumBounds_mul(): - assert AccumBounds(1, 2)*2 == AccumBounds(2, 4) - assert 2*AccumBounds(1, 2) == AccumBounds(2, 4) - assert AccumBounds(1, 2)*AccumBounds(2, 3) == AccumBounds(2, 6) - - assert AccumBounds(1, 2)*0 == 0 - assert AccumBounds(1, oo)*0 == AccumBounds(0, oo) - assert AccumBounds(-oo, 1)*0 == AccumBounds(-oo, 0) - assert AccumBounds(-oo, oo)*0 == AccumBounds(-oo, oo) - - assert AccumBounds(1, 2)*x == Mul(AccumBounds(1, 2), x, evaluate=False) - - assert AccumBounds(0, 2)*oo == AccumBounds(0, oo) - assert AccumBounds(-2, 0)*oo == AccumBounds(-oo, 0) - assert AccumBounds(0, 2)*(-oo) == AccumBounds(-oo, 0) - assert AccumBounds(-2, 0)*(-oo) == AccumBounds(0, oo) - assert AccumBounds(-1, 1)*oo == AccumBounds(-oo, oo) - assert AccumBounds(-1, 1)*(-oo) == AccumBounds(-oo, oo) - assert AccumBounds(-oo, oo)*oo == AccumBounds(-oo, oo) + assert B(1, 2)*2 == B(2, 4) + assert 2*B(1, 2) == B(2, 4) + assert B(1, 2)*B(2, 3) == B(2, 6) + assert B(0, 2)*B(2, oo) == B(0, oo) + l, r = B(-oo, oo), B(-a, a) + assert l*r == B(-oo, oo) + assert r*l == B(-oo, oo) + l, r = B(1, oo), B(-3, -2) + assert l*r == B(-oo, -2) + assert r*l == B(-oo, -2) + assert B(1, 2)*0 == 0 + assert B(1, oo)*0 == B(0, oo) + assert B(-oo, 1)*0 == B(-oo, 0) + assert B(-oo, oo)*0 == B(-oo, oo) + + assert B(1, 2)*x == Mul(B(1, 2), x, evaluate=False) + + assert B(0, 2)*oo == B(0, oo) + assert B(-2, 0)*oo == B(-oo, 0) + assert B(0, 2)*(-oo) == B(-oo, 0) + assert B(-2, 0)*(-oo) == B(0, oo) + assert B(-1, 1)*oo == B(-oo, oo) + assert B(-1, 1)*(-oo) == B(-oo, oo) + assert B(-oo, oo)*oo == B(-oo, oo) def test_AccumBounds_div(): - assert AccumBounds(-1, 3)/AccumBounds(3, 4) == AccumBounds(Rational(-1, 3), 1) - assert AccumBounds(-2, 4)/AccumBounds(-3, 4) == AccumBounds(-oo, oo) - assert AccumBounds(-3, -2)/AccumBounds(-4, 0) == AccumBounds(S.Half, oo) + assert B(-1, 3)/B(3, 4) == B(Rational(-1, 3), 1) + assert B(-2, 4)/B(-3, 4) == B(-oo, oo) + assert B(-3, -2)/B(-4, 0) == B(S.Half, oo) # these two tests can have a better answer - # after Union of AccumBounds is improved - assert AccumBounds(-3, -2)/AccumBounds(-2, 1) == AccumBounds(-oo, oo) - assert AccumBounds(2, 3)/AccumBounds(-2, 2) == AccumBounds(-oo, oo) - - assert AccumBounds(-3, -2)/AccumBounds(0, 4) == AccumBounds(-oo, Rational(-1, 2)) - assert AccumBounds(2, 4)/AccumBounds(-3, 0) == AccumBounds(-oo, Rational(-2, 3)) - assert AccumBounds(2, 4)/AccumBounds(0, 3) == AccumBounds(Rational(2, 3), oo) - - assert AccumBounds(0, 1)/AccumBounds(0, 1) == AccumBounds(0, oo) - assert AccumBounds(-1, 0)/AccumBounds(0, 1) == AccumBounds(-oo, 0) - assert AccumBounds(-1, 2)/AccumBounds(-2, 2) == AccumBounds(-oo, oo) - - assert 1/AccumBounds(-1, 2) == AccumBounds(-oo, oo) - assert 1/AccumBounds(0, 2) == AccumBounds(S.Half, oo) - assert (-1)/AccumBounds(0, 2) == AccumBounds(-oo, Rational(-1, 2)) - assert 1/AccumBounds(-oo, 0) == AccumBounds(-oo, 0) - assert 1/AccumBounds(-1, 0) == AccumBounds(-oo, -1) - assert (-2)/AccumBounds(-oo, 0) == AccumBounds(0, oo) - assert 1/AccumBounds(-oo, -1) == AccumBounds(-1, 0) - - assert AccumBounds(1, 2)/a == Mul(AccumBounds(1, 2), 1/a, evaluate=False) - - assert AccumBounds(1, 2)/0 == AccumBounds(1, 2)*zoo - assert AccumBounds(1, oo)/oo == AccumBounds(0, oo) - assert AccumBounds(1, oo)/(-oo) == AccumBounds(-oo, 0) - assert AccumBounds(-oo, -1)/oo == AccumBounds(-oo, 0) - assert AccumBounds(-oo, -1)/(-oo) == AccumBounds(0, oo) - assert AccumBounds(-oo, oo)/oo == AccumBounds(-oo, oo) - assert AccumBounds(-oo, oo)/(-oo) == AccumBounds(-oo, oo) - assert AccumBounds(-1, oo)/oo == AccumBounds(0, oo) - assert AccumBounds(-1, oo)/(-oo) == AccumBounds(-oo, 0) - assert AccumBounds(-oo, 1)/oo == AccumBounds(-oo, 0) - assert AccumBounds(-oo, 1)/(-oo) == AccumBounds(0, oo) + # after Union of B is improved + assert B(-3, -2)/B(-2, 1) == B(-oo, oo) + assert B(2, 3)/B(-2, 2) == B(-oo, oo) + + assert B(-3, -2)/B(0, 4) == B(-oo, Rational(-1, 2)) + assert B(2, 4)/B(-3, 0) == B(-oo, Rational(-2, 3)) + assert B(2, 4)/B(0, 3) == B(Rational(2, 3), oo) + + assert B(0, 1)/B(0, 1) == B(0, oo) + assert B(-1, 0)/B(0, 1) == B(-oo, 0) + assert B(-1, 2)/B(-2, 2) == B(-oo, oo) + + assert 1/B(-1, 2) == B(-oo, oo) + assert 1/B(0, 2) == B(S.Half, oo) + assert (-1)/B(0, 2) == B(-oo, Rational(-1, 2)) + assert 1/B(-oo, 0) == B(-oo, 0) + assert 1/B(-1, 0) == B(-oo, -1) + assert (-2)/B(-oo, 0) == B(0, oo) + assert 1/B(-oo, -1) == B(-1, 0) + + assert B(1, 2)/a == Mul(B(1, 2), 1/a, evaluate=False) + + assert B(1, 2)/0 == B(1, 2)*zoo + assert B(1, oo)/oo == B(0, oo) + assert B(1, oo)/(-oo) == B(-oo, 0) + assert B(-oo, -1)/oo == B(-oo, 0) + assert B(-oo, -1)/(-oo) == B(0, oo) + assert B(-oo, oo)/oo == B(-oo, oo) + assert B(-oo, oo)/(-oo) == B(-oo, oo) + assert B(-1, oo)/oo == B(0, oo) + assert B(-1, oo)/(-oo) == B(-oo, 0) + assert B(-oo, 1)/oo == B(-oo, 0) + assert B(-oo, 1)/(-oo) == B(0, oo) + def test_issue_18795(): r = Symbol('r', real=True) - a = AccumBounds(-1,1) - c = AccumBounds(7, oo) - b = AccumBounds(-oo, oo) - assert c - tan(r) == AccumBounds(7-tan(r), oo) - assert b + tan(r) == AccumBounds(-oo, oo) - assert (a + r)/a == AccumBounds(-oo, oo)*AccumBounds(r - 1, r + 1) - assert (b + a)/a == AccumBounds(-oo, oo) + a = B(-1,1) + c = B(7, oo) + b = B(-oo, oo) + assert c - tan(r) == B(7-tan(r), oo) + assert b + tan(r) == B(-oo, oo) + assert (a + r)/a == B(-oo, oo)*B(r - 1, r + 1) + assert (b + a)/a == B(-oo, oo) + def test_AccumBounds_func(): - assert (x**2 + 2*x + 1).subs(x, AccumBounds(-1, 1)) == AccumBounds(-1, 4) - assert exp(AccumBounds(0, 1)) == AccumBounds(1, E) - assert exp(AccumBounds(-oo, oo)) == AccumBounds(0, oo) - assert log(AccumBounds(3, 6)) == AccumBounds(log(3), log(6)) + assert (x**2 + 2*x + 1).subs(x, B(-1, 1)) == B(-1, 4) + assert exp(B(0, 1)) == B(1, E) + assert exp(B(-oo, oo)) == B(0, oo) + assert log(B(3, 6)) == B(log(3), log(6)) + + +@XFAIL +def test_AccumBounds_powf(): + nn = Symbol('nn', nonnegative=True) + assert B(1 + nn, 2 + nn)**B(1, 2) == B(1 + nn, (2 + nn)**2) + i = Symbol('i', integer=True, negative=True) + assert B(1, 2)**i == B(2**i, 1) def test_AccumBounds_pow(): - assert AccumBounds(0, 2)**2 == AccumBounds(0, 4) - assert AccumBounds(-1, 1)**2 == AccumBounds(0, 1) - assert AccumBounds(1, 2)**2 == AccumBounds(1, 4) - assert AccumBounds(-1, 2)**3 == AccumBounds(-1, 8) - assert AccumBounds(-1, 1)**0 == 1 - - assert AccumBounds(1, 2)**Rational(5, 2) == AccumBounds(1, 4*sqrt(2)) - assert AccumBounds(-1, 2)**Rational(1, 3) == AccumBounds(-1, 2**Rational(1, 3)) - assert AccumBounds(0, 2)**S.Half == AccumBounds(0, sqrt(2)) - - assert AccumBounds(-4, 2)**Rational(2, 3) == AccumBounds(0, 2*2**Rational(1, 3)) - - assert AccumBounds(-1, 5)**S.Half == AccumBounds(0, sqrt(5)) - assert AccumBounds(-oo, 2)**S.Half == AccumBounds(0, sqrt(2)) - assert AccumBounds(-2, 3)**Rational(-1, 4) == AccumBounds(0, oo) - - assert AccumBounds(1, 5)**(-2) == AccumBounds(Rational(1, 25), 1) - assert AccumBounds(-1, 3)**(-2) == AccumBounds(0, oo) - assert AccumBounds(0, 2)**(-2) == AccumBounds(Rational(1, 4), oo) - assert AccumBounds(-1, 2)**(-3) == AccumBounds(-oo, oo) - assert AccumBounds(-3, -2)**(-3) == AccumBounds(Rational(-1, 8), Rational(-1, 27)) - assert AccumBounds(-3, -2)**(-2) == AccumBounds(Rational(1, 9), Rational(1, 4)) - assert AccumBounds(0, oo)**S.Half == AccumBounds(0, oo) - assert AccumBounds(-oo, -1)**Rational(1, 3) == AccumBounds(-oo, -1) - assert AccumBounds(-2, 3)**(Rational(-1, 3)) == AccumBounds(-oo, oo) - assert AccumBounds(-oo, 0)**(-2) == AccumBounds(0, oo) - assert AccumBounds(-2, 0)**(-2) == AccumBounds(Rational(1, 4), oo) - - assert AccumBounds(Rational(1, 3), S.Half)**oo is S.Zero - assert AccumBounds(0, S.Half)**oo is S.Zero - assert AccumBounds(S.Half, 1)**oo == AccumBounds(0, oo) - assert AccumBounds(0, 1)**oo == AccumBounds(0, oo) - assert AccumBounds(2, 3)**oo is oo - assert AccumBounds(1, 2)**oo == AccumBounds(0, oo) - assert AccumBounds(S.Half, 3)**oo == AccumBounds(0, oo) - assert AccumBounds(Rational(-1, 3), Rational(-1, 4))**oo is S.Zero - assert AccumBounds(-1, Rational(-1, 2))**oo == AccumBounds(-oo, oo) - assert AccumBounds(-3, -2)**oo == FiniteSet(-oo, oo) - assert AccumBounds(-2, -1)**oo == AccumBounds(-oo, oo) - assert AccumBounds(-2, Rational(-1, 2))**oo == AccumBounds(-oo, oo) - assert AccumBounds(Rational(-1, 2), S.Half)**oo is S.Zero - assert AccumBounds(Rational(-1, 2), 1)**oo == AccumBounds(0, oo) - assert AccumBounds(Rational(-2, 3), 2)**oo == AccumBounds(0, oo) - assert AccumBounds(-1, 1)**oo == AccumBounds(-oo, oo) - assert AccumBounds(-1, S.Half)**oo == AccumBounds(-oo, oo) - assert AccumBounds(-1, 2)**oo == AccumBounds(-oo, oo) - assert AccumBounds(-2, S.Half)**oo == AccumBounds(-oo, oo) - - assert AccumBounds(1, 2)**x == Pow(AccumBounds(1, 2), x) - - assert AccumBounds(2, 3)**(-oo) is S.Zero - assert AccumBounds(0, 2)**(-oo) == AccumBounds(0, oo) - assert AccumBounds(-1, 2)**(-oo) == AccumBounds(-oo, oo) - - assert (tan(x)**sin(2*x)).subs(x, AccumBounds(0, pi/2)) == \ - Pow(AccumBounds(-oo, oo), AccumBounds(0, 1)) + assert B(0, 2)**2 == B(0, 4) + assert B(-1, 1)**2 == B(0, 1) + assert B(1, 2)**2 == B(1, 4) + assert B(-1, 2)**3 == B(-1, 8) + assert B(-1, 1)**0 == 1 + + assert B(1, 2)**Rational(5, 2) == B(1, 4*sqrt(2)) + assert B(0, 2)**S.Half == B(0, sqrt(2)) + + neg = Symbol('neg', negative=True) + unchanged(Pow, B(neg, 1), S.Half) + nn = Symbol('nn', nonnegative=True) + assert B(nn, nn + 1)**S.Half == B(sqrt(nn), sqrt(nn + 1)) + assert B(nn, nn + 1)**nn == B(nn**nn, (nn + 1)**nn) + unchanged(Pow, B(nn, nn + 1), x) + i = Symbol('i', integer=True) + unchanged(Pow, B(1, 2), i) + i = Symbol('i', integer=True, nonnegative=True) + assert B(1, 2)**i == B(1, 2**i) + assert B(0, 1)**i == B(0**i, 1) + + assert B(1, 5)**(-2) == B(Rational(1, 25), 1) + assert B(-1, 3)**(-2) == B(0, oo) + assert B(0, 2)**(-3) == B(Rational(1, 8), oo) + assert B(-2, 0)**(-3) == B(-oo, -Rational(1, 8)) + assert B(0, 2)**(-2) == B(Rational(1, 4), oo) + assert B(-1, 2)**(-3) == B(-oo, oo) + assert B(-3, -2)**(-3) == B(Rational(-1, 8), Rational(-1, 27)) + assert B(-3, -2)**(-2) == B(Rational(1, 9), Rational(1, 4)) + assert B(0, oo)**S.Half == B(0, oo) + assert B(-oo, 0)**(-2) == B(0, oo) + assert B(-2, 0)**(-2) == B(Rational(1, 4), oo) + + assert B(Rational(1, 3), S.Half)**oo is S.Zero + assert B(0, S.Half)**oo is S.Zero + assert B(S.Half, 1)**oo == B(0, oo) + assert B(0, 1)**oo == B(0, oo) + assert B(2, 3)**oo is oo + assert B(1, 2)**oo == B(0, oo) + assert B(S.Half, 3)**oo == B(0, oo) + assert B(Rational(-1, 3), Rational(-1, 4))**oo is S.Zero + assert B(-1, Rational(-1, 2))**oo is S.NaN + assert B(-3, -2)**oo is zoo + assert B(-2, -1)**oo is S.NaN + assert B(-2, Rational(-1, 2))**oo is S.NaN + assert B(Rational(-1, 2), S.Half)**oo is S.Zero + assert B(Rational(-1, 2), 1)**oo == B(0, oo) + assert B(Rational(-2, 3), 2)**oo == B(0, oo) + assert B(-1, 1)**oo == B(-oo, oo) + assert B(-1, S.Half)**oo == B(-oo, oo) + assert B(-1, 2)**oo == B(-oo, oo) + assert B(-2, S.Half)**oo == B(-oo, oo) + + assert B(1, 2)**x == Pow(B(1, 2), x, evaluate=False) + + assert B(2, 3)**(-oo) is S.Zero + assert B(0, 2)**(-oo) == B(0, oo) + assert B(-1, 2)**(-oo) == B(-oo, oo) + + assert (tan(x)**sin(2*x)).subs(x, B(0, pi/2)) == \ + Pow(B(-oo, oo), B(0, 1)) + + +def test_AccumBounds_exponent(): + # base is 0 + z = 0**B(a, a + S.Half) + assert z.subs(a, 0) == B(0, 1) + assert z.subs(a, 1) == 0 + p = z.subs(a, -1) + assert p.is_Pow and p.args == (0, B(-1, -S.Half)) + # base > 0 + # when base is 1 the type of bounds does not matter + assert 1**B(a, a + 1) == 1 + # otherwise we need to know if 0 is in the bounds + assert S.Half**B(-2, 2) == B(S(1)/4, 4) + assert 2**B(-2, 2) == B(S(1)/4, 4) + + # +eps may introduce +oo + # if there is a negative integer exponent + assert B(0, 1)**B(S(1)/2, 1) == B(0, 1) + assert B(0, 1)**B(0, 1) == B(0, 1) + + # positive bases have positive bounds + assert B(2, 3)**B(-3, -2) == B(S(1)/27, S(1)/4) + assert B(2, 3)**B(-3, 2) == B(S(1)/27, 9) + + # bounds generating imaginary parts unevaluated + unchanged(Pow, B(-1, 1), B(1, 2)) + assert B(0, S(1)/2)**B(1, oo) == B(0, S(1)/2) + assert B(0, 1)**B(1, oo) == B(0, oo) + assert B(0, 2)**B(1, oo) == B(0, oo) + assert B(0, oo)**B(1, oo) == B(0, oo) + assert B(S(1)/2, 1)**B(1, oo) == B(0, oo) + assert B(S(1)/2, 1)**B(-oo, -1) == B(0, oo) + assert B(S(1)/2, 1)**B(-oo, oo) == B(0, oo) + assert B(S(1)/2, 2)**B(1, oo) == B(0, oo) + assert B(S(1)/2, 2)**B(-oo, -1) == B(0, oo) + assert B(S(1)/2, 2)**B(-oo, oo) == B(0, oo) + assert B(S(1)/2, oo)**B(1, oo) == B(0, oo) + assert B(S(1)/2, oo)**B(-oo, -1) == B(0, oo) + assert B(S(1)/2, oo)**B(-oo, oo) == B(0, oo) + assert B(1, 2)**B(1, oo) == B(0, oo) + assert B(1, 2)**B(-oo, -1) == B(0, oo) + assert B(1, 2)**B(-oo, oo) == B(0, oo) + assert B(1, oo)**B(1, oo) == B(0, oo) + assert B(1, oo)**B(-oo, -1) == B(0, oo) + assert B(1, oo)**B(-oo, oo) == B(0, oo) + assert B(2, oo)**B(1, oo) == B(2, oo) + assert B(2, oo)**B(-oo, -1) == B(0, S(1)/2) + assert B(2, oo)**B(-oo, oo) == B(0, oo) def test_comparison_AccumBounds(): - assert (AccumBounds(1, 3) < 4) == S.true - assert (AccumBounds(1, 3) < -1) == S.false - assert (AccumBounds(1, 3) < 2).rel_op == '<' - assert (AccumBounds(1, 3) <= 2).rel_op == '<=' + assert (B(1, 3) < 4) == S.true + assert (B(1, 3) < -1) == S.false + assert (B(1, 3) < 2).rel_op == '<' + assert (B(1, 3) <= 2).rel_op == '<=' - assert (AccumBounds(1, 3) > 4) == S.false - assert (AccumBounds(1, 3) > -1) == S.true - assert (AccumBounds(1, 3) > 2).rel_op == '>' - assert (AccumBounds(1, 3) >= 2).rel_op == '>=' + assert (B(1, 3) > 4) == S.false + assert (B(1, 3) > -1) == S.true + assert (B(1, 3) > 2).rel_op == '>' + assert (B(1, 3) >= 2).rel_op == '>=' - assert (AccumBounds(1, 3) < AccumBounds(4, 6)) == S.true - assert (AccumBounds(1, 3) < AccumBounds(2, 4)).rel_op == '<' - assert (AccumBounds(1, 3) < AccumBounds(-2, 0)) == S.false + assert (B(1, 3) < B(4, 6)) == S.true + assert (B(1, 3) < B(2, 4)).rel_op == '<' + assert (B(1, 3) < B(-2, 0)) == S.false - assert (AccumBounds(1, 3) <= AccumBounds(4, 6)) == S.true - assert (AccumBounds(1, 3) <= AccumBounds(-2, 0)) == S.false + assert (B(1, 3) <= B(4, 6)) == S.true + assert (B(1, 3) <= B(-2, 0)) == S.false - assert (AccumBounds(1, 3) > AccumBounds(4, 6)) == S.false - assert (AccumBounds(1, 3) > AccumBounds(-2, 0)) == S.true + assert (B(1, 3) > B(4, 6)) == S.false + assert (B(1, 3) > B(-2, 0)) == S.true - assert (AccumBounds(1, 3) >= AccumBounds(4, 6)) == S.false - assert (AccumBounds(1, 3) >= AccumBounds(-2, 0)) == S.true + assert (B(1, 3) >= B(4, 6)) == S.false + assert (B(1, 3) >= B(-2, 0)) == S.true # issue 13499 - assert (cos(x) > 0).subs(x, oo) == (AccumBounds(-1, 1) > 0) + assert (cos(x) > 0).subs(x, oo) == (B(-1, 1) > 0) c = Symbol('c') - raises(TypeError, lambda: (AccumBounds(0, 1) < c)) - raises(TypeError, lambda: (AccumBounds(0, 1) <= c)) - raises(TypeError, lambda: (AccumBounds(0, 1) > c)) - raises(TypeError, lambda: (AccumBounds(0, 1) >= c)) + raises(TypeError, lambda: (B(0, 1) < c)) + raises(TypeError, lambda: (B(0, 1) <= c)) + raises(TypeError, lambda: (B(0, 1) > c)) + raises(TypeError, lambda: (B(0, 1) >= c)) def test_contains_AccumBounds(): - assert (1 in AccumBounds(1, 2)) == S.true - raises(TypeError, lambda: a in AccumBounds(1, 2)) - assert 0 in AccumBounds(-1, 0) + assert (1 in B(1, 2)) == S.true + raises(TypeError, lambda: a in B(1, 2)) + assert 0 in B(-1, 0) raises(TypeError, lambda: - (cos(1)**2 + sin(1)**2 - 1) in AccumBounds(-1, 0)) - assert (-oo in AccumBounds(1, oo)) == S.true - assert (oo in AccumBounds(-oo, 0)) == S.true + (cos(1)**2 + sin(1)**2 - 1) in B(-1, 0)) + assert (-oo in B(1, oo)) == S.true + assert (oo in B(-oo, 0)) == S.true # issue 13159 - assert Mul(0, AccumBounds(-1, 1)) == Mul(AccumBounds(-1, 1), 0) == 0 + assert Mul(0, B(-1, 1)) == Mul(B(-1, 1), 0) == 0 import itertools - for perm in itertools.permutations([0, AccumBounds(-1, 1), x]): + for perm in itertools.permutations([0, B(-1, 1), x]): assert Mul(*perm) == 0 def test_intersection_AccumBounds(): - assert AccumBounds(0, 3).intersection(AccumBounds(1, 2)) == AccumBounds(1, 2) - assert AccumBounds(0, 3).intersection(AccumBounds(1, 4)) == AccumBounds(1, 3) - assert AccumBounds(0, 3).intersection(AccumBounds(-1, 2)) == AccumBounds(0, 2) - assert AccumBounds(0, 3).intersection(AccumBounds(-1, 4)) == AccumBounds(0, 3) - assert AccumBounds(0, 1).intersection(AccumBounds(2, 3)) == S.EmptySet - raises(TypeError, lambda: AccumBounds(0, 3).intersection(1)) + assert B(0, 3).intersection(B(1, 2)) == B(1, 2) + assert B(0, 3).intersection(B(1, 4)) == B(1, 3) + assert B(0, 3).intersection(B(-1, 2)) == B(0, 2) + assert B(0, 3).intersection(B(-1, 4)) == B(0, 3) + assert B(0, 1).intersection(B(2, 3)) == S.EmptySet + raises(TypeError, lambda: B(0, 3).intersection(1)) def test_union_AccumBounds(): - assert AccumBounds(0, 3).union(AccumBounds(1, 2)) == AccumBounds(0, 3) - assert AccumBounds(0, 3).union(AccumBounds(1, 4)) == AccumBounds(0, 4) - assert AccumBounds(0, 3).union(AccumBounds(-1, 2)) == AccumBounds(-1, 3) - assert AccumBounds(0, 3).union(AccumBounds(-1, 4)) == AccumBounds(-1, 4) - raises(TypeError, lambda: AccumBounds(0, 3).union(1)) + assert B(0, 3).union(B(1, 2)) == B(0, 3) + assert B(0, 3).union(B(1, 4)) == B(0, 4) + assert B(0, 3).union(B(-1, 2)) == B(-1, 3) + assert B(0, 3).union(B(-1, 4)) == B(-1, 4) + raises(TypeError, lambda: B(0, 3).union(1)) def test_issue_16469(): diff --git a/sympy/calculus/util.py b/sympy/calculus/util.py index 15b8042aa1e5..0c63094a6ac0 100644 --- a/sympy/calculus/util.py +++ b/sympy/calculus/util.py @@ -4,7 +4,7 @@ from sympy.core.compatibility import iterable from sympy.core.expr import AtomicExpr, Expr from sympy.core.function import expand_mul -from sympy.core.numbers import _sympifyit, oo +from sympy.core.numbers import _sympifyit, oo, zoo from sympy.core.relational import is_le, is_lt, is_ge, is_gt from sympy.core.sympify import _sympify from sympy.functions.elementary.miscellaneous import Min, Max @@ -894,22 +894,21 @@ class AccumulationBounds(AtomicExpr): `X * Y = \{ x*y \mid x \in X \cap y \in Y\}` - There is, however, no consensus on Interval division. - - `X / Y = \{ z \mid \exists x \in X, y \in Y \mid y \neq 0, z = x/y\}` - - Note: According to this definition the quotient of two AccumulationBounds - may not be a AccumulationBounds object but rather a union of - AccumulationBounds. - - Note - ==== - - The main focus in the interval arithmetic is on the simplest way to - calculate upper and lower endpoints for the range of values of a - function in one or more variables. These barriers are not necessarily - the supremum or infimum, since the precise calculation of those values - can be difficult or impossible. + When an AccumBounds is raised to a negative power, if 0 is contained + between the bounds then an infinite range is returned, otherwise if an + endpoint is 0 then a semi-infinite range with consistent sign will be returned. + + AccumBounds in expressions behave a lot like Intervals but the + semantics are not necessarily the same. Division (or exponentiation + to a negative integer power) could be handled with *intervals* by + returning a union of the results obtained after splitting the + bounds between negatives and positives, but that is not done with + AccumBounds. In addition, bounds are assumed to be independent of + each other; if the same bound is used in more than one place in an + expression, the result may not be the supremum or infimum of the + expression (see below). Finally, when a boundary is ``1``, + exponentiation to the power of ``oo`` yields ``oo``, neither + ``1`` nor ``nan``. Examples ======== @@ -936,21 +935,36 @@ class AccumulationBounds(AtomicExpr): `X^n = \{ x^n \mid x \in X\}` - otherwise + >>> AccumBounds(1, 4)**(S(1)/2) + AccumBounds(1, 2) + + otherwise, an infinite or semi-infinite result is obtained: - `X^n = \{ x^n \mid x \neq 0, x \in X\} \cup \{-\infty, \infty\}` + >>> 1/AccumBounds(-1, 1) + AccumBounds(-oo, oo) + >>> 1/AccumBounds(0, 2) + AccumBounds(1/2, oo) + >>> 1/AccumBounds(-oo, 0) + AccumBounds(-oo, 0) - Here for fractional `n`, the part of `X` resulting in a complex - AccumulationBounds object is neglected. + A boundary of 1 will always generate all nonnegatives: - >>> AccumBounds(-1, 4)**(S(1)/2) - AccumBounds(0, 2) + >>> AccumBounds(1, 2)**oo + AccumBounds(0, oo) + >>> AccumBounds(0, 1)**oo + AccumBounds(0, oo) - >>> AccumBounds(1, 2)**2 - AccumBounds(1, 4) + If the exponent is itself an AccumulationBounds or is not an + integer then unevaluated results will be returned unless the base + values are positive: - >>> AccumBounds(-1, oo)**(-1) - AccumBounds(-oo, oo) + >>> AccumBounds(2, 3)**AccumBounds(-1, 2) + AccumBounds(1/3, 9) + >>> AccumBounds(-2, 3)**AccumBounds(-1, 2) + AccumBounds(-2, 3)**AccumBounds(-1, 2) + + >>> AccumBounds(-2, -1)**(S(1)/2) + sqrt(AccumBounds(-2, -1)) Note: `^2` is not same as `*` @@ -980,8 +994,9 @@ class AccumulationBounds(AtomicExpr): object. But it doesn't necessarily evaluate the AccumulationBounds for that expression. - Same expression can be evaluated to different values depending upon - the form it is used for substitution. For example: + The same expression can be evaluated to different values depending upon + the form it is used for substitution since each instance of an + AccumulationBounds is considered independent. For example: >>> (x**2 + 2*x + 1).subs(x, AccumBounds(-1, 1)) AccumBounds(-1, 4) @@ -1014,15 +1029,18 @@ def __new__(cls, min, max): if not min.is_extended_real or not max.is_extended_real: raise ValueError("Only real AccumulationBounds are supported") - # Make sure that the created AccumBounds object will be valid. - if max.is_comparable and min.is_comparable: - if max < min: - raise ValueError( - "Lower limit should be smaller than upper limit") - if max == min: return max + # Make sure that the created AccumBounds object will be valid. + if max.is_number and min.is_number: + bad = max.is_comparable and min.is_comparable and max < min + else: + bad = (max - min).is_extended_negative + if bad: + raise ValueError( + "Lower limit should be smaller than upper limit") + return Basic.__new__(cls, min, max) # setting the operation priority @@ -1159,16 +1177,18 @@ def __rsub__(self, other): @_sympifyit('other', NotImplemented) def __mul__(self, other): + if self.args == (-oo, oo): + return self if isinstance(other, Expr): if isinstance(other, AccumBounds): - return AccumBounds(Min(Mul(self.min, other.min), - Mul(self.min, other.max), - Mul(self.max, other.min), - Mul(self.max, other.max)), - Max(Mul(self.min, other.min), - Mul(self.min, other.max), - Mul(self.max, other.min), - Mul(self.max, other.max))) + if other.args == (-oo, oo): + return other + v = set() + for i in self.args: + vi = other*i + for i in vi.args or (vi,): + v.add(i) + return AccumBounds(Min(*v), Max(*v)) if other is S.Infinity: if self.min.is_zero: return AccumBounds(0, oo) @@ -1181,8 +1201,6 @@ def __mul__(self, other): return AccumBounds(0, oo) if other.is_extended_real: if other.is_zero: - if self == AccumBounds(-oo, oo): - return AccumBounds(-oo, oo) if self.max is S.Infinity: return AccumBounds(0, oo) if self.min is S.NegativeInfinity: @@ -1223,9 +1241,9 @@ def __truediv__(self, other): if other.max.is_zero: return AccumBounds(self.max / other.min, oo) if other.max.is_extended_positive: - # the actual answer is a Union of AccumBounds, - # Union(AccumBounds(-oo, self.max/other.max), - # AccumBounds(self.max/other.min, oo)) + # if we were dealing with intervals we would return + # Union(Interval(-oo, self.max/other.max), + # Interval(self.max/other.min, oo)) return AccumBounds(-oo, oo) if other.min.is_zero and other.max.is_extended_positive: @@ -1236,9 +1254,9 @@ def __truediv__(self, other): if other.max.is_zero: return AccumBounds(-oo, self.min / other.min) if other.max.is_extended_positive: - # the actual answer is a Union of AccumBounds, - # Union(AccumBounds(-oo, self.min/other.min), - # AccumBounds(self.min/other.max, oo)) + # if we were dealing with intervals we would return + # Union(Interval(-oo, self.min/other.min), + # Interval(self.min/other.max, oo)) return AccumBounds(-oo, oo) if other.min.is_zero and other.max.is_extended_positive: @@ -1290,7 +1308,6 @@ def __rtruediv__(self, other): @_sympifyit('other', NotImplemented) def __pow__(self, other): - from sympy.functions.elementary.miscellaneous import real_root if isinstance(other, Expr): if other is S.Infinity: if self.min.is_extended_nonnegative: @@ -1303,8 +1320,8 @@ def __pow__(self, other): if self.min > -1: return S.Zero if self.max < -1: - return FiniteSet(-oo, oo) - return AccumBounds(-oo, oo) + return zoo + return S.NaN else: if self.min > -1: if self.max < 1: @@ -1313,55 +1330,96 @@ def __pow__(self, other): return AccumBounds(-oo, oo) if other is S.NegativeInfinity: - return (1 / self)**oo + return (1/self)**oo - if other.is_extended_real and other.is_number: - if other.is_zero: - return S.One - - if other.is_Integer: - if self.min.is_extended_positive: - return AccumBounds( - Min(self.min ** other, self.max ** other), - Max(self.min ** other, self.max ** other)) - elif self.max.is_extended_negative: - return AccumBounds( - Min(self.max ** other, self.min ** other), - Max(self.max ** other, self.min ** other)) - - if other % 2 == 0: - if other.is_extended_negative: - if self.min.is_zero: - return AccumBounds(self.max**other, oo) - if self.max.is_zero: - return AccumBounds(self.min**other, oo) - return AccumBounds(0, oo) - return AccumBounds( - S.Zero, Max(self.min**other, self.max**other)) - else: - if other.is_extended_negative: - if self.min.is_zero: - return AccumBounds(self.max**other, oo) - if self.max.is_zero: - return AccumBounds(-oo, self.min**other) - return AccumBounds(-oo, oo) - return AccumBounds(self.min**other, self.max**other) + # generically true + if (self.max - self.min).is_nonnegative: + # well defined + if self.min.is_nonnegative: + # no 0 to worry about + if other.is_nonnegative: + # no infinity to worry about + return self.func(self.min**other, self.max**other) + + if other.is_zero: + return S.One # x**0 = 1 + if other.is_Integer or other.is_integer: + if self.min.is_extended_positive: + return AccumBounds( + Min(self.min**other, self.max**other), + Max(self.min**other, self.max**other)) + elif self.max.is_extended_negative: + return AccumBounds( + Min(self.max**other, self.min**other), + Max(self.max**other, self.min**other)) + + if other % 2 == 0: + if other.is_extended_negative: + if self.min.is_zero: + return AccumBounds(self.max**other, oo) + if self.max.is_zero: + return AccumBounds(self.min**other, oo) + return AccumBounds(0, oo) + return AccumBounds( + S.Zero, Max(self.min**other, self.max**other)) + elif other % 2 == 1: + if other.is_extended_negative: + if self.min.is_zero: + return AccumBounds(self.max**other, oo) + if self.max.is_zero: + return AccumBounds(-oo, self.min**other) + return AccumBounds(-oo, oo) + return AccumBounds(self.min**other, self.max**other) + + # non-integer exponent + # 0**neg or neg**frac yields complex + if (other.is_number or other.is_rational) and ( + self.min.is_extended_nonnegative or ( + other.is_extended_nonnegative and + self.min.is_extended_nonnegative)): num, den = other.as_numer_denom() - if num == S.One: - if den % 2 == 0: - if S.Zero in self: - if self.min.is_extended_negative: - return AccumBounds(0, real_root(self.max, den)) - return AccumBounds(real_root(self.min, den), - real_root(self.max, den)) - if den!=1: - num_pow = self**num - return num_pow**(1 / den) - return AccumBounds(-oo, oo) + if num is S.One: + return AccumBounds(*[i**(1/den) for i in self.args]) + + elif den is not S.One: # e.g. if other is not Float + return (self**num)**(1/den) # ok for non-negative base + + if isinstance(other, AccumBounds): + if (self.min.is_extended_positive or + self.min.is_extended_nonnegative and + other.min.is_extended_nonnegative): + p = [self**i for i in other.args] + if not any(i.is_Pow for i in p): + a = [j for i in p for j in i.args or (i,)] + try: + return self.func(min(a), max(a)) + except TypeError: # can't sort + pass + + return Pow(self, other, evaluate=False) return NotImplemented + @_sympifyit('other', NotImplemented) + def __rpow__(self, other): + if other.is_real and other.is_extended_nonnegative and ( + self.max - self.min).is_extended_positive: + if other is S.One: + return S.One + if other.is_extended_positive: + a, b = [other**i for i in self.args] + if min(a, b) != a: + a, b = b, a + return self.func(a, b) + if other.is_zero: + if self.min.is_zero: + return self.func(0, 1) + if self.min.is_extended_positive: + return S.Zero + + return Pow(other, self, evaluate=False) + def __abs__(self): if self.max.is_extended_negative: return self.__neg__() diff --git a/sympy/codegen/rewriting.py b/sympy/codegen/rewriting.py index 0899fb925e94..f086f1d55815 100644 --- a/sympy/codegen/rewriting.py +++ b/sympy/codegen/rewriting.py @@ -181,9 +181,10 @@ def __init__(self, func, func_m_1): self.func_m_1 = func_m_1 def _try_func_m_1(self, expr): - protected, old_new = expr.replace(self.func, lambda arg: Dummy(), map=True) + old_new = {} + protected = expr.replace(self.func, lambda arg: old_new.setdefault(arg, Dummy())) factored = protected.factor() - new_old = {v: k for k, v in old_new.items()} + new_old = {v: self.func(k) for k, v in old_new.items()} return factored.replace(_d - 1, lambda d: self.func_m_1(new_old[d].args[0])).xreplace(new_old) def __call__(self, e): diff --git a/sympy/concrete/tests/test_sums_products.py b/sympy/concrete/tests/test_sums_products.py index d4435828a629..4218842431a5 100644 --- a/sympy/concrete/tests/test_sums_products.py +++ b/sympy/concrete/tests/test_sums_products.py @@ -1327,11 +1327,15 @@ def test_issue_8016(): cos(pi*n/2)*gamma(m + 1)/gamma(n/2 + 1)/gamma(m - n/2 + 1) -@XFAIL def test_issue_14313(): assert Sum(S.Half**floor(n/2), (n, 1, oo)).is_convergent() +def test_issue_14563(): + # The assertion was failing due to no assumptions methods in Sums and Product + assert 1 % Sum(1, (x, 0, 1)) == 1 + + def test_issue_16735(): assert Sum(5**n/gamma(n+1), (n, 1, oo)).is_convergent() is S.true @@ -1443,13 +1447,9 @@ def test_summation_by_residues(): assert eval_sum_residue(1 / x**6, (x, S(1), oo)) == pi**6/945 assert eval_sum_residue(1 / (x**2 + 9), (x, -oo, oo)) == pi/(3*tanh(3*pi)) assert eval_sum_residue(1 / (x**2 + 1)**2, (x, -oo, oo)) == \ - -pi*(-1/(2*tanh(pi)) - I*(-I*pi/tanh(pi) + \ - I*pi*tanh(pi))/(4*tanh(pi)) + I*(-I*pi*tanh(pi) + \ - I*pi/tanh(pi))/(4*tanh(pi))) + -pi*(-pi/(2*tanh(pi)**2) - 1/(2*tanh(pi)) + pi/2) assert eval_sum_residue(x**2 / (x**2 + 1)**2, (x, -oo, oo)) == \ - -pi*(-1/(2*tanh(pi)) - I*(-I*pi*tanh(pi) + \ - I*pi/tanh(pi))/(4*tanh(pi)) + I*(-I*pi/tanh(pi) + \ - I*pi*tanh(pi))/(4*tanh(pi))) + -pi*(-pi/2 - 1/(2*tanh(pi)) + pi/(2*tanh(pi)**2)) assert eval_sum_residue(1 / (4*x**2 - 1), (x, -oo, oo)) == 0 assert eval_sum_residue(x**2 / (x**2 - S(1)/4)**2, (x, -oo, oo)) == pi**2/2 assert eval_sum_residue(1 / (4*x**2 - 1)**2, (x, -oo, oo)) == pi**2/8 @@ -1490,7 +1490,6 @@ def test_summation_by_residues(): assert eval_sum_residue((-1)**x / x**2, (x, S(2), oo)) == 1 - pi**2/12 -@XFAIL def test_summation_by_residues_failing(): x = Symbol('x') diff --git a/sympy/core/add.py b/sympy/core/add.py index ca9bcfc4a30d..1108f47f5e94 100644 --- a/sympy/core/add.py +++ b/sympy/core/add.py @@ -493,6 +493,7 @@ def _eval_power(self, e): c = [sign(i) if i in bigs else i/big for i in c] addpow = Add(*[c*m for c, m in zip(c, m)])**e return big**e*addpow + @cacheit def _eval_derivative(self, s): return self.func(*[a.diff(s) for a in self.args]) @@ -1004,18 +1005,18 @@ def _eval_as_leading_term(self, x, cdir=0): except TypeError: return expr - new_expr=new_expr.together() - if new_expr.is_Add: - new_expr = new_expr.simplify() - - if not new_expr: + is_zero = new_expr.is_zero + if is_zero is None: + new_expr = new_expr.trigsimp().cancel() + is_zero = new_expr.is_zero + if is_zero is True: # simple leading term analysis gave us cancelled terms but we have to send # back a term, so compute the leading term (via series) n0 = min.getn() res = Order(1) incr = S.One while res.is_Order: - res = old._eval_nseries(x, n=n0+incr, logx=None, cdir=cdir).cancel().powsimp().trigsimp() + res = old._eval_nseries(x, n=n0+incr, logx=None, cdir=cdir).cancel().trigsimp() incr *= 2 return res.as_leading_term(x, cdir=cdir) diff --git a/sympy/core/basic.py b/sympy/core/basic.py index 8a613c130c63..ffb9fac6a06b 100644 --- a/sympy/core/basic.py +++ b/sympy/core/basic.py @@ -509,6 +509,10 @@ def free_symbols(self): @property def expr_free_symbols(self): + from sympy.utilities.exceptions import SymPyDeprecationWarning + SymPyDeprecationWarning(feature="expr_free_symbols method", + issue=21494, + deprecated_since_version="1.9").warn() return set() def as_dummy(self): @@ -762,24 +766,17 @@ def subs(self, *args, **kwargs): """ Substitutes old for new in an expression after sympifying args. - This method is capable of complicated substitution by calling - ``self._eval_subs()`` method. - - ``args`` is either: - - two arguments, e.g. ``foo.subs(old, new)`` - - an ``Equation`` instance, e.g. ``foo.subs(Eqn(old, new))`` - - one iterable argument, e.g. ``foo.subs(iterable)``. The iterable may be - o an iterable container with ``(old, new)`` pairs. In this case the + `args` is either: + - two arguments, e.g. foo.subs(old, new) + - one iterable argument, e.g. foo.subs(iterable). The iterable may be + o an iterable container with (old, new) pairs. In this case the replacements are processed in the order given with successive patterns possibly affecting replacements already made. - o an iterable container with ``Equation`` instance - instances, e.g. ``Eqn(old, new)``. This is processed as an - iterable container of ``(old, new)`` pairs. o a dict or set whose key/value items correspond to old/new pairs. In this case the old/new pairs will be sorted by op count and in - case of a tie, by number of args and the ``default_sort_key``. - The resulting sorted list is then processed as an iterable - container (see previous). + case of a tie, by number of args and the default_sort_key. The + resulting sorted list is then processed as an iterable container + (see previous). If the keyword ``simultaneous`` is True, the subexpressions will not be evaluated until all the substitutions have been made. @@ -801,22 +798,10 @@ def subs(self, *args, **kwargs): >>> (x + y).subs(reversed(reps)) x**2 + 2 - ``Equation`` instance is treated as a container of ``(old, new)`` pair. - - >>> from sympy import Eqn - >>> eqn1, eqn2 = Eqn(x, pi), Eqn(y, 2) - >>> (1 + x*y).subs(eqn1) - pi*y + 1 - >>> (1 + x*y).subs([eqn1, eqn2]) - 1 + 2*pi - - Complicated substitution can be done. Here, ``x**4`` is identified as - the square of ``x**2`` and substituted. - >>> (x**2 + x**4).subs(x**2, y) y**2 + y - To replace only the ``x**2`` but not the ``x**4``, use ``xreplace``: + To replace only the x**2 but not the x**4, use xreplace: >>> (x**2 + x**4).xreplace({x**2: y}) x**4 + y @@ -883,19 +868,16 @@ def subs(self, *args, **kwargs): See Also ======== - replace: replacement capable of doing wildcard-like matching, parsing of match, and conditional replacements xreplace: exact node replacement in expr tree; also capable of using matching rules - sympy.core.evalf.EvalfMixin.evalf: calculates the given formula - to a desired level of precision + sympy.core.evalf.EvalfMixin.evalf: calculates the given formula to a desired level of precision """ from sympy.core.compatibility import _nodes, default_sort_key from sympy.core.containers import Dict from sympy.core.symbol import Dummy, Symbol - from sympy.equation import Equation from sympy.utilities.misc import filldedent unordered = False @@ -906,8 +888,6 @@ def subs(self, *args, **kwargs): elif isinstance(sequence, (Dict, Mapping)): unordered = True sequence = sequence.items() - elif isinstance(sequence, Equation): - sequence = [sequence.args] elif not iterable(sequence): raise ValueError(filldedent(""" When a single argument is passed to subs @@ -920,10 +900,6 @@ def subs(self, *args, **kwargs): sequence = list(sequence) for i, s in enumerate(sequence): - if isinstance(s, Equation): - # treat equation as a container - s = s.args - if isinstance(s[0], str): # when old is a string we prefer Symbol s = Symbol(s[0]), s[1] @@ -1928,6 +1904,17 @@ def _aresame(a, b): return True +def _ne(a, b): + # use this as a second test after `a != b` if you want to make + # sure that things are truly equal, e.g. + # a, b = 0.5, S.Half + # a !=b or _ne(a, b) -> True + from .numbers import Number + # 0.5 == S.Half + if isinstance(a, Number) and isinstance(b, Number): + return a.__class__ != b.__class__ + + def _atomic(e, recursive=False): """Return atom-like quantities as far as substitution is concerned: Derivatives, Functions and Symbols. Don't diff --git a/sympy/core/expr.py b/sympy/core/expr.py index 1c7ef70b6a73..95ed85f76985 100644 --- a/sympy/core/expr.py +++ b/sympy/core/expr.py @@ -2418,11 +2418,15 @@ def extract_additively(self, c): @property def expr_free_symbols(self): """ - Like ``free_symbols``, but returns the free symbols only if they are contained in an expression node. + Like ``free_symbols``, but returns the free symbols only if + they are contained in an expression node. Examples ======== + >>> from sympy.utilities.exceptions import SymPyDeprecationWarning + >>> import warnings + >>> warnings.simplefilter("ignore", SymPyDeprecationWarning) >>> from sympy.abc import x, y >>> (x + y).expr_free_symbols {x, y} @@ -2437,6 +2441,10 @@ def expr_free_symbols(self): >>> t.free_symbols {x, y} """ + from sympy.utilities.exceptions import SymPyDeprecationWarning + SymPyDeprecationWarning(feature="expr_free_symbols method", + issue=21494, + deprecated_since_version="1.9").warn() return {j for i in self.args for j in i.expr_free_symbols} def could_extract_minus_sign(self): @@ -3921,6 +3929,10 @@ def _eval_nseries(self, x, n, logx, cdir=0): @property def expr_free_symbols(self): + from sympy.utilities.exceptions import SymPyDeprecationWarning + SymPyDeprecationWarning(feature="expr_free_symbols method", + issue=21494, + deprecated_since_version="1.9").warn() return {self} diff --git a/sympy/core/function.py b/sympy/core/function.py index a82c8d625f57..69df976b8bb5 100644 --- a/sympy/core/function.py +++ b/sympy/core/function.py @@ -2349,6 +2349,10 @@ def free_symbols(self): @property def expr_free_symbols(self): + from sympy.utilities.exceptions import SymPyDeprecationWarning + SymPyDeprecationWarning(feature="expr_free_symbols method", + issue=21494, + deprecated_since_version="1.9").warn() return (self.expr.expr_free_symbols - set(self.variables) | set(self.point.expr_free_symbols)) @@ -2394,17 +2398,29 @@ def _eval_subs(self, old, new): def _eval_derivative(self, s): # Apply the chain rule of the derivative on the substitution variables: - val = Add.fromiter(p.diff(s) * Subs(self.expr.diff(v), self.variables, self.point).doit() for v, p in zip(self.variables, self.point)) - - # Check if there are free symbols in `self.expr`: - # First get the `expr_free_symbols`, which returns the free symbols - # that are directly contained in an expression node (i.e. stop - # searching if the node isn't an expression). At this point turn the - # expressions into `free_symbols` and check if there are common free - # symbols in `self.expr` and the deriving factor. - fs1 = {j for i in self.expr_free_symbols for j in i.free_symbols} - if len(fs1 & s.free_symbols) > 0: - val += Subs(self.expr.diff(s), self.variables, self.point).doit() + f = self.expr + vp = V, P = self.variables, self.point + val = Add.fromiter(p.diff(s)*Subs(f.diff(v), *vp).doit() + for v, p in zip(V, P)) + + # these are all the free symbols in the expr + efree = f.free_symbols + # some symbols like IndexedBase include themselves and args + # as free symbols + compound = {i for i in efree if len(i.free_symbols) > 1} + # hide them and see what independent free symbols remain + dums = {Dummy() for i in compound} + masked = f.xreplace(dict(zip(compound, dums))) + ifree = masked.free_symbols - dums + # include the compound symbols + free = ifree | compound + # remove the variables already handled + free -= set(V) + # add back any free symbols of remaining compound symbols + free |= {i for j in free & compound for i in j.free_symbols} + # if symbols of s are in free then there is more to do + if free & s.free_symbols: + val += Subs(f.diff(s), self.variables, self.point).doit() return val def _eval_nseries(self, x, n, logx, cdir=0): diff --git a/sympy/core/mul.py b/sympy/core/mul.py index db84a3a62ec8..df604abf6d96 100644 --- a/sympy/core/mul.py +++ b/sympy/core/mul.py @@ -2050,24 +2050,24 @@ def _keep_coeff(coeff, factors, clear=True, sign=False): >>> _keep_coeff(S(-1), x + y, sign=True) -(x + y) """ - if not coeff.is_Number: if factors.is_Number: factors, coeff = coeff, factors else: return coeff*factors + if factors is S.One: + return coeff if coeff is S.One: return factors elif coeff is S.NegativeOne and not sign: return -factors elif factors.is_Add: if not clear and coeff.is_Rational and coeff.q != 1: - q = S(coeff.q) - for i in factors.args: - c, t = i.as_coeff_Mul() - r = c/q - if r == int(r): - return coeff*factors + args = [i.as_coeff_Mul() for i in factors.args] + args = [(_keep_coeff(c, coeff), m) for c, m in args] + if any(c == int(c) for c, _ in args): + return Add._from_args([Mul._from_args( + i[1:] if i[0] == 1 else i) for i in args]) return Mul(coeff, factors, evaluate=False) elif factors.is_Mul: margs = list(factors.args) @@ -2079,8 +2079,10 @@ def _keep_coeff(coeff, factors, clear=True, sign=False): margs.insert(0, coeff) return Mul._from_args(margs) else: - return coeff*factors - + m = coeff*factors + if m.is_Number and not factors.is_Number: + m = Mul._from_args((coeff, factors)) + return m def expand_2arg(e): from sympy.simplify.simplify import bottom_up diff --git a/sympy/core/power.py b/sympy/core/power.py index 4c8aabc70247..2f398c69347e 100644 --- a/sympy/core/power.py +++ b/sympy/core/power.py @@ -789,7 +789,14 @@ def _eval_is_polar(self): return self.base.is_polar def _eval_subs(self, old, new): - from sympy import exp, log, Symbol + from sympy import exp, log, Symbol, AccumBounds + + if isinstance(self.exp, AccumBounds): + b = self.base.subs(old, new) + e = self.exp.subs(old, new) + if isinstance(e, AccumBounds): + return e.__rpow__(b) + return self.func(b, e) def _check(ct1, ct2, old): """Return (bool, pow, remainder_pow) where, if bool is True, then the @@ -1659,7 +1666,7 @@ def mul(d1, d2): else: raise NotImplementedError() if not d.is_positive: - g = (b - f).simplify()/f + g = g.simplify() _, d = g.leadterm(x) if not d.is_positive: raise NotImplementedError() @@ -1696,7 +1703,7 @@ def mul(d1, d2): if not (e.is_integer and e.is_positive and (e*d - n).is_nonpositive and res == _mexpand(self)): res += O(x**n, x) - return res + return _mexpand(res) def _eval_as_leading_term(self, x, cdir=0): from ..series import Order diff --git a/sympy/core/relational.py b/sympy/core/relational.py index bac0dafb3e4c..8c8188f9de7d 100644 --- a/sympy/core/relational.py +++ b/sympy/core/relational.py @@ -1410,7 +1410,7 @@ def split_real_imag(expr): # Compare e.g. zoo with 1+I*oo by comparing args arglhs = arg(lhs) argrhs = arg(rhs) - # Guard against Eq(nan, nan) -> Falsesymp + # Guard against Eq(nan, nan) -> False if not (arglhs == S.NaN and argrhs == S.NaN): return fuzzy_bool(is_eq(arglhs, argrhs, assumptions)) diff --git a/sympy/core/sympify.py b/sympy/core/sympify.py index ed5ba267d941..5dbf59b97b9d 100644 --- a/sympy/core/sympify.py +++ b/sympy/core/sympify.py @@ -435,14 +435,7 @@ def sympify(a, locals=None, convert_xor=True, strict=False, rational=False, if iterable(a): try: return type(a)([sympify(x, locals=locals, convert_xor=convert_xor, - rational=rational) for x in a]) - except TypeError: - # Not all iterables are rebuildable with their type. - pass - if isinstance(a, dict): - try: - return type(a)([sympify(x, locals=locals, convert_xor=convert_xor, - rational=rational) for x in a.items()]) + rational=rational, evaluate=evaluate) for x in a]) except TypeError: # Not all iterables are rebuildable with their type. pass diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py index 221113ddbc19..8eb52785ebd7 100644 --- a/sympy/core/tests/test_args.py +++ b/sympy/core/tests/test_args.py @@ -1890,16 +1890,6 @@ def test_sympy__core__symbol__Wild(): assert _test_args(Wild('x', exclude=[x])) -@SKIP("abstract class") -def test_sympy__equation__equation__SymbolicRelation(): - pass - - -def test_sympy__equation__equation__Equation(): - from sympy.equation import Equation - assert _test_args(Equation(1,1)) - - @SKIP("abstract class") def test_sympy__functions__combinatorial__factorials__CombinatorialFunction(): pass diff --git a/sympy/core/tests/test_basic.py b/sympy/core/tests/test_basic.py index 328b00a070fb..09947a0b1c10 100644 --- a/sympy/core/tests/test_basic.py +++ b/sympy/core/tests/test_basic.py @@ -15,7 +15,6 @@ from sympy.functions.elementary.exponential import exp from sympy.testing.pytest import raises from sympy.core import I, pi -from sympy.equation import Eqn b1 = Basic() b2 = Basic(b1) @@ -143,14 +142,6 @@ def test_subs_with_unicode_symbols(): assert replaced.name == 'x' -def test_subs_with_Eqn(): - assert b21.subs(Eqn(b2, b1)) == Basic(b1, b1) - assert b21.subs(Eqn(b2, b21)) == Basic(b21, b1) - assert b3.subs(Eqn(b2, b1)) == b2 - - assert b21.subs([Eqn(b2, b1), Eqn(b1, b2)]) == Basic(b2, b2) - - def test_atoms(): assert b21.atoms() == {Basic()} diff --git a/sympy/core/tests/test_expr.py b/sympy/core/tests/test_expr.py index e60a163da74c..c9762eb0c45e 100644 --- a/sympy/core/tests/test_expr.py +++ b/sympy/core/tests/test_expr.py @@ -441,6 +441,19 @@ def test_as_leading_term(): raises(ValueError, lambda: (x + 1).as_leading_term(1)) + # https://github.com/sympy/sympy/issues/21177 + f = -3*x + (x + Rational(3, 2) - sqrt(3)*S.ImaginaryUnit/2)**2\ + - Rational(3, 2) + 3*sqrt(3)*S.ImaginaryUnit/2 + assert f.as_leading_term(x) == \ + (3*sqrt(3)*x - 3*S.ImaginaryUnit*x)/(sqrt(3) + 3*S.ImaginaryUnit) + + # https://github.com/sympy/sympy/issues/21245 + f = 1 - x - x**2 + fi = (1 + sqrt(5))/2 + assert f.subs(x, y + 1/fi).as_leading_term(y) == \ + (-36*sqrt(5)*y - 80*y)/(16*sqrt(5) + 36) + + def test_leadterm2(): assert (x*cos(1)*cos(1 + sin(1)) + sin(1 + sin(1))).leadterm(x) == \ (sin(1 + sin(1)), 0) @@ -457,7 +470,7 @@ def test_as_leading_term2(): def test_as_leading_term3(): assert (2 + pi + x).as_leading_term(x) == 2 + pi - assert (2*x + pi*x + x**2).as_leading_term(x) == (2 + pi)*x + assert (2*x + pi*x + x**2).as_leading_term(x) == 2*x + pi*x def test_as_leading_term4(): @@ -2102,6 +2115,7 @@ def test_ExprBuilder(): eb.args.extend([x, x]) assert eb.build() == x**2 + def test_non_string_equality(): # Expressions should not compare equal to strings x = symbols('x') @@ -2122,3 +2136,9 @@ def __repr__(self): assert (x == BadRepr()) is False assert (x != BadRepr()) is True + + +def test_21494(): + from sympy.testing.pytest import warns_deprecated_sympy + with warns_deprecated_sympy(): + assert x.expr_free_symbols == {x} diff --git a/sympy/core/tests/test_power.py b/sympy/core/tests/test_power.py index 0721bfbfd038..d7b15a074e77 100644 --- a/sympy/core/tests/test_power.py +++ b/sympy/core/tests/test_power.py @@ -333,7 +333,7 @@ def test_issue_6990(): a = Symbol('a') b = Symbol('b') assert (sqrt(a + b*x + x**2)).series(x, 0, 3).removeO() == \ - sqrt(a)*x**2*(1/(2*a) - b**2/(8*a**2)) + sqrt(a) + b*x/(2*sqrt(a)) + sqrt(a) + x**2*(1/(2*sqrt(a)) - b**2/(8*a**(S(3)/2))) + b*x/(2*sqrt(a)) def test_issue_6068(): diff --git a/sympy/core/tests/test_sympify.py b/sympy/core/tests/test_sympify.py index 55bb1ddc6bf8..d88a76bf1d0a 100644 --- a/sympy/core/tests/test_sympify.py +++ b/sympy/core/tests/test_sympify.py @@ -769,3 +769,29 @@ def test_issue_14706(): raises(SympifyError, lambda: sympify(numpy.array([1]), strict=True)) raises(SympifyError, lambda: sympify(z1, strict=True)) raises(SympifyError, lambda: sympify(z2, strict=True)) + + +def test_issue_21536(): + #test to check evaluate=False in case of iterable input + u = sympify("x+3*x+2", evaluate=False) + v = sympify("2*x+4*x+2+4", evaluate=False) + + assert u.is_Add and set(u.args) == {x, 3*x, 2} + assert v.is_Add and set(v.args) == {2*x, 4*x, 2, 4} + assert sympify(["x+3*x+2", "2*x+4*x+2+4"], evaluate=False) == [u, v] + + #test to check evaluate=True in case of iterable input + u = sympify("x+3*x+2", evaluate=True) + v = sympify("2*x+4*x+2+4", evaluate=True) + + assert u.is_Add and set(u.args) == {4*x, 2} + assert v.is_Add and set(v.args) == {6*x, 6} + assert sympify(["x+3*x+2", "2*x+4*x+2+4"], evaluate=True) == [u, v] + + #test to check evaluate with no input in case of iterable input + u = sympify("x+3*x+2") + v = sympify("2*x+4*x+2+4") + + assert u.is_Add and set(u.args) == {4*x, 2} + assert v.is_Add and set(v.args) == {6*x, 6} + assert sympify(["x+3*x+2", "2*x+4*x+2+4"]) == [u, v] diff --git a/sympy/diffgeom/diffgeom.py b/sympy/diffgeom/diffgeom.py index ba82beebbab6..f4f1aa444d41 100644 --- a/sympy/diffgeom/diffgeom.py +++ b/sympy/diffgeom/diffgeom.py @@ -193,8 +193,9 @@ class CoordSystem(Basic): relations : dict, optional Key is a tuple of two strings, who are the names of the systems where - the coordinates transform from and transform to. Value is a tuple of - transformed coordinates. + the coordinates transform from and transform to. + Value is a tuple of the symbols before transformation and a tuple of + the expressions after transformation. Examples ======== @@ -209,8 +210,8 @@ class CoordSystem(Basic): >>> x, y = symbols('x y', real=True) >>> r, theta = symbols('r theta', nonnegative=True) >>> relation_dict = { - ... ('Car2D', 'Pol'): (sqrt(x**2 + y**2), atan2(y, x)), - ... ('Pol', 'Car2D'): (r*cos(theta), r*sin(theta)) + ... ('Car2D', 'Pol'): [(x, y), (sqrt(x**2 + y**2), atan2(y, x))], + ... ('Pol', 'Car2D'): [(r, theta), (r*cos(theta), r*sin(theta))] ... } >>> Car2D = CoordSystem('Car2D', p, (x, y), relation_dict) >>> Pol = CoordSystem('Pol', p, (r, theta), relation_dict) @@ -314,8 +315,12 @@ def __new__(cls, name, patch, symbols=None, relations={}, **kwargs): if not isinstance(s2, Str): s2 = Str(s2) key = Tuple(s1, s2) + + # Old version used Lambda as a value. if isinstance(v, Lambda): - v = tuple(v(*symbols)) + v = (tuple(v.signature), tuple(v.expr)) + else: + v = (tuple(v[0]), tuple(v[1])) rel_temp[key] = v relations = Dict(rel_temp) @@ -396,7 +401,7 @@ def transformation(self, sys): if self == sys: expr = Matrix(self.symbols) elif key in self.relations: - expr = Matrix(self.relations[key]) + expr = Matrix(self.relations[key][1]) elif key[::-1] in self.relations: expr = Matrix(self._inverse_transformation(sys, self)) else: @@ -404,44 +409,58 @@ def transformation(self, sys): return Lambda(signature, expr) @staticmethod - def _inverse_transformation(sys1, sys2): - # Find the transformation relation from sys2 to sys1 - forward_transform_expressions = sys1.transform(sys2) - - inv_results = solve( - [t[0] - t[1] for t in zip(sys2.symbols, forward_transform_expressions)], - list(sys1.symbols), dict=True) - if len(inv_results) == 0: - raise NotImplementedError( - "Cannot solve inverse of transformation from {} to {}".format(sys1, sys2)) - elif len(inv_results) > 1: - raise ValueError( - "Obtained multiple results for inverse of transformation from {} to {}".format(sys1, sys2) - ) + def _solve_inverse(sym1, sym2, exprs, sys1_name, sys2_name): + ret = solve( + [t[0] - t[1] for t in zip(sym2, exprs)], + list(sym1), dict=True) + + if len(ret) == 0: + temp = "Cannot solve inverse relation from {} to {}." + raise NotImplementedError(temp.format(sys1_name, sys2_name)) + elif len(ret) > 1: + temp = "Obtained multiple inverse relation from {} to {}." + raise ValueError(temp.format(sys1_name, sys2_name)) + + return ret[0] - inv_results = inv_results[0] + @classmethod + def _inverse_transformation(cls, sys1, sys2): + # Find the transformation relation from sys2 to sys1 + forward = sys1.transform(sys2) + inv_results = cls._solve_inverse(sys1.symbols, sys2.symbols, forward, + sys1.name, sys2.name) signature = tuple(sys1.symbols) return [inv_results[s] for s in signature] @classmethod @cacheit def _indirect_transformation(cls, sys1, sys2): - # Find the transformation relation between two indirectly connected coordinate systems + # Find the transformation relation between two indirectly connected + # coordinate systems + rel = sys1.relations path = cls._dijkstra(sys1, sys2) - Lambdas = [] - for i in range(len(path) - 1): - s1, s2 = path[i], path[i + 1] - Lambdas.append(s1.transformation(s2)) - syms = Lambdas[-1].signature - expr = syms - for l in reversed(Lambdas): - expr = l(*expr) - return Lambda(syms, expr) + + transforms = [] + for s1, s2 in zip(path, path[1:]): + if (s1, s2) in rel: + transforms.append(rel[(s1, s2)]) + else: + sym2, inv_exprs = rel[(s2, s1)] + sym1 = tuple(Dummy() for i in sym2) + ret = cls._solve_inverse(sym2, sym1, inv_exprs, s2, s1) + ret = tuple(ret[s] for s in sym2) + transforms.append((sym1, ret)) + syms = sys1.args[2] + exprs = syms + for newsyms, newexprs in transforms: + exprs = tuple(e.subs(zip(newsyms, exprs)) for e in newexprs) + return exprs @staticmethod def _dijkstra(sys1, sys2): # Use Dijkstra algorithm to find the shortest path between two indirectly-connected # coordinate systems + # return value is the list of the names of the systems. relations = sys1.relations graph = {} for s1, s2 in relations.keys(): @@ -465,7 +484,7 @@ def visit(sys): path_dict[newsys][1] = [i for i in path_dict[sys][1]] path_dict[newsys][1].append(sys) - visit(sys1) + visit(sys1.name) while True: min_distance = max(path_dict.values(), key=lambda x:x[0])[0] @@ -478,10 +497,10 @@ def visit(sys): break visit(newsys) - result = path_dict[sys2][1] - result.append(sys2) + result = path_dict[sys2.name][1] + result.append(sys2.name) - if result == [sys2]: + if result == [sys2.name]: raise KeyError("Two coordinate systems are not connected.") return result diff --git a/sympy/diffgeom/rn.py b/sympy/diffgeom/rn.py index 635e85f7365a..c6254734f4da 100644 --- a/sympy/diffgeom/rn.py +++ b/sympy/diffgeom/rn.py @@ -30,8 +30,8 @@ r, theta = symbols('rho theta', nonnegative=True) relations_2d = { - ('rectangular', 'polar'): (sqrt(x**2 + y**2), atan2(y, x)), - ('polar', 'rectangular'): (r*cos(theta), r*sin(theta)), + ('rectangular', 'polar'): [(x, y), (sqrt(x**2 + y**2), atan2(y, x))], + ('polar', 'rectangular'): [(r, theta), (r*cos(theta), r*sin(theta))], } R2_r = CoordSystem('rectangular', R2_origin, (x, y), relations_2d) # type: Any @@ -74,18 +74,24 @@ rho, psi, r, theta, phi = symbols('rho psi r theta phi', nonnegative=True) relations_3d = { - ('rectangular', 'cylindrical'): (sqrt(x**2 + y**2), atan2(y, x), z), - ('cylindrical', 'rectangular'): (rho*cos(psi), rho*sin(psi), z), - ('rectangular', 'spherical'): (sqrt(x**2 + y**2 + z**2), - acos(z/sqrt(x**2 + y**2 + z**2)), - atan2(y, x)), - ('spherical', 'rectangular'): (r*sin(theta)*cos(phi), - r*sin(theta)*sin(phi), - r*cos(theta)), - ('cylindrical', 'spherical'): (sqrt(rho**2 + z**2), - acos(z/sqrt(rho**2 + z**2)), - psi), - ('spherical', 'cylindrical'): (r*sin(theta), phi, r*cos(theta)), + ('rectangular', 'cylindrical'): [(x, y, z), + (sqrt(x**2 + y**2), atan2(y, x), z)], + ('cylindrical', 'rectangular'): [(rho, psi, z), + (rho*cos(psi), rho*sin(psi), z)], + ('rectangular', 'spherical'): [(x, y, z), + (sqrt(x**2 + y**2 + z**2), + acos(z/sqrt(x**2 + y**2 + z**2)), + atan2(y, x))], + ('spherical', 'rectangular'): [(r, theta, phi), + (r*sin(theta)*cos(phi), + r*sin(theta)*sin(phi), + r*cos(theta))], + ('cylindrical', 'spherical'): [(rho, psi, z), + (sqrt(rho**2 + z**2), + acos(z/sqrt(rho**2 + z**2)), + psi)], + ('spherical', 'cylindrical'): [(r, theta, phi), + (r*sin(theta), phi, r*cos(theta))], } R3_r = CoordSystem('rectangular', R3_origin, (x, y, z), relations_3d) # type: Any diff --git a/sympy/diffgeom/tests/test_diffgeom.py b/sympy/diffgeom/tests/test_diffgeom.py index f5c1209702d4..dc36547e87eb 100644 --- a/sympy/diffgeom/tests/test_diffgeom.py +++ b/sympy/diffgeom/tests/test_diffgeom.py @@ -1,10 +1,10 @@ +from sympy.core import Lambda, Symbol, symbols from sympy.diffgeom.rn import R2, R2_p, R2_r, R3_r, R3_c, R3_s, R2_origin from sympy.diffgeom import (CoordSystem, Commutator, Differential, TensorProduct, WedgeProduct, BaseCovarDerivativeOp, CovarDerivativeOp, LieDerivative, covariant_order, contravariant_order, twoform_to_matrix, metric_to_Christoffel_1st, metric_to_Christoffel_2nd, metric_to_Riemann_components, metric_to_Ricci_components, intcurve_diffequ, intcurve_series) -from sympy.core import Symbol, symbols from sympy.simplify import trigsimp, simplify from sympy.functions import sqrt, atan2, sin from sympy.matrices import Matrix @@ -14,6 +14,78 @@ TP = TensorProduct +def test_coordsys_transform(): + # test inverse transforms + p, q, r, s = symbols('p q r s') + rel = {('first', 'second'): [(p, q), (q, -p)]} + R2_pq = CoordSystem('first', R2_origin, [p, q], rel) + R2_rs = CoordSystem('second', R2_origin, [r, s], rel) + r, s = R2_rs.symbols + assert R2_rs.transform(R2_pq) == Matrix([[-s], [r]]) + + # inverse transform impossible case + a, b = symbols('a b', positive=True) + rel = {('first', 'second'): [(a,), (-a,)]} + R2_a = CoordSystem('first', R2_origin, [a], rel) + R2_b = CoordSystem('second', R2_origin, [b], rel) + # This transformation is uninvertible because there is no positive a, b satisfying a = -b + with raises(NotImplementedError): + R2_b.transform(R2_a) + + # inverse transform ambiguous case + c, d = symbols('c d') + rel = {('first', 'second'): [(c,), (c**2,)]} + R2_c = CoordSystem('first', R2_origin, [c], rel) + R2_d = CoordSystem('second', R2_origin, [d], rel) + # The transform method should throw if it finds multiple inverses for a coordinate transformation. + with raises(ValueError): + R2_d.transform(R2_c) + + # test indirect transformation + a, b, c, d, e, f = symbols('a, b, c, d, e, f') + rel = {('C1', 'C2'): [(a, b), (2*a, 3*b)], + ('C2', 'C3'): [(c, d), (3*c, 2*d)]} + C1 = CoordSystem('C1', R2_origin, (a, b), rel) + C2 = CoordSystem('C2', R2_origin, (c, d), rel) + C3 = CoordSystem('C3', R2_origin, (e, f), rel) + a, b = C1.symbols + c, d = C2.symbols + e, f = C3.symbols + assert C2.transform(C1) == Matrix([c/2, d/3]) + assert C1.transform(C3) == Matrix([6*a, 6*b]) + assert C3.transform(C1) == Matrix([e/6, f/6]) + assert C3.transform(C2) == Matrix([e/3, f/2]) + + a, b, c, d, e, f = symbols('a, b, c, d, e, f') + rel = {('C1', 'C2'): [(a, b), (2*a, 3*b + 1)], + ('C3', 'C2'): [(e, f), (-e - 2, 2*f)]} + C1 = CoordSystem('C1', R2_origin, (a, b), rel) + C2 = CoordSystem('C2', R2_origin, (c, d), rel) + C3 = CoordSystem('C3', R2_origin, (e, f), rel) + a, b = C1.symbols + c, d = C2.symbols + e, f = C3.symbols + assert C2.transform(C1) == Matrix([c/2, (d - 1)/3]) + assert C1.transform(C3) == Matrix([-2*a - 2, (3*b + 1)/2]) + assert C3.transform(C1) == Matrix([-e/2 - 1, (2*f - 1)/3]) + assert C3.transform(C2) == Matrix([-e - 2, 2*f]) + + # old signature uses Lambda + a, b, c, d, e, f = symbols('a, b, c, d, e, f') + rel = {('C1', 'C2'): Lambda((a, b), (2*a, 3*b + 1)), + ('C3', 'C2'): Lambda((e, f), (-e - 2, 2*f))} + C1 = CoordSystem('C1', R2_origin, (a, b), rel) + C2 = CoordSystem('C2', R2_origin, (c, d), rel) + C3 = CoordSystem('C3', R2_origin, (e, f), rel) + a, b = C1.symbols + c, d = C2.symbols + e, f = C3.symbols + assert C2.transform(C1) == Matrix([c/2, (d - 1)/3]) + assert C1.transform(C3) == Matrix([-2*a - 2, (3*b + 1)/2]) + assert C3.transform(C1) == Matrix([-e/2 - 1, (2*f - 1)/3]) + assert C3.transform(C2) == Matrix([-e - 2, 2*f]) + + def test_R2(): x0, y0, r0, theta0 = symbols('x0, y0, r0, theta0', real=True) point_r = R2_r.point([x0, y0]) @@ -37,45 +109,6 @@ def test_R2(): R2_r, R2_r.coord_tuple_transform_to(R2_p, m)).applyfunc(simplify) -def test_inverse_transformations(): - p, q, r, s = symbols('p q r s') - - relations_quarter_rotation = { - ('first', 'second'): (q, -p) - } - - R2_pq = CoordSystem('first', R2_origin, [p, q], relations_quarter_rotation) - R2_rs = CoordSystem('second', R2_origin, [r, s], relations_quarter_rotation) - - # The transform method should derive the inverse transformation if not given explicitly - assert R2_rs.transform(R2_pq) == Matrix([[-R2_rs.symbols[1]], [R2_rs.symbols[0]]]) - - a, b = symbols('a b', positive=True) - relations_uninvertible_transformation = { - ('first', 'second'): (-a,) - } - - R2_a = CoordSystem('first', R2_origin, [a], relations_uninvertible_transformation) - R2_b = CoordSystem('second', R2_origin, [b], relations_uninvertible_transformation) - - # The transform method should throw if it cannot invert the coordinate transformation. - # This transformation is uninvertible because there is no positive a, b satisfying a = -b - with raises(NotImplementedError): - R2_b.transform(R2_a) - - c, d = symbols('c d') - relations_ambiguous_inverse = { - ('first', 'second'): (c**2,) - } - - R2_c = CoordSystem('first', R2_origin, [c], relations_ambiguous_inverse) - R2_d = CoordSystem('second', R2_origin, [d], relations_ambiguous_inverse) - - # The transform method should throw if it finds multiple inverses for a coordinate transformation. - with raises(ValueError): - R2_d.transform(R2_c) - - def test_R3(): a, b, c = symbols('a b c', positive=True) m = Matrix([[a], [b], [c]]) diff --git a/sympy/equation/__init__.py b/sympy/equation/__init__.py deleted file mode 100644 index edf4087a0428..000000000000 --- a/sympy/equation/__init__.py +++ /dev/null @@ -1,7 +0,0 @@ -""" -Module for symbolic relations. -""" - -from .equation import Equation, Eqn - -__all__ = ["Equation", "Eqn",] diff --git a/sympy/equation/equation.py b/sympy/equation/equation.py deleted file mode 100644 index f5e5f58c46e7..000000000000 --- a/sympy/equation/equation.py +++ /dev/null @@ -1,117 +0,0 @@ -from sympy.core import Basic -from sympy.core.sympify import _sympify - -from .sideproxy import SideProxy - - -class SymbolicRelation(Basic): - """ - Base class for all symbolic binary relations. - - Explanation - =========== - - Unlike boolean relation, symbolic relation behaves as a container for - the arguments. Its truth value is never evaluated, and all features are - aimed for symbolic manipulation of the arguments. - - See the docstring of :obj:`~.Equation` for the examples. - - See Also - ======== - - sympy.core.relational.Relational : Boolean relation - - """ - - def __new__(cls, lhs, rhs, **kwargs): - lhs = _sympify(lhs) - rhs = _sympify(rhs) - return super().__new__(cls, lhs, rhs) - - @property - def lhs(self): - """The left-hand side of the relation.""" - return self.args[0] - - @property - def rhs(self): - """The right-hand side of the relation.""" - return self.args[1] - - @property - def apply(self): - """Proxy object to apply operation on both sides.""" - return SideProxy(self, "both") - - @property - def applylhs(self): - """Proxy object to apply operation on left hand sides.""" - return SideProxy(self, "lhs") - - @property - def applyrhs(self): - """Proxy object to apply operation on right hand sides.""" - return SideProxy(self, "rhs") - - -class Equation(SymbolicRelation): - """ - Symbolic equation. - - Examples - ======== - - Symbolic equation is not reduced to boolean value. - - >>> from sympy import Eqn - >>> Eqn(1, 1).simplify() - Eqn(1, 1) - - Arguments can be manipulated by ``applylhs``, ``applyrhs``, or ``apply`` - properties. Note that this is purely structural manipulation, and does not - guarantee mathematical correctness. - - >>> from sympy import cos, sin, gamma, trigsimp - >>> from sympy.abc import x - >>> eqn = Eqn(sin(x)**2 + cos(x)**2, gamma(x)/gamma(x-2)) - >>> eqn.apply.simplify() # apply simplify method on both sides - Eqn(1, (x - 2)*(x - 1)) - >>> eqn.applyrhs.simplify() # apply simplify method on right hand side - Eqn(sin(x)**2 + cos(x)**2, (x - 2)*(x - 1)) - >>> eqn.applylhs(trigsimp) # apply trigsimp function on left hand side - Eqn(1, gamma(x)/gamma(x - 2)) - - Equation is a valid argument for ``subs()`` method. - - >>> (x**2).subs(Eqn(x, 2)) - 4 - - See Also - ======== - - sympy.core.relational.Equality : Boolean equality - - """ - rel_op = "==" - - @property - def reversed(self): - """ - Return the equation with sides reversed. - - Examples - ======== - - >>> from sympy import Eqn - >>> from sympy.abc import x - >>> Eqn(x, 1) - Eqn(x, 1) - >>> _.reversed - Eqn(1, x) - - """ - return self.func(self.rhs, self.lhs) - - -Eqn = Equation diff --git a/sympy/equation/sideproxy.py b/sympy/equation/sideproxy.py deleted file mode 100644 index 293bcfe31af7..000000000000 --- a/sympy/equation/sideproxy.py +++ /dev/null @@ -1,120 +0,0 @@ -""" -Module to implement proxy object for symbolic manipulation of relations. -""" - -from functools import partial - - -class SideProxy: - """ - Proxy object to apply methods on relation. - - This object passes the expression manipulation down to the arguments - of applied relation object, and returns a new one. It is implemented - to distinguish the operation on relation itself, and operation on - the arguments. - - Applying attributes on proxy object calls the attributes of relation's - arguments and assembles a new relation. Calling the proxy object and - passing the function as first argument applies the function on - relation's arguments and assembles a new one. - - This class is intended to be accessed by ``.apply``, ``.applylhs``, - and ``.applyrhs`` properties of ``SymbolicRelation`` object. Do not - construct this object directly. - - Parameters - ========== - - rel : SymbolicRelation - - side : "both", "lhs" or "rhs" - - Examples - ======== - - ``SideProxy`` can be accessed via ``apply`` attribute of - ``SymbolicRelation`` object. - - >>> from sympy import exp, sin, Function, Eqn - >>> from sympy.abc import x - >>> f = Function('f') - >>> eqn = Eqn(f(x), sin(x)) - >>> eqn.apply - - >>> eqn.applyrhs - - - Applying methods on ``SideProxy`` applies the operation down to the - corresponding argument(s). - - >>> eqn.apply.diff(x) - Eqn(Derivative(f(x), x), cos(x)) - >>> eqn.applyrhs.rewrite(exp) - Eqn(f(x), -I*(exp(I*x) - exp(-I*x))/2) - - Passing a function or a class to ``SideProxy`` applies it down to the - corresponding argument(s). Note that keyword arguments can be passed, but - the function or the class must be unary. - - >>> Eqn(x, 2).apply(lambda x: x**2) - Eqn(x**2, 4) - >>> Eqn(x, 0).apply(sin, evaluate=False) - Eqn(sin(x), sin(0)) - """ - def __init__(self, rel, side): - if side not in ("both", "lhs", "rhs"): - raise ValueError("Only 'both', 'lhs' or 'rhs' is allowed for side.") - self.rel = rel - self.side = side - - def __repr__(self): - return "" % (self.rel, self.side) - - def _apply_func(self, func, *args, **kwargs): - lhs, rhs = self.rel.args - if self.side in ("both", "lhs"): - lhs = func(lhs, *args, **kwargs) - if self.side in ("both", "rhs"): - rhs = func(rhs, *args, **kwargs) - return self.rel.func(lhs, rhs) - - def _apply_attr(self, attrname): - lhs, rhs = self.rel.args - if self.side in ("both", "lhs"): - lhs = getattr(lhs, attrname) - if self.side in ("both", "rhs"): - rhs = getattr(rhs, attrname) - return self.rel.func(lhs, rhs) - - def _apply_method(self, methodname, *args, **kwargs): - lhs, rhs = self.rel.args - if self.side in ("both", "lhs"): - lhs = getattr(lhs, methodname)(*args, **kwargs) - if self.side in ("both", "rhs"): - rhs = getattr(rhs, methodname)(*args, **kwargs) - return self.rel.func(lhs, rhs) - - def __getattr__(self, attrname): - lhs, rhs = self.rel.args - - attrs = [] - if self.side in ("both", "lhs"): - attrs.append(getattr(lhs, attrname)) - if self.side in ("both", "rhs"): - attrs.append(getattr(rhs, attrname)) - - if not all(callable(attr) for attr in attrs): - return partial(self._apply_attr, attrname) - elif all(callable(attr) for attr in attrs): - return partial(self._apply_method, attrname) - else: - raise AttributeError("Inconsistent methods are called.") - - def __call__(self, func, **kwargs): - lhs, rhs = self.rel.args - if self.side in ("both", "lhs"): - lhs = func(lhs, **kwargs) - if self.side in ("both", "rhs"): - rhs = func(rhs, **kwargs) - return self.rel.func(lhs, rhs) diff --git a/sympy/equation/tests/__init__.py b/sympy/equation/tests/__init__.py deleted file mode 100644 index e69de29bb2d1..000000000000 diff --git a/sympy/equation/tests/test_equation.py b/sympy/equation/tests/test_equation.py deleted file mode 100644 index e21c971e3c54..000000000000 --- a/sympy/equation/tests/test_equation.py +++ /dev/null @@ -1,9 +0,0 @@ -from sympy import Eqn, Equation, symbols - -x, y = symbols('x y') - - -def test_Eqn(): - assert Eqn(x, y) == Equation(x, y) - assert Eqn(x, y).lhs == x - assert Eqn(x, y).rhs == y diff --git a/sympy/functions/elementary/hyperbolic.py b/sympy/functions/elementary/hyperbolic.py index 1ad392579c9d..892e19cabb90 100644 --- a/sympy/functions/elementary/hyperbolic.py +++ b/sympy/functions/elementary/hyperbolic.py @@ -10,8 +10,6 @@ from sympy.core.logic import fuzzy_or, fuzzy_and - - def _rewrite_hyperbolics_as_exp(expr): expr = sympify(expr) return expr.xreplace({h: h.rewrite(exp) @@ -877,6 +875,27 @@ def _eval_as_leading_term(self, x, cdir=0): else: return self.func(arg) + def _eval_expand_trig(self, **hints): + arg = self.args[0] + if arg.is_Add: + from sympy import symmetric_poly + CX = [coth(x, evaluate=False)._eval_expand_trig() for x in arg.args] + p = [[], []] + n = len(arg.args) + for i in range(n, -1, -1): + p[(n - i) % 2].append(symmetric_poly(i, CX)) + return Add(*p[0])/Add(*p[1]) + elif arg.is_Mul: + from sympy import binomial + coeff, x = arg.as_coeff_Mul(rational=True) + if coeff.is_Integer and coeff > 1: + c = coth(x, evaluate=False) + p = [[], []] + for i in range(coeff, -1, -1): + p[(coeff - i) % 2].append(binomial(coeff, i)*c**i) + return Add(*p[0])/Add(*p[1]) + return coth(arg) + class ReciprocalHyperbolicFunction(HyperbolicFunction): """Base class for reciprocal functions of hyperbolic functions. """ @@ -939,6 +958,9 @@ def _eval_expand_complex(self, deep=True, **hints): re_part, im_part = self.as_real_imag(deep=True, **hints) return re_part + S.ImaginaryUnit*im_part + def _eval_expand_trig(self, **hints): + return self._calculate_reciprocal("_eval_expand_trig", **hints) + def _eval_as_leading_term(self, x, cdir=0): return (1/self._reciprocal_of(self.args[0]))._eval_as_leading_term(x) diff --git a/sympy/functions/elementary/miscellaneous.py b/sympy/functions/elementary/miscellaneous.py index 85b606e754bc..b458bd2af7c6 100644 --- a/sympy/functions/elementary/miscellaneous.py +++ b/sympy/functions/elementary/miscellaneous.py @@ -1,4 +1,5 @@ from sympy.core import Function, S, sympify +from sympy.utilities.iterables import sift from sympy.core.add import Add from sympy.core.containers import Tuple from sympy.core.compatibility import ordered @@ -432,7 +433,6 @@ def _collapse_arguments(cls, args, **assumptions): >>> Min(a, Max(b, Min(c, d, Max(a, e)))) Min(a, Max(b, Min(a, c, d))) - """ from sympy.utilities.iterables import ordered from sympy.simplify.simplify import walk @@ -518,28 +518,33 @@ def do(ai, a): # easy case where all functions contain something in common; # trying to find some optimal subset of args to modify takes # too long + + def factor_minmax(args): + is_other = lambda arg: isinstance(arg, other) + other_args, remaining_args = sift(args, is_other, binary=True) + if not other_args: + return args + + # Min(Max(x, y, z), Max(x, y, u, v)) -> {x,y}, ({z}, {u,v}) + arg_sets = [set(arg.args) for arg in other_args] + common = set.intersection(*arg_sets) + if not common: + return args + + new_other_args = list(common) + arg_sets_diff = [arg_set - common for arg_set in arg_sets] + + # If any set is empty after removing common then all can be + # discarded e.g. Min(Max(a, b, c), Max(a, b)) -> Max(a, b) + if all(arg_sets_diff): + other_args_diff = [other(*s, evaluate=False) for s in arg_sets_diff] + new_other_args.append(cls(*other_args_diff, evaluate=False)) + + other_args_factored = other(*new_other_args, evaluate=False) + return remaining_args + [other_args_factored] + if len(args) > 1: - common = None - remove = [] - sets = [] - for i in range(len(args)): - a = args[i] - if not isinstance(a, other): - continue - s = set(a.args) - common = s if common is None else (common & s) - if not common: - break - sets.append(s) - remove.append(i) - if common: - sets = filter(None, [s - common for s in sets]) - sets = [other(*s, evaluate=False) for s in sets] - for i in reversed(remove): - args.pop(i) - oargs = [cls(*sets)] if sets else [] - oargs.extend(common) - args.append(other(*oargs, evaluate=False)) + args = factor_minmax(args) return args diff --git a/sympy/functions/elementary/piecewise.py b/sympy/functions/elementary/piecewise.py index 91bfceafaa50..303b94e4beb6 100644 --- a/sympy/functions/elementary/piecewise.py +++ b/sympy/functions/elementary/piecewise.py @@ -590,7 +590,12 @@ def handler(ipw): elif hi.is_Symbol: pos = pos.xreplace({hi: lo + p}).xreplace({p: hi - lo}) neg = neg.xreplace({hi: lo - p}).xreplace({p: lo - hi}) - + # evaluate limits that may have unevaluate Min/Max + touch = lambda _: _.replace( + lambda x: isinstance(x, (Min, Max)), + lambda x: x.func(*x.args)) + neg = touch(neg) + pos = touch(pos) # assemble return expression; make the first condition be Lt # b/c then the first expression will look the same whether # the lo or hi limit is symbolic diff --git a/sympy/functions/elementary/tests/test_hyperbolic.py b/sympy/functions/elementary/tests/test_hyperbolic.py index 582edd42e78c..24c2fae25dd4 100644 --- a/sympy/functions/elementary/tests/test_hyperbolic.py +++ b/sympy/functions/elementary/tests/test_hyperbolic.py @@ -1,7 +1,7 @@ from sympy import (symbols, Symbol, sinh, nan, oo, zoo, pi, asinh, acosh, log, sqrt, coth, I, cot, E, tanh, tan, cosh, cos, S, sin, Rational, atanh, acoth, Integer, O, exp, sech, sec, csch, asech, acsch, acos, asin, expand_mul, - AccumBounds, im, re) + AccumBounds, im, re, expand_trig) from sympy.core.expr import unchanged from sympy.core.function import ArgumentIndexError @@ -327,6 +327,11 @@ def test_coth(): x = Symbol('x', extended_real=True) assert coth(x).as_real_imag(deep=False) == (coth(x), 0) + assert expand_trig(coth(2*x)) == (coth(x)**2 + 1)/(2*coth(x)) + assert expand_trig(coth(3*x)) == (coth(x)**3 + 3*coth(x))/(1 + 3*coth(x)**2) + + assert expand_trig(coth(x + y)) == (1 + coth(x)*coth(y))/(coth(x) + coth(y)) + def test_coth_series(): x = Symbol('x') @@ -397,6 +402,8 @@ def test_csch(): assert csch(n).is_real is True + assert expand_trig(csch(x + y)) == 1/(sinh(x)*cosh(y) + cosh(x)*sinh(y)) + def test_csch_series(): x = Symbol('x') @@ -465,6 +472,8 @@ def test_sech(): assert sech(n).is_real is True + assert expand_trig(sech(x + y)) == 1/(cosh(x)*cosh(y) + sinh(x)*sinh(y)) + def test_sech_series(): x = Symbol('x') diff --git a/sympy/functions/elementary/tests/test_miscellaneous.py b/sympy/functions/elementary/tests/test_miscellaneous.py index 0e1a5b1d08ce..bc6f2e2762ee 100644 --- a/sympy/functions/elementary/tests/test_miscellaneous.py +++ b/sympy/functions/elementary/tests/test_miscellaneous.py @@ -410,6 +410,10 @@ def test_issue_12638(): assert Min(a, b, Max(a, b, c)) == Min(a, b) assert Min(a, b, Max(a, c)) == Min(a, b) +def test_issue_21399(): + from sympy.abc import a, b, c + assert Max(Min(a, b), Min(a, b, c)) == Min(a, b) + def test_instantiation_evaluation(): from sympy.abc import v, w, x, y, z diff --git a/sympy/functions/elementary/tests/test_piecewise.py b/sympy/functions/elementary/tests/test_piecewise.py index 1833939f5ee2..935bf613c567 100644 --- a/sympy/functions/elementary/tests/test_piecewise.py +++ b/sympy/functions/elementary/tests/test_piecewise.py @@ -233,10 +233,8 @@ def test_piecewise_integrate1ca(): assert g.integrate((x, 1, -2)) == g1y.subs(y, -2) assert g.integrate((x, 0, 1)) == gy1.subs(y, 0) assert g.integrate((x, 1, 0)) == g1y.subs(y, 0) - # XXX Make test pass without simplify - assert g.integrate((x, 2, 1)) == gy1.subs(y, 2).simplify() - assert g.integrate((x, 1, 2)) == g1y.subs(y, 2).simplify() - + assert g.integrate((x, 2, 1)) == gy1.subs(y, 2) + assert g.integrate((x, 1, 2)) == g1y.subs(y, 2) assert piecewise_fold(gy1.rewrite(Piecewise)) == \ Piecewise( (1, y <= -1), @@ -249,22 +247,19 @@ def test_piecewise_integrate1ca(): (y**2/2 + y - S.Half, y <= 0), (-y**2/2 + y - S.Half, y < 1), (0, True)) - - # g1y and gy1 should simplify if the condition that y < 1 - # is applied, e.g. Min(1, Max(-1, y)) --> Max(-1, y) - # XXX Make test pass without simplify - assert gy1.simplify() == Piecewise( + assert gy1 == Piecewise( ( -Min(1, Max(-1, y))**2/2 - Min(1, Max(-1, y)) + Min(1, Max(0, y))**2 + S.Half, y < 1), (0, True) ) - assert g1y.simplify() == Piecewise( + assert g1y == Piecewise( ( Min(1, Max(-1, y))**2/2 + Min(1, Max(-1, y)) - Min(1, Max(0, y))**2 - S.Half, y < 1), (0, True)) + @slow def test_piecewise_integrate1cb(): y = symbols('y', real=True) diff --git a/sympy/functions/elementary/tests/test_trigonometric.py b/sympy/functions/elementary/tests/test_trigonometric.py index 61b63a319899..07e44ab75c39 100644 --- a/sympy/functions/elementary/tests/test_trigonometric.py +++ b/sympy/functions/elementary/tests/test_trigonometric.py @@ -3,7 +3,7 @@ cosh, atan2, exp, log, asinh, acoth, atanh, O, cancel, Matrix, re, im, Float, Pow, gcd, sec, csc, cot, diff, simplify, Heaviside, arg, conjugate, series, FiniteSet, asec, acsc, Mul, sinc, jn, - AccumBounds, Interval, ImageSet, Lambda, besselj, Add) + AccumBounds, Interval, ImageSet, Lambda, besselj, Add, limit) from sympy.core.expr import unchanged from sympy.core.function import ArgumentIndexError from sympy.core.relational import Ne, Eq @@ -196,6 +196,17 @@ def test_sin_rewrite(): assert sin(cos(x)).rewrite(Pow) == sin(cos(x)) +def _test_extrig(f, i, e): + from sympy.core.expr import unchanged + from sympy.core.function import expand_trig + assert unchanged(f, i) + assert expand_trig(f(i)) == f(i) + # testing directly instead of with .expand(trig=True) + # because the other expansions undo the unevaluated Mul + assert expand_trig(f(Mul(i, 1, evaluate=False))) == e + assert abs(f(i) - e).n() < 1e-10 + + def test_sin_expansion(): # Note: these formulas are not unique. The ones here come from the # Chebyshev formulas. @@ -205,8 +216,8 @@ def test_sin_expansion(): assert sin(2*x).expand(trig=True) == 2*sin(x)*cos(x) assert sin(3*x).expand(trig=True) == -4*sin(x)**3 + 3*sin(x) assert sin(4*x).expand(trig=True) == -8*sin(x)**3*cos(x) + 4*sin(x)*cos(x) - assert sin(2).expand(trig=True) == 2*sin(1)*cos(1) - assert sin(3).expand(trig=True) == -4*sin(1)**3 + 3*sin(1) + _test_extrig(sin, 2, 2*sin(1)*cos(1)) + _test_extrig(sin, 3, -4*sin(1)**3 + 3*sin(1)) def test_sin_AccumBounds(): @@ -410,8 +421,8 @@ def test_cos_expansion(): assert cos(2*x).expand(trig=True) == 2*cos(x)**2 - 1 assert cos(3*x).expand(trig=True) == 4*cos(x)**3 - 3*cos(x) assert cos(4*x).expand(trig=True) == 8*cos(x)**4 - 8*cos(x)**2 + 1 - assert cos(2).expand(trig=True) == 2*cos(1)**2 - 1 - assert cos(3).expand(trig=True) == 4*cos(1)**3 - 3*cos(1) + _test_extrig(cos, 2, 2*cos(1)**2 - 1) + _test_extrig(cos, 3, 4*cos(1)**3 - 3*cos(1)) def test_cos_AccumBounds(): @@ -529,6 +540,10 @@ def test_tan(): assert tan(r).is_finite is None assert tan(I*r).is_finite is True + # https://github.com/sympy/sympy/issues/21177 + f = tan(pi*(x + S(3)/2))/(3*x) + assert f.as_leading_term(x) == -1/(3*pi*x**2) + def test_tan_series(): assert tan(x).series(x, 0, 9) == \ @@ -585,6 +600,8 @@ def test_tan_expansion(): assert 0 == tan(2*x).expand(trig=True).rewrite(tan).subs([(tan(x), Rational(1, 7))])*24 - 7 assert 0 == tan(3*x).expand(trig=True).rewrite(tan).subs([(tan(x), Rational(1, 5))])*55 - 37 assert 0 == tan(4*x - pi/4).expand(trig=True).rewrite(tan).subs([(tan(x), Rational(1, 5))])*239 - 1 + _test_extrig(tan, 2, 2*tan(1)/(1 - tan(1)**2)) + _test_extrig(tan, 3, (-tan(1)**3 + 3*tan(1))/(1 - 3*tan(1)**2)) def test_tan_AccumBounds(): @@ -684,6 +701,10 @@ def test_cot(): assert cot(x).subs(x, 3*pi) is zoo + # https://github.com/sympy/sympy/issues/21177 + f = cot(pi*(x + 4))/(3*x) + assert f.as_leading_term(x) == 1/(3*pi*x**2) + def test_tan_cot_sin_cos_evalf(): assert abs((tan(pi*Rational(8, 15))*cos(pi*Rational(8, 15))/sin(pi*Rational(8, 15)) - 1).evalf()) < 1e-14 @@ -745,15 +766,23 @@ def test_cot_subs(): def test_cot_expansion(): - assert cot(x + y).expand(trig=True) == ((cot(x)*cot(y) - 1)/(cot(x) + cot(y))).expand() - assert cot(x - y).expand(trig=True) == (-(cot(x)*cot(y) + 1)/(cot(x) - cot(y))).expand() - assert cot(x + y + z).expand(trig=True) == ( + assert cot(x + y).expand(trig=True).together() == ( + (cot(x)*cot(y) - 1)/(cot(x) + cot(y))) + assert cot(x - y).expand(trig=True).together() == ( + cot(x)*cot(-y) - 1)/(cot(x) + cot(-y)) + assert cot(x + y + z).expand(trig=True).together() == ( (cot(x)*cot(y)*cot(z) - cot(x) - cot(y) - cot(z))/ - (-1 + cot(x)*cot(y) + cot(x)*cot(z) + cot(y)*cot(z))).expand() - assert cot(3*x).expand(trig=True) == ((cot(x)**3 - 3*cot(x))/(3*cot(x)**2 - 1)).expand() - assert 0 == cot(2*x).expand(trig=True).rewrite(cot).subs([(cot(x), Rational(1, 3))])*3 + 4 - assert 0 == cot(3*x).expand(trig=True).rewrite(cot).subs([(cot(x), Rational(1, 5))])*55 - 37 - assert 0 == cot(4*x - pi/4).expand(trig=True).rewrite(cot).subs([(cot(x), Rational(1, 7))])*863 + 191 + (-1 + cot(x)*cot(y) + cot(x)*cot(z) + cot(y)*cot(z))) + assert cot(3*x).expand(trig=True).together() == ( + (cot(x)**2 - 3)*cot(x)/(3*cot(x)**2 - 1)) + assert cot(2*x).expand(trig=True) == cot(x)/2 - 1/(2*cot(x)) + assert cot(3*x).expand(trig=True).together() == ( + cot(x)**2 - 3)*cot(x)/(3*cot(x)**2 - 1) + assert cot(4*x - pi/4).expand(trig=True).cancel() == ( + -tan(x)**4 + 4*tan(x)**3 + 6*tan(x)**2 - 4*tan(x) - 1 + )/(tan(x)**4 + 4*tan(x)**3 - 6*tan(x)**2 - 4*tan(x) + 1) + _test_extrig(cot, 2, (-1 + cot(1)**2)/(2*cot(1))) + _test_extrig(cot, 3, (-3*cot(1) + cot(1)**3)/(-1 + 3*cot(1)**2)) def test_cot_AccumBounds(): @@ -787,11 +816,20 @@ def test_sinc(): assert sinc(-x) == sinc(x) - assert sinc(x).diff() == Piecewise(((x*cos(x) - sin(x)) / x**2, Ne(x, 0)), (0, True)) - - assert sinc(x).diff(x).equals(sinc(x).rewrite(sin).diff(x)) - - assert sinc(x).diff().subs(x, 0) is S.Zero + assert sinc(x).diff(x) == cos(x)/x - sin(x)/x**2 + assert sinc(x).diff(x) == (sin(x)/x).diff(x) + assert sinc(x).diff(x, x) == (-sin(x) - 2*cos(x)/x + 2*sin(x)/x**2)/x + assert sinc(x).diff(x, x) == (sin(x)/x).diff(x, x) + assert limit(sinc(x).diff(x), x, 0) == 0 + assert limit(sinc(x).diff(x, x), x, 0) == -S(1)/3 + + # https://github.com/sympy/sympy/issues/11402 + # + # assert sinc(x).diff(x) == Piecewise(((x*cos(x) - sin(x)) / x**2, Ne(x, 0)), (0, True)) + # + # assert sinc(x).diff(x).equals(sinc(x).rewrite(sin).diff(x)) + # + # assert sinc(x).diff(x).subs(x, 0) is S.Zero assert sinc(x).series() == 1 - x**2/6 + x**4/120 + O(x**6) diff --git a/sympy/functions/elementary/trigonometric.py b/sympy/functions/elementary/trigonometric.py index 2146fd263a60..5be402678c11 100644 --- a/sympy/functions/elementary/trigonometric.py +++ b/sympy/functions/elementary/trigonometric.py @@ -460,7 +460,7 @@ def _eval_expand_trig(self, **hints): cx = cos(x, evaluate=False)._eval_expand_trig() cy = cos(y, evaluate=False)._eval_expand_trig() return sx*cy + sy*cx - else: + elif arg.is_Mul: n, x = arg.as_coeff_Mul(rational=True) if n.is_Integer: # n will be positive because of .eval # canonicalization @@ -918,7 +918,7 @@ def _eval_expand_trig(self, **hints): cx = cos(x, evaluate=False)._eval_expand_trig() cy = cos(y, evaluate=False)._eval_expand_trig() return cx*cy - sx*sy - else: + elif arg.is_Mul: coeff, terms = arg.as_coeff_Mul(rational=True) if coeff.is_Integer: return chebyshevt(coeff, cos(terms)) @@ -1202,7 +1202,7 @@ def _eval_expand_trig(self, **hints): p[1 - i % 2] += symmetric_poly(i, Y)*(-1)**((i % 4)//2) return (p[0]/p[1]).subs(list(zip(Y, TX))) - else: + elif arg.is_Mul: coeff, terms = arg.as_coeff_Mul(rational=True) if coeff.is_Integer and coeff > 1: I = S.ImaginaryUnit @@ -1254,14 +1254,12 @@ def _eval_rewrite_as_sqrt(self, arg, **kwargs): def _eval_as_leading_term(self, x, cdir=0): arg = self.args[0] - x0 = arg.subs(x, 0) - n = x0/S.Pi + x0 = arg.subs(x, 0).cancel() + n = 2*x0/S.Pi if n.is_integer: - lt = (arg - n*S.Pi).as_leading_term(x) + lt = (arg - n*S.Pi/2).as_leading_term(x) return lt if n.is_even else -1/lt - if not x0.is_finite: - return self - return self.func(x0) + return self.func(x0) if x0.is_finite else self def _eval_is_extended_real(self): # FIXME: currently tan(pi/2) return zoo @@ -1537,13 +1535,13 @@ def _eval_rewrite_as_sqrt(self, arg, **kwargs): return y def _eval_as_leading_term(self, x, cdir=0): - from sympy import Order - arg = self.args[0].as_leading_term(x) - - if x in arg.free_symbols and Order(1, x).contains(arg): - return 1/arg - else: - return self.func(arg) + arg = self.args[0] + x0 = arg.subs(x, 0).cancel() + n = 2*x0/S.Pi + if n.is_integer: + lt = (arg - n*S.Pi/2).as_leading_term(x) + return 1/lt if n.is_even else -lt + return self.func(x0) if x0.is_finite else self def _eval_is_extended_real(self): return self.args[0].is_extended_real @@ -1567,14 +1565,14 @@ def _eval_expand_trig(self, **hints): for i in range(n, -1, -1): p[(n - i) % 2] += symmetric_poly(i, Y)*(-1)**(((n - i) % 4)//2) return (p[0]/p[1]).subs(list(zip(Y, CX))) - else: + elif arg.is_Mul: coeff, terms = arg.as_coeff_Mul(rational=True) if coeff.is_Integer and coeff > 1: I = S.ImaginaryUnit z = Symbol('dummy', real=True) P = ((z + I)**coeff).expand() return (re(P)/im(P)).subs([(z, cot(terms))]) - return cot(arg) + return cot(arg) # XXX sec and csc return 1/cos and 1/sin def _eval_is_finite(self): arg = self.args[0] @@ -1933,7 +1931,7 @@ class sinc(Function): * Differentiation >>> sinc(x).diff() - Piecewise(((x*cos(x) - sin(x))/x**2, Ne(x, 0)), (0, True)) + cos(x)/x - sin(x)/x**2 * Series Expansion @@ -1961,7 +1959,13 @@ class sinc(Function): def fdiff(self, argindex=1): x = self.args[0] if argindex == 1: - return Piecewise(((x*cos(x) - sin(x))/x**2, Ne(x, S.Zero)), (S.Zero, S.true)) + # We would like to return the Piecewise here, but Piecewise.diff + # currently can't handle removable singularities, meaning things + # like sinc(x).diff(x, 2) give the wrong answer at x = 0. See + # https://github.com/sympy/sympy/issues/11402. + # + # return Piecewise(((x*cos(x) - sin(x))/x**2, Ne(x, S.Zero)), (S.Zero, S.true)) + return cos(x)/x - sin(x)/x**2 else: raise ArgumentIndexError(self, argindex) diff --git a/sympy/functions/special/tests/test_hyper.py b/sympy/functions/special/tests/test_hyper.py index eb14631aab8a..3d3dc5067fe9 100644 --- a/sympy/functions/special/tests/test_hyper.py +++ b/sympy/functions/special/tests/test_hyper.py @@ -345,6 +345,10 @@ def test_limits(): assert limit(meijerg((), (), (1,), (0,), -x), x, 0) == \ meijerg(((), ()), ((1,), (0,)), 0) # issue 6052 + # https://github.com/sympy/sympy/issues/11465 + assert limit(1/hyper((1, ), (1, ), x), x, 0) == 1 + + def test_appellf1(): a, b1, b2, c, x, y = symbols('a b1 b2 c x y') assert appellf1(a, b2, b1, c, y, x) == appellf1(a, b1, b2, c, x, y) diff --git a/sympy/integrals/tests/test_transforms.py b/sympy/integrals/tests/test_transforms.py index 69d6f450fb81..218ee7ec7576 100644 --- a/sympy/integrals/tests/test_transforms.py +++ b/sympy/integrals/tests/test_transforms.py @@ -451,7 +451,7 @@ def mysimp(expr): @slow def test_laplace_transform(): - from sympy import fresnels, fresnelc + from sympy import fresnels, fresnelc, DiracDelta LT = laplace_transform a, b, c, = symbols('a b c', positive=True) t = symbols('t') @@ -505,6 +505,28 @@ def test_laplace_transform(): ((s - 1)/((s - 1)**2 + 1), -oo), ] + # DiracDelta function: standard cases + assert LT(DiracDelta(t), t, s) == (1, -oo, True) + assert LT(DiracDelta(a*t), t, s) == (1/a, -oo, True) + assert LT(DiracDelta(t/42), t, s) == (42, -oo, True) + assert LT(DiracDelta(t+42), t, s) == (0, -oo, True) + assert LT(DiracDelta(t)+DiracDelta(t-42), t, s) == \ + (1 + exp(-42*s), -oo, True) + assert LT(DiracDelta(t)-a*exp(-a*t), t, s) == (-a/(a + s) + 1, 0, True) + assert LT(exp(-t)*(DiracDelta(t)+DiracDelta(t-42)), t, s) == \ + (exp(-42*s - 42) + 1, -oo, True) + # Collection of cases that cannot be fully evaluated and/or would catch + # some common implementation errors + assert LT(DiracDelta(t**2), t, s) == LaplaceTransform(DiracDelta(t**2), t, s) + assert LT(DiracDelta(t**2 - 1), t, s) == (exp(-s)/2, -oo, True) + assert LT(DiracDelta(t*(1 - t)), t, s) == \ + LaplaceTransform(DiracDelta(-t**2 + t), t, s) + assert LT((DiracDelta(t) + 1)*(DiracDelta(t - 1) + 1), t, s) == \ + (LaplaceTransform(DiracDelta(t)*DiracDelta(t - 1), t, s) + \ + 1 + exp(-s) + 1/s, 0, True) + assert LT(DiracDelta(2*t - 2*exp(a)), t, s) == \ + (exp(-s*exp(a))/2, -oo, True) + # Fresnel functions assert laplace_transform(fresnels(t), t, s) == \ ((-sin(s**2/(2*pi))*fresnels(s/pi) + sin(s**2/(2*pi))/2 - @@ -552,7 +574,8 @@ def test_issue_8368_7173(): def test_inverse_laplace_transform(): - from sympy import sinh, cosh, besselj, besseli, simplify, factor_terms + from sympy import sinh, cosh, besselj, besseli, simplify, factor_terms,\ + DiracDelta ILT = inverse_laplace_transform a, b, c, = symbols('a b c', positive=True) t = symbols('t') @@ -560,28 +583,39 @@ def test_inverse_laplace_transform(): def simp_hyp(expr): return factor_terms(expand_mul(expr)).rewrite(sin) - # just test inverses of all of the above + assert ILT(1, s, t) == DiracDelta(t) assert ILT(1/s, s, t) == Heaviside(t) + assert ILT(a/(a + s), s, t) == a*exp(-a*t)*Heaviside(t) + assert ILT(s/(a + s), s, t) == -a*exp(-a*t)*Heaviside(t) + DiracDelta(t) + assert ILT((a + s)**(-2), s, t) == t*exp(-a*t)*Heaviside(t) + assert ILT((a + s)**(-5), s, t) == t**4*exp(-a*t)*Heaviside(t)/24 + assert ILT(a/(a**2 + s**2), s, t) == sin(a*t)*Heaviside(t) + assert ILT(s/(s**2 + a**2), s, t) == cos(a*t)*Heaviside(t) + assert ILT(b/(b**2 + (a + s)**2), s, t) == exp(-a*t)*sin(b*t)*Heaviside(t) + assert ILT(b*s/(b**2 + (a + s)**2), s, t) +\ + (a*sin(b*t) - b*cos(b*t))*exp(-a*t)*Heaviside(t) == 0 + assert ILT(exp(-a*s)/s, s, t) == Heaviside(-a + t) + assert ILT(exp(-a*s)/(b + s), s, t) == exp(b*(a - t))*Heaviside(-a + t) + assert ILT((b + s)/(a**2 + (b + s)**2), s, t) == \ + exp(-b*t)*cos(a*t)*Heaviside(t) + assert ILT(exp(-a*s)/s**b, s, t) == \ + (-a + t)**(b - 1)*Heaviside(-a + t)/gamma(b) + assert ILT(exp(-a*s)/sqrt(s**2 + 1), s, t) == \ + Heaviside(-a + t)*besselj(0, a - t) + assert ILT(1/(s*sqrt(s + 1)), s, t) == Heaviside(t)*erf(sqrt(t)) + assert ILT(1/(s**2*(s**2 + 1)), s, t) == (t - sin(t))*Heaviside(t) + assert ILT(s**2/(s**2 + 1), s, t) == -sin(t)*Heaviside(t) + DiracDelta(t) + assert ILT(1 - 1/(s**2 + 1), s, t) == -sin(t)*Heaviside(t) + DiracDelta(t) assert ILT(1/s**2, s, t) == t*Heaviside(t) assert ILT(1/s**5, s, t) == t**4*Heaviside(t)/24 - assert ILT(exp(-a*s)/s, s, t) == Heaviside(t - a) - assert ILT(exp(-a*s)/(s + b), s, t) == exp(b*(a - t))*Heaviside(-a + t) - assert ILT(a/(s**2 + a**2), s, t) == sin(a*t)*Heaviside(t) - assert ILT(s/(s**2 + a**2), s, t) == cos(a*t)*Heaviside(t) - # TODO is there a way around simp_hyp? assert simp_hyp(ILT(a/(s**2 - a**2), s, t)) == sinh(a*t)*Heaviside(t) assert simp_hyp(ILT(s/(s**2 - a**2), s, t)) == cosh(a*t)*Heaviside(t) - assert ILT(a/((s + b)**2 + a**2), s, t) == exp(-b*t)*sin(a*t)*Heaviside(t) - assert ILT( - (s + b)/((s + b)**2 + a**2), s, t) == exp(-b*t)*cos(a*t)*Heaviside(t) # TODO sinh/cosh shifted come out a mess. also delayed trig is a mess # TODO should this simplify further? assert ILT(exp(-a*s)/s**b, s, t) == \ (t - a)**(b - 1)*Heaviside(t - a)/gamma(b) - assert ILT(exp(-a*s)/sqrt(1 + s**2), s, t) == \ Heaviside(t - a)*besselj(0, a - t) # note: besselj(0, x) is even - # XXX ILT turns these branch factor into trig functions ... assert simplify(ILT(a**b*(s + sqrt(s**2 - a**2))**(-b)/sqrt(s**2 - a**2), s, t).rewrite(exp)) == \ @@ -598,7 +632,6 @@ def simp_hyp(expr): assert ILT( (s * eye(2) - Matrix([[1, 0], [0, 2]])).inv(), s, t) ==\ Matrix([[exp(t)*Heaviside(t), 0], [0, exp(2*t)*Heaviside(t)]]) - def test_inverse_laplace_transform_delta(): from sympy import DiracDelta ILT = inverse_laplace_transform diff --git a/sympy/integrals/transforms.py b/sympy/integrals/transforms.py index 0f98d89db120..35b16bcf2d12 100644 --- a/sympy/integrals/transforms.py +++ b/sympy/integrals/transforms.py @@ -7,6 +7,7 @@ from sympy.core.relational import _canonical, Ge, Gt from sympy.core.numbers import oo from sympy.core.symbol import Dummy +from sympy.functions import DiracDelta from sympy.functions.elementary.miscellaneous import Max from sympy.integrals import integrate, Integral from sympy.integrals.meijerint import _dummy @@ -15,6 +16,7 @@ from sympy.utilities import default_sort_key from sympy.utilities.exceptions import SymPyDeprecationWarning from sympy.matrices.matrices import MatrixBase +from sympy.polys.matrices.linsolve import _lin_eq2dict, PolyNonlinearError ########################################################################## @@ -1005,14 +1007,39 @@ def repl(ex, *args): expr = repl(expr, Unequality, replue) return S(expr) +def expand_dirac_delta(expr): + """ + Expand an expression involving DiractDelta to get it as a linear + combination of DiracDelta functions. + """ + return _lin_eq2dict(expr, expr.atoms(DiracDelta)) @_noconds def _laplace_transform(f, t, s_, simplify=True): """ The backend function for Laplace transforms. """ from sympy import (re, Max, exp, pi, Min, periodic_argument as arg_, - arg, cos, Wild, symbols, polar_lift) + arg, cos, Wild, symbols, polar_lift, Add) s = Dummy('s') - F = integrate(exp(-s*t) * f, (t, 0, oo)) + a = Wild('a', exclude=[t]) + deltazero = [] + deltanonzero = [] + try: + integratable, deltadict = expand_dirac_delta(f) + except PolyNonlinearError: + raise IntegralTransformError( + 'Laplace', f, 'could not expand DiracDelta expressions') + for dirac_func, dirac_coeff in deltadict.items(): + p = dirac_func.match(DiracDelta(a*t)) + if p: + deltazero.append(dirac_coeff.subs(t,0)/p[a]) + else: + if dirac_func.args[0].subs(t,0).is_zero: + raise IntegralTransformError('Laplace', f,\ + 'not implemented yet.') + else: + deltanonzero.append(dirac_func*dirac_coeff) + F = Add(integrate(exp(-s*t) * Add(integratable, *deltanonzero), (t, 0, oo)), + Add(*deltazero)) if not F.has(Integral): return _simplify(F.subs(s, s_), simplify), -oo, S.true @@ -1150,18 +1177,25 @@ def laplace_transform(f, t, s, legacy_matrix=True, **hints): r""" Compute the Laplace Transform `F(s)` of `f(t)`, - .. math :: F(s) = \int_0^\infty e^{-st} f(t) \mathrm{d}t. + .. math :: F(s) = \int_{0^{-}}^\infty e^{-st} f(t) \mathrm{d}t. Explanation =========== - For all "sensible" functions, this converges absolutely in a + For all sensible functions, this converges absolutely in a half plane `a < \operatorname{Re}(s)`. - This function returns ``(F, a, cond)`` - where ``F`` is the Laplace transform of ``f``, `\operatorname{Re}(s) > a` is the half-plane + This function returns ``(F, a, cond)`` where ``F`` is the Laplace + transform of ``f``, `\operatorname{Re}(s) > a` is the half-plane of convergence, and ``cond`` are auxiliary convergence conditions. + The lower bound is `0^{-}`, meaning that this bound should be approached + from the lower side. This is only necessary if distributions are involved. + At present, it is only done if `f(t)` contains ``DiracDelta``, in which + case the Laplace transform is computed as + + .. math :: F(s) = \lim_{\tau\to 0^{-}} \int_{\tau}^\infty e^{-st} f(t) \mathrm{d}t. + If the integral cannot be computed in closed form, this function returns an unevaluated :class:`LaplaceTransform` object. @@ -1182,8 +1216,11 @@ def laplace_transform(f, t, s, legacy_matrix=True, **hints): >>> from sympy.integrals import laplace_transform >>> from sympy.abc import t, s, a + >>> from sympy.functions import DiracDelta, exp >>> laplace_transform(t**a, t, s) (gamma(a + 1)/(s*s**a), 0, re(a) > -1) + >>> laplace_transform(DiracDelta(t)-a*exp(-a*t),t,s) + (-a/(a + s) + 1, 0, Abs(arg(a)) <= pi/2) See Also ======== @@ -1220,7 +1257,8 @@ def laplace_transform(f, t, s, legacy_matrix=True, **hints): @_noconds_(True) def _inverse_laplace_transform(F, s, t_, plane, simplify=True): """ The backend function for inverse Laplace transforms. """ - from sympy import exp, Heaviside, log, expand_complex, Integral, Piecewise + from sympy import exp, Heaviside, log, expand_complex, Integral,\ + Piecewise, Add from sympy.integrals.meijerint import meijerint_inversion, _get_coeff_exp # There are two strategies we can try: # 1) Use inverse mellin transforms - related by a simple change of variables. @@ -1240,6 +1278,14 @@ def pw_simp(*args): return Heaviside(1/abs(coeff) - t**exponent)*e1 \ + Heaviside(t**exponent - 1/abs(coeff))*e2 + if F.is_rational_function(s): + F = F.apart(s) + + if F.is_Add: + f = Add(*[_inverse_laplace_transform(X, s, t, plane, simplify)\ + for X in F.args]) + return _simplify(f.subs(t, t_), simplify), True + try: f, cond = inverse_mellin_transform(F, s, exp(-t), (None, oo), needeval=True, noconds=False) @@ -1268,18 +1314,20 @@ def pw_simp(*args): def simp_heaviside(arg, H0=S.Half): a = arg.subs(exp(-t), u) if a.has(t): - return Heaviside(arg,H0) + return Heaviside(arg, H0) rel = _solve_inequality(a > 0, u) if rel.lts == u: k = log(rel.gts) - return Heaviside(t + k,H0) + return Heaviside(t + k, H0) else: k = log(rel.lts) - return Heaviside(-(t + k),H0) + return Heaviside(-(t + k), H0) + f = f.replace(Heaviside, simp_heaviside) def simp_exp(arg): return expand_complex(exp(arg)) + f = f.replace(exp, simp_exp) # TODO it would be nice to fix cosh and sinh ... simplify messes these diff --git a/sympy/logic/algorithms/minisat22_wrapper.py b/sympy/logic/algorithms/minisat22_wrapper.py new file mode 100644 index 000000000000..3db63bb30c55 --- /dev/null +++ b/sympy/logic/algorithms/minisat22_wrapper.py @@ -0,0 +1,46 @@ +from sympy.assumptions.cnf import EncodedCNF + +def minisat22_satisfiable(expr, all_models=False, minimal=False): + + if not isinstance(expr, EncodedCNF): + exprs = EncodedCNF() + exprs.add_prop(expr) + expr = exprs + + from pysat.solvers import Minisat22 + + # Return UNSAT when False (encoded as 0) is present in the CNF + if {0} in expr.data: + if all_models: + return (f for f in [False]) + return False + + r = Minisat22(expr.data) + + if minimal: + r.set_phases([-(i+1) for i in range(r.nof_vars())]) + + if not r.solve(): + return False + + if not all_models: + return {expr.symbols[abs(lit) - 1]: lit > 0 for lit in r.get_model()} + + else: + # Make solutions sympy compatible by creating a generator + def _gen(results): + satisfiable = False + while results.solve(): + sol = results.get_model() + yield {expr.symbols[abs(lit) - 1]: lit > 0 for lit in sol} + if minimal: + results.add_clause([-i for i in sol if i>0]) + else: + results.add_clause([-i for i in sol]) + satisfiable = True + if not satisfiable: + yield False + raise StopIteration + + + return _gen(r) diff --git a/sympy/logic/boolalg.py b/sympy/logic/boolalg.py index 5de530b0a476..fde747b538d9 100644 --- a/sympy/logic/boolalg.py +++ b/sympy/logic/boolalg.py @@ -369,7 +369,7 @@ class when they evaluate to false. Notes ====== - See note in :py:class`sympy.logic.boolalg.BooleanTrue` + See the notes section in :py:class:`sympy.logic.boolalg.BooleanTrue` Examples ======== diff --git a/sympy/logic/inference.py b/sympy/logic/inference.py index e0628f3ec3ab..329186bd0849 100644 --- a/sympy/logic/inference.py +++ b/sympy/logic/inference.py @@ -35,7 +35,7 @@ def literal_symbol(literal): raise ValueError("Argument must be a boolean literal.") -def satisfiable(expr, algorithm=None, all_models=False): +def satisfiable(expr, algorithm=None, all_models=False, minimal=False): """ Check satisfiability of a propositional sentence. Returns a model when it succeeds. @@ -88,6 +88,11 @@ def satisfiable(expr, algorithm=None, all_models=False): # is not installed algorithm = "dpll2" + if algorithm=="minisat22": + pysat = import_module('pysat') + if pysat is None: + algorithm = "dpll2" + if algorithm == "dpll": from sympy.logic.algorithms.dpll import dpll_satisfiable return dpll_satisfiable(expr) @@ -97,6 +102,9 @@ def satisfiable(expr, algorithm=None, all_models=False): elif algorithm == "pycosat": from sympy.logic.algorithms.pycosat_wrapper import pycosat_satisfiable return pycosat_satisfiable(expr, all_models) + elif algorithm == "minisat22": + from sympy.logic.algorithms.minisat22_wrapper import minisat22_satisfiable + return minisat22_satisfiable(expr, all_models, minimal) raise NotImplementedError diff --git a/sympy/logic/tests/test_inference.py b/sympy/logic/tests/test_inference.py index e87567dfea8c..cde7c65bb5c8 100644 --- a/sympy/logic/tests/test_inference.py +++ b/sympy/logic/tests/test_inference.py @@ -122,6 +122,50 @@ def test_dpll2_satisfiable(): assert dpll2_satisfiable( Equivalent(A, B) & ~A ) == {A: False, B: False} +def test_minisat22_satisfiable(): + A, B, C = symbols('A,B,C') + minisat22_satisfiable = lambda expr: satisfiable(expr, algorithm="minisat22") + assert minisat22_satisfiable( A & ~A ) is False + assert minisat22_satisfiable( A & ~B ) == {A: True, B: False} + assert minisat22_satisfiable( + A | B ) in ({A: True}, {B: False}, {A: False, B: True}, {A: True, B: True}, {A: True, B: False}) + assert minisat22_satisfiable( + (~A | B) & (~B | A) ) in ({A: True, B: True}, {A: False, B: False}) + assert minisat22_satisfiable( (A | B) & (~B | C) ) in ({A: True, B: False, C: True}, + {A: True, B: True, C: True}, {A: False, B: True, C: True}, {A: True, B: False, C: False}) + assert minisat22_satisfiable( A & B & C ) == {A: True, B: True, C: True} + assert minisat22_satisfiable( (A | B) & (A >> B) ) in ({B: True, A: False}, + {B: True, A: True}) + assert minisat22_satisfiable( Equivalent(A, B) & A ) == {A: True, B: True} + assert minisat22_satisfiable( Equivalent(A, B) & ~A ) == {A: False, B: False} + +def test_minisat22_minimal_satisfiable(): + A, B, C = symbols('A,B,C') + minisat22_satisfiable = lambda expr, minimal=True: satisfiable(expr, algorithm="minisat22", minimal=True) + assert minisat22_satisfiable( A & ~A ) is False + assert minisat22_satisfiable( A & ~B ) == {A: True, B: False} + assert minisat22_satisfiable( + A | B ) in ({A: True}, {B: False}, {A: False, B: True}, {A: True, B: True}, {A: True, B: False}) + assert minisat22_satisfiable( + (~A | B) & (~B | A) ) in ({A: True, B: True}, {A: False, B: False}) + assert minisat22_satisfiable( (A | B) & (~B | C) ) in ({A: True, B: False, C: True}, + {A: True, B: True, C: True}, {A: False, B: True, C: True}, {A: True, B: False, C: False}) + assert minisat22_satisfiable( A & B & C ) == {A: True, B: True, C: True} + assert minisat22_satisfiable( (A | B) & (A >> B) ) in ({B: True, A: False}, + {B: True, A: True}) + assert minisat22_satisfiable( Equivalent(A, B) & A ) == {A: True, B: True} + assert minisat22_satisfiable( Equivalent(A, B) & ~A ) == {A: False, B: False} + g = satisfiable((A | B | C),algorithm="minisat22",minimal=True,all_models=True) + sol = next(g) + first_solution = {key for key, value in sol.items() if value} + sol=next(g) + second_solution = {key for key, value in sol.items() if value} + sol=next(g) + third_solution = {key for key, value in sol.items() if value} + assert not first_solution <= second_solution + assert not second_solution <= third_solution + assert not first_solution <= third_solution + def test_satisfiable(): A, B, C = symbols('A,B,C') assert satisfiable(A & (A >> B) & ~B) is False diff --git a/sympy/parsing/latex/_parse_latex_antlr.py b/sympy/parsing/latex/_parse_latex_antlr.py index 29908e6a7fc2..2e4a45ad20e0 100644 --- a/sympy/parsing/latex/_parse_latex_antlr.py +++ b/sympy/parsing/latex/_parse_latex_antlr.py @@ -119,7 +119,8 @@ def convert_add(add): elif add.SUB(): lh = convert_add(add.additive(0)) rh = convert_add(add.additive(1)) - return sympy.Add(lh, -1 * rh, evaluate=False) + return sympy.Add(lh, sympy.Mul(-1, rh, evaluate=False), + evaluate=False) else: return convert_mp(add.mp()) @@ -163,11 +164,8 @@ def convert_unary(unary): return convert_unary(nested_unary) elif unary.SUB(): numabs = convert_unary(nested_unary) - if numabs == 1: - # Use Integer(-1) instead of Mul(-1, 1) - return -numabs - else: - return sympy.Mul(-1, convert_unary(nested_unary), evaluate=False) + # Use Integer(-n) instead of Mul(-1, n) + return -numabs elif postfix: return convert_postfix_list(postfix) diff --git a/sympy/parsing/tests/test_latex.py b/sympy/parsing/tests/test_latex.py index 984820935e7c..f5c61cde13dd 100644 --- a/sympy/parsing/tests/test_latex.py +++ b/sympy/parsing/tests/test_latex.py @@ -76,8 +76,8 @@ def test_import(): GOOD_PAIRS = [ (r"0", 0), (r"1", 1), - (r"-3.14", _Mul(-1, 3.14)), - (r"(-7.13)(1.5)", _Mul(_Mul(-1, 7.13), 1.5)), + (r"-3.14", -3.14), + (r"(-7.13)(1.5)", _Mul(-7.13, 1.5)), (r"x", x), (r"2x", 2*x), (r"x^2", x**2), @@ -242,12 +242,14 @@ def test_import(): (r"\int x \, dx", Integral(x, x)), (r"\log_2 x", _log(x, 2)), (r"\log_a x", _log(x, a)), + (r"5^0 - 4^0", _Add(_Pow(5, 0), _Mul(-1, _Pow(4, 0)))), ] + def test_parseable(): from sympy.parsing.latex import parse_latex for latex_str, sympy_expr in GOOD_PAIRS: - assert parse_latex(latex_str) == sympy_expr + assert parse_latex(latex_str) == sympy_expr, latex_str # These bad LaTeX strings should raise a LaTeXParsingError when parsed BAD_STRINGS = [ diff --git a/sympy/physics/vector/frame.py b/sympy/physics/vector/frame.py index aa9e21e4b561..e181a7d36701 100644 --- a/sympy/physics/vector/frame.py +++ b/sympy/physics/vector/frame.py @@ -4,6 +4,8 @@ from sympy.physics.vector.vector import Vector, _check_vector from sympy.utilities.misc import translate +from warnings import warn + __all__ = ['CoordinateSym', 'ReferenceFrame'] @@ -554,6 +556,25 @@ def _dcm(self, parent, parent_orient): self._dcm_dict = self._dlist[0] = {} # Reset the _dcm_cache self._dcm_cache = {} + + else: + #Check for loops and raise warning accordingly. + visited = [] + queue = list(frames) + cont = True #Flag to control queue loop. + while queue and cont: + node = queue.pop(0) + if node not in visited: + visited.append(node) + neighbors = node._dcm_dict.keys() + for neighbor in neighbors: + if neighbor == parent: + warn('Loops are defined among the orientation of frames.' + \ + ' This is likely not desired and may cause errors in your calculations.') + cont = False + break + queue.append(neighbor) + # Add the dcm relationship to _dcm_dict self._dcm_dict.update({parent: parent_orient.T}) parent._dcm_dict.update({self: parent_orient}) @@ -579,6 +600,12 @@ def orient_axis(self, parent, axis, angle): angle : sympifiable Angle in radians by which it the frame is to be rotated. + Warns + ====== + + UserWarning + If the orientation creates a kinematic loop. + Examples ======== @@ -657,6 +684,12 @@ def orient_explicit(self, parent, dcm): Direction cosine matrix that specifies the relative rotation between the two reference frames. + Warns + ====== + + UserWarning + If the orientation creates a kinematic loop. + Examples ======== @@ -761,6 +794,12 @@ def orient_body_fixed(self, parent, angles, rotation_order): Tait-Bryan): zxz, xyx, yzy, zyz, xzx, yxy, xyz, yzx, zxy, xzy, zyx, and yxz. + Warns + ====== + + UserWarning + If the orientation creates a kinematic loop. + Examples ======== @@ -773,6 +812,7 @@ def orient_body_fixed(self, parent, angles, rotation_order): >>> B = ReferenceFrame('B') >>> B1 = ReferenceFrame('B1') >>> B2 = ReferenceFrame('B2') + >>> B3 = ReferenceFrame('B3') For example, a classic Euler Angle rotation can be done by: @@ -790,8 +830,8 @@ def orient_body_fixed(self, parent, angles, rotation_order): >>> B1.orient_axis(N, N.x, q1) >>> B2.orient_axis(B1, B1.y, q2) - >>> B.orient_axis(B2, B2.x, q3) - >>> B.dcm(N) + >>> B3.orient_axis(B2, B2.x, q3) + >>> B3.dcm(N) Matrix([ [ cos(q2), sin(q1)*sin(q2), -sin(q2)*cos(q1)], [sin(q2)*sin(q3), -sin(q1)*sin(q3)*cos(q2) + cos(q1)*cos(q3), sin(q1)*cos(q3) + sin(q3)*cos(q1)*cos(q2)], @@ -871,6 +911,12 @@ def orient_space_fixed(self, parent, angles, rotation_order): ``'131'``, or the integer ``131``. There are 12 unique valid rotation orders. + Warns + ====== + + UserWarning + If the orientation creates a kinematic loop. + Examples ======== @@ -883,6 +929,7 @@ def orient_space_fixed(self, parent, angles, rotation_order): >>> B = ReferenceFrame('B') >>> B1 = ReferenceFrame('B1') >>> B2 = ReferenceFrame('B2') + >>> B3 = ReferenceFrame('B3') >>> B.orient_space_fixed(N, (q1, q2, q3), '312') >>> B.dcm(N) @@ -895,8 +942,8 @@ def orient_space_fixed(self, parent, angles, rotation_order): >>> B1.orient_axis(N, N.z, q1) >>> B2.orient_axis(B1, N.x, q2) - >>> B.orient_axis(B2, N.y, q3) - >>> B.dcm(N).simplify() + >>> B3.orient_axis(B2, N.y, q3) + >>> B3.dcm(N).simplify() Matrix([ [ sin(q1)*sin(q2)*sin(q3) + cos(q1)*cos(q3), sin(q1)*cos(q2), sin(q1)*sin(q2)*cos(q3) - sin(q3)*cos(q1)], [-sin(q1)*cos(q3) + sin(q2)*sin(q3)*cos(q1), cos(q1)*cos(q2), sin(q1)*sin(q3) + sin(q2)*cos(q1)*cos(q3)], @@ -992,6 +1039,12 @@ def orient_quaternion(self, parent, numbers): The four quaternion scalar numbers as defined above: ``q0``, ``q1``, ``q2``, ``q3``. + Warns + ====== + + UserWarning + If the orientation creates a kinematic loop. + Examples ======== @@ -1098,6 +1151,12 @@ def orient(self, parent, rot_type, amounts, rot_order=''): ``'123'`` and integer ``123`` are equivalent, for example. Required for ``'Body'`` and ``'Space'``. + Warns + ====== + + UserWarning + If the orientation creates a kinematic loop. + """ _check_frame(parent) diff --git a/sympy/physics/vector/tests/test_frame.py b/sympy/physics/vector/tests/test_frame.py index 61fe2030bbe2..819ba0cb54b2 100644 --- a/sympy/physics/vector/tests/test_frame.py +++ b/sympy/physics/vector/tests/test_frame.py @@ -6,6 +6,7 @@ from sympy.physics.vector.frame import _check_frame from sympy.physics.vector.vector import VectorTypeError from sympy.testing.pytest import raises +import warnings Vector.simp = True @@ -472,6 +473,22 @@ def test_orient_quaternion(): B.orient_quaternion(A, (0,0,0,0)) assert B.dcm(A) == Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) +def test_looped_frame_warning(): + A = ReferenceFrame('A') + B = ReferenceFrame('B') + C = ReferenceFrame('C') + + a, b, c = symbols('a b c') + B.orient_axis(A, A.x, a) + C.orient_axis(B, B.x, b) + + with warnings.catch_warnings(record = True) as w: + warnings.simplefilter("always") + A.orient_axis(C, C.x, c) + assert issubclass(w[-1].category, UserWarning) + assert 'Loops are defined among the orientation of frames. ' + \ + 'This is likely not desired and may cause errors in your calculations.' in str(w[-1].message) + def test_frame_dict(): A = ReferenceFrame('A') B = ReferenceFrame('B') diff --git a/sympy/polys/constructor.py b/sympy/polys/constructor.py index adcfd13dc3cc..3ae1babce9c0 100644 --- a/sympy/polys/constructor.py +++ b/sympy/polys/constructor.py @@ -48,7 +48,7 @@ def _construct_simple(coeffs, opt): float_numbers.append(x) if y.is_Float: float_numbers.append(y) - if is_algebraic(coeff): + elif is_algebraic(coeff): if floats: # there are both algebraics and reals -> EX return False diff --git a/sympy/polys/domains/expressiondomain.py b/sympy/polys/domains/expressiondomain.py index be0fb61f7afc..ee8d7886e6dd 100644 --- a/sympy/polys/domains/expressiondomain.py +++ b/sympy/polys/domains/expressiondomain.py @@ -64,10 +64,14 @@ def _to_ex(f, g): def __add__(f, g): g = f._to_ex(g) - if g is not None: - return f.simplify(f.ex + g.ex) - else: + if g is None: return NotImplemented + elif g == EX.zero: + return f + elif f == EX.zero: + return g + else: + return f.simplify(f.ex + g.ex) def __radd__(f, g): return f.simplify(f.__class__(g).ex + f.ex) @@ -75,10 +79,14 @@ def __radd__(f, g): def __sub__(f, g): g = f._to_ex(g) - if g is not None: - return f.simplify(f.ex - g.ex) - else: + if g is None: return NotImplemented + elif g == EX.zero: + return f + elif f == EX.zero: + return -g + else: + return f.simplify(f.ex - g.ex) def __rsub__(f, g): return f.simplify(f.__class__(g).ex - f.ex) @@ -86,11 +94,16 @@ def __rsub__(f, g): def __mul__(f, g): g = f._to_ex(g) - if g is not None: - return f.simplify(f.ex*g.ex) - else: + if g is None: return NotImplemented + if EX.zero in (f, g): + return EX.zero + elif f.ex.is_Number and g.ex.is_Number: + return f.__class__(f.ex*g.ex) + + return f.simplify(f.ex*g.ex) + def __rmul__(f, g): return f.simplify(f.__class__(g).ex*f.ex) diff --git a/sympy/polys/matrices/ddm.py b/sympy/polys/matrices/ddm.py index fb002f9dc360..d1770b21b832 100644 --- a/sympy/polys/matrices/ddm.py +++ b/sympy/polys/matrices/ddm.py @@ -30,7 +30,7 @@ >>> ddm_idet([[0, 1], [-1, 0]], QQ) 1 >>> A - [[-1, 0], [0, 1]] + [[-1, 0], [0, -1]] Note that ddm_idet modifies the input matrix in-place. It is recommended to use the DDM.det method as a friendlier interface to this instead which takes @@ -284,7 +284,9 @@ def applyfunc(self, func, domain): def rref(a): """Reduced-row echelon form of a and list of pivots""" b = a.copy() - pivots = ddm_irref(b) + K = a.domain + partial_pivot = K.is_RealField or K.is_ComplexField + pivots = ddm_irref(b, _partial_pivot=partial_pivot) return b, pivots def nullspace(a): diff --git a/sympy/polys/matrices/dense.py b/sympy/polys/matrices/dense.py index 50f4024594f2..dfaff42a2969 100644 --- a/sympy/polys/matrices/dense.py +++ b/sympy/polys/matrices/dense.py @@ -85,7 +85,7 @@ def ddm_imatmul(a, b, c): ai[j] = sum(map(mul, bi, cTj), ai[j]) -def ddm_irref(a): +def ddm_irref(a, _partial_pivot=False): """a <-- rref(a)""" # a is (m x n) m = len(a) @@ -97,6 +97,15 @@ def ddm_irref(a): pivots = [] for j in range(n): + # Proper pivoting should be used for all domains for performance + # reasons but it is only strictly needed for RR and CC (and possibly + # other domains like RR(x)). This path is used by DDM.rref() if the + # domain is RR or CC. It uses partial (row) pivoting based on the + # absolute value of the pivot candidates. + if _partial_pivot: + ip = max(range(i, m), key=lambda ip: abs(a[ip][j])) + a[i], a[ip] = a[ip], a[i] + # pivot aij = a[i][j] @@ -140,7 +149,7 @@ def ddm_irref(a): def ddm_idet(a, K): """a <-- echelon(a); return det""" - # Fraction-free Gaussian elimination + # Bareiss algorithm # https://www.math.usm.edu/perry/Research/Thesis_DRL.pdf # a is (m x n) @@ -149,44 +158,27 @@ def ddm_idet(a, K): return K.one n = len(a[0]) - is_field = K.is_Field - # uf keeps track of the effect of row swaps and multiplies + exquo = K.exquo + # uf keeps track of the sign change from row swaps uf = K.one - for j in range(n-1): - # if zero on the diagonal need to swap - if not a[j][j]: - for l in range(j+1, n): - if a[l][j]: - a[j], a[l] = a[l], a[j] + + for k in range(n-1): + if not a[k][k]: + for i in range(k+1, n): + if a[i][k]: + a[k], a[i] = a[i], a[k] uf = -uf break else: - # unable to swap: det = 0 return K.zero - for i in range(j+1, n): - if a[i][j]: - if not is_field: - d = K.gcd(a[j][j], a[i][j]) - b = a[j][j] // d - c = a[i][j] // d - else: - b = a[j][j] - c = a[i][j] - # account for multiplying row i by b - uf = b * uf - for k in range(j+1, n): - a[i][k] = b*a[i][k] - c*a[j][k] - - # triangular det is product of diagonal - prod = K.one - for i in range(n): - prod = prod * a[i][i] - # incorporate swaps and multiplies - if not is_field: - D = prod // uf - else: - D = prod / uf - return D + + akkm1 = a[k-1][k-1] if k else K.one + + for i in range(k+1, n): + for j in range(k+1, n): + a[i][j] = exquo(a[i][j]*a[k][k] - a[i][k]*a[k][j], akkm1) + + return uf * a[-1][-1] def ddm_iinv(ainv, a, K): diff --git a/sympy/polys/matrices/linsolve.py b/sympy/polys/matrices/linsolve.py index 62e6d90aae36..337a9939ae39 100644 --- a/sympy/polys/matrices/linsolve.py +++ b/sympy/polys/matrices/linsolve.py @@ -73,6 +73,12 @@ def _linsolve(eqs, syms): Aaug = sympy_dict_to_dm(eqsdict, rhs, syms) K = Aaug.domain + # sdm_irref has issues with float matrices. This uses the ddm_rref() + # function. When sdm_rref() can handle float matrices reasonably this + # should be removed... + if K.is_RealField or K.is_ComplexField: + Aaug = Aaug.to_ddm().rref()[0].to_sdm() + # Compute reduced-row echelon form (RREF) Arref, pivots, nzcols = sdm_irref(Aaug) diff --git a/sympy/polys/matrices/sdm.py b/sympy/polys/matrices/sdm.py index cfa624185a5d..7c4ad43660cc 100644 --- a/sympy/polys/matrices/sdm.py +++ b/sympy/polys/matrices/sdm.py @@ -904,6 +904,8 @@ def sdm_irref(A): Ajnz = set(Aj) for k in Ajnz - Ainz: Ai[k] = - Aij * Aj[k] + Ai.pop(j) + Ainz.remove(j) for k in Ajnz & Ainz: Aik = Ai[k] - Aij * Aj[k] if Aik: @@ -938,6 +940,8 @@ def sdm_irref(A): for l in Ainz - Aknz: Ak[l] = - Akj * Ai[l] nonzero_columns[l].add(k) + Ak.pop(j) + Aknz.remove(j) for l in Ainz & Aknz: Akl = Ak[l] - Akj * Ai[l] if Akl: diff --git a/sympy/polys/matrices/tests/test_linsolve.py b/sympy/polys/matrices/tests/test_linsolve.py index eda4cdbdf3fa..a77beebdbab7 100644 --- a/sympy/polys/matrices/tests/test_linsolve.py +++ b/sympy/polys/matrices/tests/test_linsolve.py @@ -7,7 +7,7 @@ from sympy.testing.pytest import raises from sympy import S, Eq, I -from sympy.abc import x, y +from sympy.abc import x, y, z from sympy.polys.matrices.linsolve import _linsolve from sympy.polys.solvers import PolyNonlinearError @@ -23,6 +23,83 @@ def test__linsolve(): raises(PolyNonlinearError, lambda: _linsolve([x*(1 + x)], [x])) +def test__linsolve_float(): + + # This should give the exact answer: + eqs = [ + y - x, + y - 0.0216 * x + ] + sol = {x:0.0, y:0.0} + assert _linsolve(eqs, (x, y)) == sol + + # Other cases should be close to eps + + def all_close(sol1, sol2, eps=1e-15): + close = lambda a, b: abs(a - b) < eps + assert sol1.keys() == sol2.keys() + return all(close(sol1[s], sol2[s]) for s in sol1) + + eqs = [ + 0.8*x + 0.8*z + 0.2, + 0.9*x + 0.7*y + 0.2*z + 0.9, + 0.7*x + 0.2*y + 0.2*z + 0.5 + ] + sol_exact = {x:-29/42, y:-11/21, z:37/84} + sol_linsolve = _linsolve(eqs, [x,y,z]) + assert all_close(sol_exact, sol_linsolve) + + eqs = [ + 0.9*x + 0.3*y + 0.4*z + 0.6, + 0.6*x + 0.9*y + 0.1*z + 0.7, + 0.4*x + 0.6*y + 0.9*z + 0.5 + ] + sol_exact = {x:-88/175, y:-46/105, z:-1/25} + sol_linsolve = _linsolve(eqs, [x,y,z]) + assert all_close(sol_exact, sol_linsolve) + + eqs = [ + 0.4*x + 0.3*y + 0.6*z + 0.7, + 0.4*x + 0.3*y + 0.9*z + 0.9, + 0.7*x + 0.9*y, + ] + sol_exact = {x:-9/5, y:7/5, z:-2/3} + sol_linsolve = _linsolve(eqs, [x,y,z]) + assert all_close(sol_exact, sol_linsolve) + + eqs = [ + x*(0.7 + 0.6*I) + y*(0.4 + 0.7*I) + z*(0.9 + 0.1*I) + 0.5, + 0.2*I*x + 0.2*I*y + z*(0.9 + 0.2*I) + 0.1, + x*(0.9 + 0.7*I) + y*(0.9 + 0.7*I) + z*(0.9 + 0.4*I) + 0.4, + ] + sol_exact = { + x:-6157/7995 - 411/5330*I, + y:8519/15990 + 1784/7995*I, + z:-34/533 + 107/1599*I, + } + sol_linsolve = _linsolve(eqs, [x,y,z]) + assert all_close(sol_exact, sol_linsolve) + + # XXX: This system for x and y over RR(z) is problematic. + # + # eqs = [ + # x*(0.2*z + 0.9) + y*(0.5*z + 0.8) + 0.6, + # 0.1*x*z + y*(0.1*z + 0.6) + 0.9, + # ] + # + # linsolve(eqs, [x, y]) + # The solution for x comes out as + # + # -3.9e-5*z**2 - 3.6e-5*z - 8.67361737988404e-20 + # x = ---------------------------------------------- + # 3.0e-6*z**3 - 1.3e-5*z**2 - 5.4e-5*z + # + # The 8e-20 in the numerator should be zero which would allow z to cancel + # from top and bottom. It should be possible to avoid this somehow because + # the inverse of the matrix only has a quadratic factor (the determinant) + # in the denominator. + + def test__linsolve_deprecated(): assert _linsolve([Eq(x**2, x**2+y)], [x, y]) == {x:x, y:S.Zero} assert _linsolve([(x+y)**2-x**2], [x]) == {x:-y/2} diff --git a/sympy/polys/tests/test_constructor.py b/sympy/polys/tests/test_constructor.py index 0f2c747202fe..0feb681d26ba 100644 --- a/sympy/polys/tests/test_constructor.py +++ b/sympy/polys/tests/test_constructor.py @@ -27,6 +27,9 @@ def test_construct_domain(): assert isinstance(result[0], ComplexField) assert result[1] == [CC(3.14), CC(1.0j), CC(0.5)] + assert construct_domain([1.0+I]) == (CC, [CC(1.0, 1.0)]) + assert construct_domain([2.0+3.0*I]) == (CC, [CC(2.0, 3.0)]) + assert construct_domain([1, I]) == (ZZ_I, [ZZ_I(1, 0), ZZ_I(0, 1)]) assert construct_domain([1, I/2]) == (QQ_I, [QQ_I(1, 0), QQ_I(0, S.Half)]) diff --git a/sympy/polys/tests/test_polytools.py b/sympy/polys/tests/test_polytools.py index c6b7a5c126c5..762e88f720cf 100644 --- a/sympy/polys/tests/test_polytools.py +++ b/sympy/polys/tests/test_polytools.py @@ -51,15 +51,18 @@ from sympy.polys.fields import field from sympy.polys.domains import FF, ZZ, QQ, ZZ_I, QQ_I, RR, EX from sympy.polys.domains.realfield import RealField +from sympy.polys.domains.complexfield import ComplexField from sympy.polys.orderings import lex, grlex, grevlex from sympy import ( S, Integer, Rational, Float, Mul, Symbol, sqrt, Piecewise, Derivative, exp, sin, tanh, expand, oo, I, pi, re, im, rootof, Eq, Tuple, Expr, diff) +from sympy.core.add import Add from sympy.core.basic import _aresame from sympy.core.compatibility import iterable from sympy.core.mul import _keep_coeff +from sympy.core.power import Pow from sympy.testing.pytest import raises, warns_deprecated_sympy from sympy.abc import a, b, c, d, p, q, t, w, x, y, z @@ -387,6 +390,7 @@ def test_Poly__new__(): modulus=65537, symmetric=False) assert isinstance(Poly(x**2 + x + 1.0).get_domain(), RealField) + assert isinstance(Poly(x**2 + x + I + 1.0).get_domain(), ComplexField) def test_Poly__args(): @@ -3324,6 +3328,14 @@ def test_keep_coeff(): assert _keep_coeff(S(2), x + 1) == u assert _keep_coeff(x, 1/x) == 1 assert _keep_coeff(x + 1, S(2)) == u + assert _keep_coeff(S.Half, S.One) == S.Half + p = Pow(2, 3, evaluate=False) + assert _keep_coeff(S(-1), p) == Mul(-1, p, evaluate=False) + a = Add(2, p, evaluate=False) + assert _keep_coeff(S.Half, a, clear=True + ) == Mul(S.Half, a, evaluate=False) + assert _keep_coeff(S.Half, a, clear=False + ) == Add(1, Mul(S.Half, p, evaluate=False), evaluate=False) def test_poly_matching_consistency(): diff --git a/sympy/printing/latex.py b/sympy/printing/latex.py index 13ed00660904..eb663eec7fa6 100644 --- a/sympy/printing/latex.py +++ b/sympy/printing/latex.py @@ -1601,9 +1601,6 @@ def _print_Relational(self, expr): return "%s %s %s" % (self._print(expr.lhs), charmap[expr.rel_op], self._print(expr.rhs)) - def _print_SymbolicRelation(self, expr): - return self._print_Relational(expr) - def _print_Piecewise(self, expr): ecpairs = [r"%s & \text{for}\: %s" % (self._print(e), self._print(c)) for e, c in expr.args[:-1]] diff --git a/sympy/printing/precedence.py b/sympy/printing/precedence.py index 35e96e93815c..a54f2ffce279 100644 --- a/sympy/printing/precedence.py +++ b/sympy/printing/precedence.py @@ -48,7 +48,6 @@ "KroneckerProduct": PRECEDENCE["Mul"], "Equality": PRECEDENCE["Mul"], "Unequality": PRECEDENCE["Mul"], - "Equation": PRECEDENCE["Mul"], } # Sometimes it's not enough to assign a fixed precedence value to a diff --git a/sympy/printing/pretty/pretty.py b/sympy/printing/pretty/pretty.py index e634b90d00be..06bd23478734 100644 --- a/sympy/printing/pretty/pretty.py +++ b/sympy/printing/pretty/pretty.py @@ -220,14 +220,12 @@ def _print_binomial(self, e): def _print_Relational(self, e): op = prettyForm(' ' + xsym(e.rel_op) + ' ') + l = self._print(e.lhs) r = self._print(e.rhs) pform = prettyForm(*stringPict.next(l, op, r)) return pform - def _print_SymbolicRelation(self, expr): - return self._print_Relational(expr) - def _print_Not(self, e): from sympy import Equivalent, Implies if self._use_unicode: diff --git a/sympy/printing/pretty/tests/test_pretty.py b/sympy/printing/pretty/tests/test_pretty.py index b04bde42ae2e..42fae775d346 100644 --- a/sympy/printing/pretty/tests/test_pretty.py +++ b/sympy/printing/pretty/tests/test_pretty.py @@ -9,7 +9,7 @@ SeqPer, SeqFormula, SeqAdd, SeqMul, fourier_series, fps, ITE, Complement, Interval, Intersection, Union, EulerGamma, GoldenRatio, LambertW, airyai, airybi, airyaiprime, airybiprime, fresnelc, fresnels, - Heaviside, dirichlet_eta, diag, MatrixSlice, Eqn,) + Heaviside, dirichlet_eta, diag, MatrixSlice) from sympy.codegen.ast import (Assignment, AddAugmentedAssignment, SubAugmentedAssignment, MulAugmentedAssignment, DivAugmentedAssignment, ModAugmentedAssignment) @@ -1164,11 +1164,6 @@ def test_pretty_relational(): assert pretty(expr) in [ascii_str_1, ascii_str_2] assert upretty(expr) in [ucode_str_1, ucode_str_2] - -def test_SymbolicRelation(): - assert pretty(Eqn(x, y)) == "x = y" - - def test_Assignment(): expr = Assignment(x, y) ascii_str = \ diff --git a/sympy/printing/str.py b/sympy/printing/str.py index e163b4851539..c3fdcdd4358f 100644 --- a/sympy/printing/str.py +++ b/sympy/printing/str.py @@ -762,19 +762,6 @@ def _print_Relational(self, expr): self._relationals.get(expr.rel_op) or expr.rel_op, self.parenthesize(expr.rhs, precedence(expr))) - def _print_SymbolicRelation(self, expr): - charmap = { - "==": "Eqn", - } - - if expr.rel_op in charmap: - return '%s(%s, %s)' % (charmap[expr.rel_op], self._print(expr.lhs), - self._print(expr.rhs)) - - return '%s %s %s' % (self.parenthesize(expr.lhs, precedence(expr)), - self._relationals.get(expr.rel_op) or expr.rel_op, - self.parenthesize(expr.rhs, precedence(expr))) - def _print_ComplexRootOf(self, expr): return "CRootOf(%s, %d)" % (self._print_Add(expr.expr, order='lex'), expr.index) diff --git a/sympy/printing/tests/test_latex.py b/sympy/printing/tests/test_latex.py index dfb9c134550e..f04d04c404a9 100644 --- a/sympy/printing/tests/test_latex.py +++ b/sympy/printing/tests/test_latex.py @@ -21,8 +21,7 @@ SeqAdd, SeqMul, fourier_series, pi, ConditionSet, ComplexRegion, fps, AccumBounds, reduced_totient, primenu, primeomega, SingularityFunction, stieltjes, mathieuc, mathieus, mathieucprime, mathieusprime, - UnevaluatedExpr, Quaternion, I, KroneckerProduct, LambertW, - Ne, Gt, Ge, Lt, Le, Eqn,) + UnevaluatedExpr, Quaternion, I, KroneckerProduct, LambertW) from sympy.ntheory.factor_ import udivisor_sigma @@ -2610,19 +2609,6 @@ def test_issue_17092(): assert latex(Derivative(x_star, x_star,2)) == r'\frac{d^{2}}{d \left(x^{*}\right)^{2}} x^{*}' -def test_Relational(): - assert latex(Eq(x, y)) == "x = y" - assert latex(Ne(x, y)) == r"x \neq y" - assert latex(Lt(x, y)) == "x < y" - assert latex(Gt(x, y)) == "x > y" - assert latex(Le(x, y)) == r"x \leq y" - assert latex(Ge(x, y)) == r"x \geq y" - - -def test_SymbolicRelation(): - assert latex(Eqn(x, y)) == "x = y" - - def test_latex_decimal_separator(): x, y, z, t = symbols('x y z t') diff --git a/sympy/printing/tests/test_str.py b/sympy/printing/tests/test_str.py index 94b0a13e9f3c..690b1a8bbfb6 100644 --- a/sympy/printing/tests/test_str.py +++ b/sympy/printing/tests/test_str.py @@ -5,7 +5,7 @@ symbols, Wild, WildFunction, zeta, zoo, Dummy, Dict, Tuple, FiniteSet, factor, subfactorial, true, false, Equivalent, Xor, Complement, SymmetricDifference, AccumBounds, UnevaluatedExpr, Eq, Ne, Quaternion, Subs, MatrixSymbol, MatrixSlice, - Q, Eqn,) + Q) from sympy.core import Expr, Mul from sympy.core.parameters import _exp_is_pow from sympy.external import import_module @@ -587,10 +587,6 @@ def test_Relational(): assert str(Ne(x, 1) & Ne(x, 2)) == "Ne(x, 1) & Ne(x, 2)" -def test_SymbolicRelation(): - assert str(Eqn(x, y)) == 'Eqn(x, y)' - - def test_AppliedBinaryRelation(): assert str(Q.eq(x, y)) == "Q.eq(x, y)" assert str(Q.ne(x, y)) == "Q.ne(x, y)" diff --git a/sympy/series/tests/test_demidovich.py b/sympy/series/tests/test_demidovich.py index 5c9db347f8d1..728a9eb7b510 100644 --- a/sympy/series/tests/test_demidovich.py +++ b/sympy/series/tests/test_demidovich.py @@ -110,7 +110,7 @@ def test_f1b(): assert limit(x*sin(pi/x), x, oo) == pi # 220 assert limit((1 - cos(x))/x**2, x, 0) == S.Half # 221 assert limit(x*sin(1/x), x, oo) == 1 # 227b - assert limit((cos(m*x) - cos(n*x))/x**2, x, 0) == ((n**2 - m**2)/2) # 232 + assert limit((cos(m*x) - cos(n*x))/x**2, x, 0) == -m**2/2 + n**2/2 # 232 assert limit((tan(x) - sin(x))/x**3, x, 0) == S.Half # 233 assert limit((x - sin(2*x))/(x + sin(3*x)), x, 0) == -Rational(1, 4) # 237 assert limit((1 - sqrt(cos(x)))/x**2, x, 0) == Rational(1, 4) # 239 diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py index 3d5e31373ee7..f7ac655dfe74 100644 --- a/sympy/series/tests/test_limits.py +++ b/sympy/series/tests/test_limits.py @@ -112,7 +112,7 @@ def eval(cls, arg): def test_issue_3885(): - assert limit(x*y + x*z, z, 2) == x*(y + 2) + assert limit(x*y + x*z, z, 2) == x*y + 2*x def test_Limit(): @@ -229,7 +229,7 @@ def test_doit(): assert l.doit() is oo -def test_AccumBounds(): +def test_series_AccumBounds(): assert limit(sin(k) - sin(k + 1), k, oo) == AccumBounds(-2, 2) assert limit(cos(k) - cos(k + 1) + 1, k, oo) == AccumBounds(-1, 3) @@ -243,9 +243,8 @@ def test_AccumBounds(): t2 = Mul(S.Half, Add(AccumBounds(-2, 2), sin(1)), 1/(-cos(1) + 1)) assert limit(simplify(Sum(sin(n).rewrite(exp), (n, 0, k)).doit().rewrite(sin)), k, oo) == t2 - assert limit(frac(x)**x, x, oo) == AccumBounds(0, oo) - assert limit(((sin(x) + 1)/2)**x, x, oo) == AccumBounds(0, oo) - # Possible improvement: AccumBounds(0, 1) + assert limit(frac(x)**x, x, oo) == AccumBounds(0, oo) # wolfram gives (0, 1) + assert limit(((sin(x) + 1)/2)**x, x, oo) == AccumBounds(0, oo) # wolfram says 0 @XFAIL @@ -422,7 +421,7 @@ def test_issue_5740(): def test_issue_6366(): n = Symbol('n', integer=True, positive=True) r = (n + 1)*x**(n + 1)/(x**(n + 1) - 1) - x/(x - 1) - assert limit(r, x, 1) == n/2 + assert limit(r, x, 1).cancel() == n/2 def test_factorial(): @@ -441,7 +440,7 @@ def test_factorial(): def test_issue_6560(): e = (5*x**3/4 - x*Rational(3, 4) + (y*(3*x**2/2 - S.Half) + 35*x**4/8 - 15*x**2/4 + Rational(3, 8))/(2*(y + 1))) - assert limit(e, y, oo) == (5*x**3 + 3*x**2 - 3*x - 1)/4 + assert limit(e, y, oo) == 5*x**3/4 + 3*x**2/4 - 3*x/4 - Rational(1, 4) @XFAIL def test_issue_5172(): @@ -515,7 +514,7 @@ def test_branch_cuts(): def test_issue_6364(): a = Symbol('a') e = z/(1 - sqrt(1 + z)*sin(a)**2 - sqrt(1 - z)*cos(a)**2) - assert limit(e, z, 0).simplify() == 2/cos(2*a) + assert limit(e, z, 0) == 2/(2*cos(a)**2 - 1) def test_issue_4099(): @@ -639,7 +638,7 @@ def test_issue_12769(): 2*F0**(b + 1)*K**(2*b)*a*r*s0*(b**2 - 2*b + 1) - \ 2*F0**(b + 1)*K**(b + 1)*a**2*(b - 1))/((b - 1)*(b**2 - 2*b + 1))))*(b*r - b - r + 1) - assert fx.subs(K, F0).cancel().together() == limit(fx, K, F0).together() + assert fx.subs(K, F0).factor(deep=True) == limit(fx, K, F0).factor(deep=True) def test_issue_13332(): @@ -893,7 +892,7 @@ def test_issue_19586(): def test_issue_13715(): n = Symbol('n') p = Symbol('p', zero=True) - assert limit(n + p, n, 0) == p + assert limit(n + p, n, 0) == 0 def test_issue_15055(): @@ -943,3 +942,7 @@ def test_issue_20704(): def test_issue_21038(): assert limit(sin(pi*x)/(3*x - 12), x, 4) == pi/3 + +def test_issue_21550(): + r = (sqrt(5) - 1)/2 + assert limit((x - r)/(x**2 + x - 1), x, r) == (-1 + sqrt(5))/(5 - sqrt(5)) diff --git a/sympy/series/tests/test_nseries.py b/sympy/series/tests/test_nseries.py index ba7d25316485..bfaf699ca034 100644 --- a/sympy/series/tests/test_nseries.py +++ b/sympy/series/tests/test_nseries.py @@ -382,7 +382,7 @@ def test_series3(): def test_bug4(): w = Symbol("w") e = x/(w**4 + x**2*w**4 + 2*x*w**4)*w**4 - assert e.nseries(w, n=2).removeO() in [x/(1 + 2*x + x**2), + assert e.nseries(w, n=2).removeO().expand() in [x/(1 + 2*x + x**2), 1/(1 + x/2 + 1/x/2)/2, 1/x/(1 + 2/x + x**(-2))] @@ -481,8 +481,8 @@ def test_issue_4441(): assert f.series(x, 0, 5) == 1 - a*x + a**2*x**2 - a**3*x**3 + \ a**4*x**4 + O(x**5) f = 1/(1 + (a + b)*x) - assert f.series(x, 0, 3) == 1 + x*(-a - b) + \ - x**2*(a + b)**2 + O(x**3) + assert f.series(x, 0, 3) == 1 + x*(-a - b)\ + + x**2*(a**2 + 2*a*b + b**2) + O(x**3) def test_issue_4329(): diff --git a/sympy/series/tests/test_residues.py b/sympy/series/tests/test_residues.py index 0a599f6a3e2d..842d31869db9 100644 --- a/sympy/series/tests/test_residues.py +++ b/sympy/series/tests/test_residues.py @@ -1,5 +1,5 @@ -from sympy import (residue, Symbol, Function, sin, I, exp, log, pi, - factorial, sqrt, Rational, cot) +from sympy import (residue, Symbol, Function, sin, I, exp, log, pi, S, + factorial, sqrt, Rational, tan, cot, tanh) from sympy.testing.pytest import XFAIL, raises from sympy.abc import x, z, a, s @@ -78,4 +78,14 @@ def test_issue_14037(): def test_issue_21176(): - assert residue(x**2*cot(pi*x)/(x**4 + 1), x, -sqrt(2)/2 - sqrt(2)*I/2) == 0 + f = x**2*cot(pi*x)/(x**4 + 1) + assert residue(f, x, -sqrt(2)/2 - sqrt(2)*I/2).together() == \ + sqrt(2)*(1 - I)/(8*tan(sqrt(2)*pi/2 + sqrt(2)*I*pi/2)) + + +def test_issue_21177(): + r = -sqrt(3)*tanh(sqrt(3)*pi/2)/3 + a = residue(cot(pi*x)/((x - 1)*(x - 2) + 1), x, S(3)/2 - sqrt(3)*I/2) + b = residue(cot(pi*x)/(x**2 - 3*x + 3), x, S(3)/2 - sqrt(3)*I/2) + assert a == r + assert (b - a).cancel() == 0 diff --git a/sympy/series/tests/test_series.py b/sympy/series/tests/test_series.py index 9bcc5781be04..e8cc90754d14 100644 --- a/sympy/series/tests/test_series.py +++ b/sympy/series/tests/test_series.py @@ -219,14 +219,12 @@ def test_issue_12791(): expr = (-beta**2*varphi*sin(theta) + beta**2*cos(theta) + \ beta*varphi*sin(theta) - beta*cos(theta) - beta + 1)/(beta*cos(theta) - 1)**2 - sol = 0.5/(0.5*cos(theta) - 1)**2 - 0.25*cos(theta)/(0.5*cos(theta) - 1)**2 \ - + (beta - 0.5)*(-0.5*varphi*sin(theta)*cos(theta)/((0.5*cos(theta) - 1) \ - **2*(0.5*cos(theta) - 1.0)) - 1/(0.5*cos(theta) - 1)**2 + 0.5*cos(theta) \ - **2/((0.5*cos(theta) - 1)**2*(0.5*cos(theta) - 1.0)) - 1.0*cos(theta) \ - /((0.5*cos(theta) - 1)**2*(0.5*cos(theta) - 1.0))) + 0.25*varphi* \ - sin(theta)/(0.5*cos(theta) - 1)**2 + O((beta - 0.5)**2, (beta, 0.5)) + sol = 0.5/(0.5*cos(theta) - 1.0)**2 - 0.25*cos(theta)/(0.5*cos(theta)\ + - 1.0)**2 + (beta - 0.5)*(-0.25*varphi*sin(2*theta) - 1.5*cos(theta)\ + + 0.25*cos(2*theta) + 1.25)/(0.5*cos(theta) - 1.0)**3\ + + 0.25*varphi*sin(theta)/(0.5*cos(theta) - 1.0)**2 + O((beta - 0.5)**2, (beta, 0.5)) - assert expr.series(beta, 0.5, 2) == sol + assert expr.series(beta, 0.5, 2).trigsimp() == sol def test_issue_14885(): @@ -323,3 +321,19 @@ def test_issue_20551(): expr = (exp(x)/x).series(x, n=None) terms = [ next(expr) for i in range(3) ] assert terms == [1/x, 1, x/2] + + +def test_issue_20697(): + p_0, p_1, p_2, p_3, b_0, b_1, b_2 = symbols('p_0 p_1 p_2 p_3 b_0 b_1 b_2') + Q = (p_0 + (p_1 + (p_2 + p_3/y)/y)/y)/(1 + ((p_3/(b_0*y) + (b_0*p_2\ + - b_1*p_3)/b_0**2)/y + (b_0**2*p_1 - b_0*b_1*p_2 - p_3*(b_0*b_2\ + - b_1**2))/b_0**3)/y) + assert Q.series(y, n=3) == b_2*y**2 + b_1*y + b_0 + O(y**3) + + +def test_issue_21245(): + fi = (1 + sqrt(5))/2 + assert (1/(1 - x - x**2)).series(x, 1/fi, 1) == \ + (36/(-36*sqrt(5) - 80) + 16*sqrt(5)/(-36*sqrt(5) - 80))/(x - 1/(S.Half\ + + sqrt(5)/2)) - 1220*sqrt(5)/(-6100*sqrt(5) - 13640) - 2728\ + /(-6100*sqrt(5) - 13640) + O(x - 2/(1 + sqrt(5)), (x, 2/(1 + sqrt(5)))) diff --git a/sympy/simplify/cse_main.py b/sympy/simplify/cse_main.py index d4708e808900..8412f699ff2e 100644 --- a/sympy/simplify/cse_main.py +++ b/sympy/simplify/cse_main.py @@ -199,6 +199,8 @@ def get_common_arg_candidates(self, argset, min_func_i=0): """ from collections import defaultdict count_map = defaultdict(lambda: 0) + if not argset: + return count_map funcsets = [self.arg_to_funcset[arg] for arg in argset] # As an optimization below, we handle the largest funcset separately from diff --git a/sympy/simplify/tests/test_cse.py b/sympy/simplify/tests/test_cse.py index f076d8f533ae..60affd09159e 100644 --- a/sympy/simplify/tests/test_cse.py +++ b/sympy/simplify/tests/test_cse.py @@ -546,6 +546,12 @@ def test_unevaluated_mul(): eq = Mul(x + y, x + y, evaluate=False) assert cse(eq) == ([(x0, x + y)], [x0**2]) + def test_issue_18991(): A = MatrixSymbol('A', 2, 2) assert signsimp(-A * A - A) == -A * A - A + + +def test_unevaluated_Mul(): + m = [Mul(1, 2, evaluate=False)] + assert cse(m) == ([], m) diff --git a/sympy/solvers/diophantine/diophantine.py b/sympy/solvers/diophantine/diophantine.py index b9a47ede3128..22967b709f64 100644 --- a/sympy/solvers/diophantine/diophantine.py +++ b/sympy/solvers/diophantine/diophantine.py @@ -470,7 +470,7 @@ class BinaryQuadratic(DiophantineEquationType): .. [1] Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0, [online], Available: http://www.alpertron.com.ar/METHODS.HTM .. [2] Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online], - Available: http://www.jpr2718.org/ax2p.pdf + Available: https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf """ @@ -570,7 +570,7 @@ def solve(self, parameters=None, limit=None) -> DiophantineSolutionSet: # (3) Method used when B**2 - 4*A*C is a square, is described in p. 6 of the below paper # by John P. Robertson. - # http://www.jpr2718.org/ax2p.pdf + # https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf elif is_square(discr): if A != 0: @@ -1923,7 +1923,7 @@ def diop_quadratic(eq, param=symbols("t", integer=True)): .. [1] Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0, [online], Available: http://www.alpertron.com.ar/METHODS.HTM .. [2] Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online], - Available: http://www.jpr2718.org/ax2p.pdf + Available: https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf See Also ======== @@ -2006,7 +2006,7 @@ def diop_DN(D, N, t=symbols("t", integer=True)): .. [1] Solving the generalized Pell equation x**2 - D*y**2 = N, John P. Robertson, July 31, 2004, Pages 16 - 17. [online], Available: - http://www.jpr2718.org/pell.pdf + https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf """ if D < 0: if N == 0: @@ -2349,7 +2349,7 @@ def PQa(P_0, Q_0, D): ========== .. [1] Solving the generalized Pell equation x^2 - Dy^2 = N, John P. - Robertson, July 31, 2004, Pages 4 - 8. http://www.jpr2718.org/pell.pdf + Robertson, July 31, 2004, Pages 4 - 8. https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf """ A_i_2 = B_i_1 = 0 A_i_1 = B_i_2 = 1 @@ -2420,7 +2420,7 @@ def diop_bf_DN(D, N, t=symbols("t", integer=True)): ========== .. [1] Solving the generalized Pell equation x**2 - D*y**2 = N, John P. - Robertson, July 31, 2004, Page 15. http://www.jpr2718.org/pell.pdf + Robertson, July 31, 2004, Page 15. https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf """ D = as_int(D) N = as_int(N) @@ -2501,7 +2501,7 @@ def equivalent(u, v, r, s, D, N): ========== .. [1] Solving the generalized Pell equation x**2 - D*y**2 = N, John P. - Robertson, July 31, 2004, Page 12. http://www.jpr2718.org/pell.pdf + Robertson, July 31, 2004, Page 12. https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf """ return divisible(u*r - D*v*s, N) and divisible(u*s - v*r, N) @@ -2627,7 +2627,7 @@ def transformation_to_DN(eq): .. [1] Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0, John P.Robertson, May 8, 2003, Page 7 - 11. - http://www.jpr2718.org/ax2p.pdf + https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf """ var, coeff, diop_type = classify_diop(eq, _dict=False) @@ -2723,7 +2723,7 @@ def find_DN(eq): .. [1] Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0, John P.Robertson, May 8, 2003, Page 7 - 11. - http://www.jpr2718.org/ax2p.pdf + https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf """ var, coeff, diop_type = classify_diop(eq, _dict=False) if diop_type == BinaryQuadratic.name: diff --git a/sympy/solvers/diophantine/tests/test_diophantine.py b/sympy/solvers/diophantine/tests/test_diophantine.py index 1f97568928ab..26a9280710c2 100644 --- a/sympy/solvers/diophantine/tests/test_diophantine.py +++ b/sympy/solvers/diophantine/tests/test_diophantine.py @@ -195,7 +195,7 @@ def test_quadratic_non_perfect_slow(): def test_DN(): # Most of the test cases were adapted from, # Solving the generalized Pell equation x**2 - D*y**2 = N, John P. Robertson, July 31, 2004. - # http://www.jpr2718.org/pell.pdf + # https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf # others are verified using Wolfram Alpha. # Covers cases where D <= 0 or D > 0 and D is a square or N = 0 diff --git a/sympy/solvers/ode/ode.py b/sympy/solvers/ode/ode.py index 86911dc9e3a0..a66baaf9c367 100644 --- a/sympy/solvers/ode/ode.py +++ b/sympy/solvers/ode/ode.py @@ -991,7 +991,6 @@ class in it. Note that a hint may do this anyway if a = Wild('a', exclude=[f(x)]) d = Wild('d', exclude=[df, f(x).diff(x, 2)]) e = Wild('e', exclude=[df]) - k = Wild('k', exclude=[df]) n = Wild('n', exclude=[x, f(x), df]) c1 = Wild('c1', exclude=[x]) a3 = Wild('a3', exclude=[f(x), df, f(x).diff(x, 2)]) @@ -1055,6 +1054,9 @@ class in it. Note that a hint may do this anyway if Factorable: ('factorable',), RiccatiSpecial: ('Riccati_special_minus2',), SecondNonlinearAutonomousConserved: ('2nd_nonlinear_autonomous_conserved',), + Liouville: ('Liouville',), + Separable: ('separable',), + SeparableReduced: ('separable_reduced',), } for solvercls in solvers: solver = solvercls(ode) @@ -1128,13 +1130,6 @@ class in it. Note that a hint may do this anyway if ## Separable Case: y' == P(y)*Q(x) r[d] = separatevars(r[d]) r[e] = separatevars(r[e]) - # m1[coeff]*m1[x]*m1[y] + m2[coeff]*m2[x]*m2[y]*y' - m1 = separatevars(r[d], dict=True, symbols=(x, y)) - m2 = separatevars(r[e], dict=True, symbols=(x, y)) - if m1 and m2: - r1 = {'m1': m1, 'm2': m2, 'y': y} - matching_hints["separable"] = r1 - matching_hints["separable_Integral"] = r1 ## First order equation with homogeneous coefficients: # dy/dx == F(y/x) or dy/dx == F(x/y) @@ -1187,87 +1182,7 @@ class in it. Note that a hint may do this anyway if matching_hints["linear_coefficients"] = r2 matching_hints["linear_coefficients_Integral"] = r2 - ## Equation of the form y' + (y/x)*H(x^n*y) = 0 - # that can be reduced to separable form - - factor = simplify(x/f(x)*num/den) - - # Try representing factor in terms of x^n*y - # where n is lowest power of x in factor; - # first remove terms like sqrt(2)*3 from factor.atoms(Mul) - num, dem = factor.as_numer_denom() - num = expand(num) - dem = expand(dem) - def _degree(expr, x): - # Made this function to calculate the degree of - # x in an expression. If expr will be of form - # x**p*y, (wheare p can be variables/rationals) then it - # will return p. - for val in expr: - if val.has(x): - if isinstance(val, Pow) and val.as_base_exp()[0] == x: - return (val.as_base_exp()[1]) - elif val == x: - return (val.as_base_exp()[1]) - else: - return _degree(val.args, x) - return 0 - - def _powers(expr): - # this function will return all the different relative power of x w.r.t f(x). - # expr = x**p * f(x)**q then it will return {p/q}. - pows = set() - if isinstance(expr, Add): - exprs = expr.atoms(Add) - elif isinstance(expr, Mul): - exprs = expr.atoms(Mul) - elif isinstance(expr, Pow): - exprs = expr.atoms(Pow) - else: - exprs = {expr} - - for arg in exprs: - if arg.has(x): - _, u = arg.as_independent(x, f(x)) - pow = _degree((u.subs(f(x), y), ), x)/_degree((u.subs(f(x), y), ), y) - pows.add(pow) - return pows - - pows = _powers(num) - pows.update(_powers(dem)) - pows = list(pows) - if(len(pows)==1) and pows[0]!=zoo: - t = Dummy('t') - r2 = {'t': t} - num = num.subs(x**pows[0]*f(x), t) - dem = dem.subs(x**pows[0]*f(x), t) - test = num/dem - free = test.free_symbols - if len(free) == 1 and free.pop() == t: - r2.update({'power' : pows[0], 'u' : test}) - matching_hints['separable_reduced'] = r2 - matching_hints["separable_reduced_Integral"] = r2 - - elif order == 2: - # Liouville ODE in the form - # f(x).diff(x, 2) + g(f(x))*(f(x).diff(x))**2 + h(x)*f(x).diff(x) - # See Goldstein and Braun, "Advanced Methods for the Solution of - # Differential Equations", pg. 98 - - s = d*f(x).diff(x, 2) + e*df**2 + k*df - r = reduced_eq.collect(f(x).diff(x)).match(s) - if r and r[d] != 0: - y = Dummy('y') - g = simplify(r[e]/r[d]).subs(f(x), y) - h = simplify(r[k]/r[d]).subs(f(x), y) - if y in h.free_symbols or x in g.free_symbols: - pass - else: - r = {'g': g, 'h': h, 'y': y} - matching_hints["Liouville"] = r - matching_hints["Liouville_Integral"] = r - # Homogeneous second order differential equation of the form # a3*f(x).diff(x, 2) + b3*f(x).diff(x) + c3 # It has a definite power series solution at point x0 if, b3/a3 and c3/a3 @@ -3161,78 +3076,6 @@ def homogeneous_order(eq, *symbols): return e -def ode_Liouville(eq, func, order, match): - r""" - Solves 2nd order Liouville differential equations. - - The general form of a Liouville ODE is - - .. math:: \frac{d^2 y}{dx^2} + g(y) \left(\! - \frac{dy}{dx}\!\right)^2 + h(x) - \frac{dy}{dx}\text{.} - - The general solution is: - - >>> from sympy import Function, dsolve, Eq, pprint, diff - >>> from sympy.abc import x - >>> f, g, h = map(Function, ['f', 'g', 'h']) - >>> genform = Eq(diff(f(x),x,x) + g(f(x))*diff(f(x),x)**2 + - ... h(x)*diff(f(x),x), 0) - >>> pprint(genform) - 2 2 - /d \ d d - g(f(x))*|--(f(x))| + h(x)*--(f(x)) + ---(f(x)) = 0 - \dx / dx 2 - dx - >>> pprint(dsolve(genform, f(x), hint='Liouville_Integral')) - f(x) - / / - | | - | / | / - | | | | - | - | h(x) dx | | g(y) dy - | | | | - | / | / - C1 + C2* | e dx + | e dy = 0 - | | - / / - - Examples - ======== - - >>> from sympy import Function, dsolve, Eq, pprint - >>> from sympy.abc import x - >>> f = Function('f') - >>> pprint(dsolve(diff(f(x), x, x) + diff(f(x), x)**2/f(x) + - ... diff(f(x), x)/x, f(x), hint='Liouville')) - ________________ ________________ - [f(x) = -\/ C1 + C2*log(x) , f(x) = \/ C1 + C2*log(x) ] - - References - ========== - - - Goldstein and Braun, "Advanced Methods for the Solution of Differential - Equations", pp. 98 - - http://www.maplesoft.com/support/help/Maple/view.aspx?path=odeadvisor/Liouville - - # indirect doctest - - """ - # Liouville ODE: - # f(x).diff(x, 2) + g(f(x))*(f(x).diff(x, 2))**2 + h(x)*f(x).diff(x) - # See Goldstein and Braun, "Advanced Methods for the Solution of - # Differential Equations", pg. 98, as well as - # http://www.maplesoft.com/support/help/view.aspx?path=odeadvisor/Liouville - x = func.args[0] - f = func.func - r = match # f(x).diff(x, 2) + g*f(x).diff(x)**2 + h*f(x).diff(x) - y = r['y'] - C1, C2 = get_numbered_constants(eq, num=2) - int = Integral(exp(Integral(r['g'], y)), (y, None, f(x))) - sol = Eq(int + C1*Integral(exp(-Integral(r['h'], x)), x) + C2, 0) - return sol - - def ode_2nd_power_series_ordinary(eq, func, order, match): r""" Gives a power series solution to a second order homogeneous differential @@ -4280,82 +4123,6 @@ def ode_linear_coefficients(eq, func, order, match): return ode_1st_homogeneous_coeff_best(eq, func, order, match) -def ode_separable_reduced(eq, func, order, match): - r""" - Solves a differential equation that can be reduced to the separable form. - - The general form of this equation is - - .. math:: y' + (y/x) H(x^n y) = 0\text{}. - - This can be solved by substituting `u(y) = x^n y`. The equation then - reduces to the separable form `\frac{u'}{u (\mathrm{power} - H(u))} - - \frac{1}{x} = 0`. - - The general solution is: - - >>> from sympy import Function, dsolve, pprint - >>> from sympy.abc import x, n - >>> f, g = map(Function, ['f', 'g']) - >>> genform = f(x).diff(x) + (f(x)/x)*g(x**n*f(x)) - >>> pprint(genform) - / n \ - d f(x)*g\x *f(x)/ - --(f(x)) + --------------- - dx x - >>> pprint(dsolve(genform, hint='separable_reduced')) - n - x *f(x) - / - | - | 1 - | ------------ dy = C1 + log(x) - | y*(n - g(y)) - | - / - - See Also - ======== - :meth:`sympy.solvers.ode.ode.ode_separable` - - Examples - ======== - - >>> from sympy import Function, pprint - >>> from sympy.solvers.ode.ode import dsolve - >>> from sympy.abc import x - >>> f = Function('f') - >>> d = f(x).diff(x) - >>> eq = (x - x**2*f(x))*d - f(x) - >>> dsolve(eq, hint='separable_reduced') - [Eq(f(x), (1 - sqrt(C1*x**2 + 1))/x), Eq(f(x), (sqrt(C1*x**2 + 1) + 1)/x)] - >>> pprint(dsolve(eq, hint='separable_reduced')) - ___________ ___________ - / 2 / 2 - 1 - \/ C1*x + 1 \/ C1*x + 1 + 1 - [f(x) = ------------------, f(x) = ------------------] - x x - - References - ========== - - - Joel Moses, "Symbolic Integration - The Stormy Decade", Communications - of the ACM, Volume 14, Number 8, August 1971, pp. 558 - """ - - # Arguments are passed in a way so that they are coherent with the - # ode_separable function - x = func.args[0] - f = func.func - y = Dummy('y') - u = match['u'].subs(match['t'], y) - ycoeff = 1/(y*(match['power'] - u)) - m1 = {y: 1, x: -1/x, 'coeff': 1} - m2 = {y: ycoeff, x: 1, 'coeff': 1} - r = {'m1': m1, 'm2': m2, 'y': y, 'hint': x**match['power']*f(x)} - return ode_separable(eq, func, order, r) - - def ode_1st_power_series(eq, func, order, match): r""" The power series solution is a method which gives the Taylor series expansion @@ -5015,71 +4782,6 @@ def _solve_variation_of_parameters(eq, func, order, match): return Eq(f(x), gsol.rhs + psol) -def ode_separable(eq, func, order, match): - r""" - Solves separable 1st order differential equations. - - This is any differential equation that can be written as `P(y) - \tfrac{dy}{dx} = Q(x)`. The solution can then just be found by - rearranging terms and integrating: `\int P(y) \,dy = \int Q(x) \,dx`. - This hint uses :py:meth:`sympy.simplify.simplify.separatevars` as its back - end, so if a separable equation is not caught by this solver, it is most - likely the fault of that function. - :py:meth:`~sympy.simplify.simplify.separatevars` is - smart enough to do most expansion and factoring necessary to convert a - separable equation `F(x, y)` into the proper form `P(x)\cdot{}Q(y)`. The - general solution is:: - - >>> from sympy import Function, dsolve, Eq, pprint - >>> from sympy.abc import x - >>> a, b, c, d, f = map(Function, ['a', 'b', 'c', 'd', 'f']) - >>> genform = Eq(a(x)*b(f(x))*f(x).diff(x), c(x)*d(f(x))) - >>> pprint(genform) - d - a(x)*b(f(x))*--(f(x)) = c(x)*d(f(x)) - dx - >>> pprint(dsolve(genform, f(x), hint='separable_Integral')) - f(x) - / / - | | - | b(y) | c(x) - | ---- dy = C1 + | ---- dx - | d(y) | a(x) - | | - / / - - Examples - ======== - - >>> from sympy import Function, dsolve, Eq - >>> from sympy.abc import x - >>> f = Function('f') - >>> pprint(dsolve(Eq(f(x)*f(x).diff(x) + x, 3*x*f(x)**2), f(x), - ... hint='separable', simplify=False)) - / 2 \ 2 - log\3*f (x) - 1/ x - ---------------- = C1 + -- - 6 2 - - References - ========== - - - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations", - Dover 1963, pp. 52 - - # indirect doctest - - """ - x = func.args[0] - f = func.func - C1 = get_numbered_constants(eq, num=1) - r = match # {'m1':m1, 'm2':m2, 'y':y} - u = r.get('hint', f(x)) # get u from separable_reduced else get f(x) - return Eq(Integral(r['m2']['coeff']*r['m2'][r['y']]/r['m1'][r['y']], - (r['y'], None, u)), Integral(-r['m1']['coeff']*r['m1'][x]/ - r['m2'][x], x) + C1) - - def checkinfsol(eq, infinitesimals, func=None, order=None): r""" This function is used to check if the given infinitesimals are the @@ -6972,4 +6674,4 @@ def _nonlinear_3eq_order1_type5(x, y, z, t, eq): #This import is written at the bottom to avoid circular imports. from .single import (NthAlgebraic, Factorable, FirstLinear, AlmostLinear, Bernoulli, SingleODEProblem, SingleODESolver, RiccatiSpecial, - SecondNonlinearAutonomousConserved, FirstExact) + SecondNonlinearAutonomousConserved, FirstExact, Liouville, Separable, SeparableReduced) diff --git a/sympy/solvers/ode/single.py b/sympy/solvers/ode/single.py index d1f9eb1b818d..1aa13abeb873 100644 --- a/sympy/solvers/ode/single.py +++ b/sympy/solvers/ode/single.py @@ -9,19 +9,18 @@ from typing import Dict, Type from typing import Iterator, List, Optional - -from sympy.core import S +from sympy.core import Add, S, Pow from sympy.core.exprtools import factor_terms from sympy.core.expr import Expr from sympy.core.function import AppliedUndef, Derivative, Function, expand, Subs -from sympy.core.numbers import Float +from sympy.core.numbers import Float, zoo from sympy.core.relational import Equality, Eq from sympy.core.symbol import Symbol, Dummy, Wild from sympy.core.mul import Mul from sympy.functions import exp, sqrt, tan, log from sympy.integrals import Integral from sympy.polys.polytools import cancel, factor -from sympy.simplify import simplify +from sympy.simplify import simplify, separatevars from sympy.simplify.radsimp import fraction from sympy.utilities import numbered_symbols @@ -708,7 +707,7 @@ class Bernoulli(SinglePatternODESolver): Note that the equation is separable when `n = 1` (see the docstring of - :py:meth:`~sympy.solvers.ode.ode.ode_separable`). + :py:meth:`~sympy.solvers.ode.single.Separable`). >>> pprint(dsolve(Eq(f(x).diff(x) + P(x)*f(x), Q(x)*f(x)), f(x), ... hint='separable_Integral')) @@ -980,5 +979,337 @@ def _get_general_solution(self, *, simplify: bool = True): return [Eq(lhs, C2 + x), Eq(lhs, C2 - x)] +class Liouville(SinglePatternODESolver): + r""" + Solves 2nd order Liouville differential equations. + + The general form of a Liouville ODE is + + .. math:: \frac{d^2 y}{dx^2} + g(y) \left(\! + \frac{dy}{dx}\!\right)^2 + h(x) + \frac{dy}{dx}\text{.} + + The general solution is: + + >>> from sympy import Function, dsolve, Eq, pprint, diff + >>> from sympy.abc import x + >>> f, g, h = map(Function, ['f', 'g', 'h']) + >>> genform = Eq(diff(f(x),x,x) + g(f(x))*diff(f(x),x)**2 + + ... h(x)*diff(f(x),x), 0) + >>> pprint(genform) + 2 2 + /d \ d d + g(f(x))*|--(f(x))| + h(x)*--(f(x)) + ---(f(x)) = 0 + \dx / dx 2 + dx + >>> pprint(dsolve(genform, f(x), hint='Liouville_Integral')) + f(x) + / / + | | + | / | / + | | | | + | - | h(x) dx | | g(y) dy + | | | | + | / | / + C1 + C2* | e dx + | e dy = 0 + | | + / / + + Examples + ======== + + >>> from sympy import Function, dsolve, Eq, pprint + >>> from sympy.abc import x + >>> f = Function('f') + >>> pprint(dsolve(diff(f(x), x, x) + diff(f(x), x)**2/f(x) + + ... diff(f(x), x)/x, f(x), hint='Liouville')) + ________________ ________________ + [f(x) = -\/ C1 + C2*log(x) , f(x) = \/ C1 + C2*log(x) ] + + References + ========== + + - Goldstein and Braun, "Advanced Methods for the Solution of Differential + Equations", pp. 98 + - http://www.maplesoft.com/support/help/Maple/view.aspx?path=odeadvisor/Liouville + + # indirect doctest + + """ + hint = "Liouville" + has_integral = True + order = [2] + + def _wilds(self, f, x, order): + d = Wild('d', exclude=[f(x).diff(x), f(x).diff(x, 2)]) + e = Wild('e', exclude=[f(x).diff(x)]) + k = Wild('k', exclude=[f(x).diff(x)]) + return d, e, k + + def _equation(self, fx, x, order): + # Liouville ODE in the form + # f(x).diff(x, 2) + g(f(x))*(f(x).diff(x))**2 + h(x)*f(x).diff(x) + # See Goldstein and Braun, "Advanced Methods for the Solution of + # Differential Equations", pg. 98 + d, e, k = self.wilds() + return d*fx.diff(x, 2) + e*fx.diff(x)**2 + k*fx.diff(x) + + def _verify(self, fx): + d, e, k = self.wilds_match() + self.y = Dummy('y') + x = self.ode_problem.sym + self.g = simplify(e/d).subs(fx, self.y) + self.h = simplify(k/d).subs(fx, self.y) + if self.y in self.h.free_symbols or x in self.g.free_symbols: + return False + return True + + def _get_general_solution(self, *, simplify: bool = True): + d, e, k = self.wilds_match() + fx = self.ode_problem.func + x = self.ode_problem.sym + C1, C2 = self.ode_problem.get_numbered_constants(num=2) + int = Integral(exp(Integral(self.g, self.y)), (self.y, None, fx)) + gen_sol = Eq(int + C1*Integral(exp(-Integral(self.h, x)), x) + C2, 0) + + return [gen_sol] + + +class Separable(SinglePatternODESolver): + r""" + Solves separable 1st order differential equations. + + This is any differential equation that can be written as `P(y) + \tfrac{dy}{dx} = Q(x)`. The solution can then just be found by + rearranging terms and integrating: `\int P(y) \,dy = \int Q(x) \,dx`. + This hint uses :py:meth:`sympy.simplify.simplify.separatevars` as its back + end, so if a separable equation is not caught by this solver, it is most + likely the fault of that function. + :py:meth:`~sympy.simplify.simplify.separatevars` is + smart enough to do most expansion and factoring necessary to convert a + separable equation `F(x, y)` into the proper form `P(x)\cdot{}Q(y)`. The + general solution is:: + + >>> from sympy import Function, dsolve, Eq, pprint + >>> from sympy.abc import x + >>> a, b, c, d, f = map(Function, ['a', 'b', 'c', 'd', 'f']) + >>> genform = Eq(a(x)*b(f(x))*f(x).diff(x), c(x)*d(f(x))) + >>> pprint(genform) + d + a(x)*b(f(x))*--(f(x)) = c(x)*d(f(x)) + dx + >>> pprint(dsolve(genform, f(x), hint='separable_Integral')) + f(x) + / / + | | + | b(y) | c(x) + | ---- dy = C1 + | ---- dx + | d(y) | a(x) + | | + / / + + Examples + ======== + + >>> from sympy import Function, dsolve, Eq + >>> from sympy.abc import x + >>> f = Function('f') + >>> pprint(dsolve(Eq(f(x)*f(x).diff(x) + x, 3*x*f(x)**2), f(x), + ... hint='separable', simplify=False)) + / 2 \ 2 + log\3*f (x) - 1/ x + ---------------- = C1 + -- + 6 2 + + References + ========== + + - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations", + Dover 1963, pp. 52 + + # indirect doctest + + """ + hint = "separable" + has_integral = True + order = [1] + + def _wilds(self, f, x, order): + d = Wild('d', exclude=[f(x).diff(x), f(x).diff(x, 2)]) + e = Wild('e', exclude=[f(x).diff(x)]) + return d, e + + def _equation(self, fx, x, order): + d, e = self.wilds() + return d + e*fx.diff(x) + + def _verify(self, fx): + d, e = self.wilds_match() + self.y = Dummy('y') + x = self.ode_problem.sym + d = separatevars(d.subs(fx, self.y)) + e = separatevars(e.subs(fx, self.y)) + # m1[coeff]*m1[x]*m1[y] + m2[coeff]*m2[x]*m2[y]*y' + self.m1 = separatevars(d, dict=True, symbols=(x, self.y)) + self.m2 = separatevars(e, dict=True, symbols=(x, self.y)) + if self.m1 and self.m2: + return True + return False + + def _get_match_object(self): + fx = self.ode_problem.func + x = self.ode_problem.sym + return self.m1, self.m2, x, fx + + def _get_general_solution(self, *, simplify: bool = True): + m1, m2, x, fx = self._get_match_object() + (C1, ) = self.ode_problem.get_numbered_constants(num=1) + int = Integral(m2['coeff']*m2[self.y]/m1[self.y], + (self.y, None, fx)) + gen_sol = Eq(int, Integral(-m1['coeff']*m1[x]/ + m2[x], x) + C1) + return [gen_sol] + + +class SeparableReduced(Separable): + r""" + Solves a differential equation that can be reduced to the separable form. + + The general form of this equation is + + .. math:: y' + (y/x) H(x^n y) = 0\text{}. + + This can be solved by substituting `u(y) = x^n y`. The equation then + reduces to the separable form `\frac{u'}{u (\mathrm{power} - H(u))} - + \frac{1}{x} = 0`. + + The general solution is: + + >>> from sympy import Function, dsolve, pprint + >>> from sympy.abc import x, n + >>> f, g = map(Function, ['f', 'g']) + >>> genform = f(x).diff(x) + (f(x)/x)*g(x**n*f(x)) + >>> pprint(genform) + / n \ + d f(x)*g\x *f(x)/ + --(f(x)) + --------------- + dx x + >>> pprint(dsolve(genform, hint='separable_reduced')) + n + x *f(x) + / + | + | 1 + | ------------ dy = C1 + log(x) + | y*(n - g(y)) + | + / + + See Also + ======== + :meth:`sympy.solvers.ode.single.Separable` + + Examples + ======== + + >>> from sympy import Function, pprint + >>> from sympy.solvers.ode.ode import dsolve + >>> from sympy.abc import x + >>> f = Function('f') + >>> d = f(x).diff(x) + >>> eq = (x - x**2*f(x))*d - f(x) + >>> dsolve(eq, hint='separable_reduced') + [Eq(f(x), (1 - sqrt(C1*x**2 + 1))/x), Eq(f(x), (sqrt(C1*x**2 + 1) + 1)/x)] + >>> pprint(dsolve(eq, hint='separable_reduced')) + ___________ ___________ + / 2 / 2 + 1 - \/ C1*x + 1 \/ C1*x + 1 + 1 + [f(x) = ------------------, f(x) = ------------------] + x x + + References + ========== + + - Joel Moses, "Symbolic Integration - The Stormy Decade", Communications + of the ACM, Volume 14, Number 8, August 1971, pp. 558 + """ + hint = "separable_reduced" + has_integral = True + order = [1] + + def _degree(self, expr, x): + # Made this function to calculate the degree of + # x in an expression. If expr will be of form + # x**p*y, (wheare p can be variables/rationals) then it + # will return p. + for val in expr: + if val.has(x): + if isinstance(val, Pow) and val.as_base_exp()[0] == x: + return (val.as_base_exp()[1]) + elif val == x: + return (val.as_base_exp()[1]) + else: + return self._degree(val.args, x) + return 0 + + def _powers(self, expr): + # this function will return all the different relative power of x w.r.t f(x). + # expr = x**p * f(x)**q then it will return {p/q}. + pows = set() + fx = self.ode_problem.func + x = self.ode_problem.sym + self.y = Dummy('y') + if isinstance(expr, Add): + exprs = expr.atoms(Add) + elif isinstance(expr, Mul): + exprs = expr.atoms(Mul) + elif isinstance(expr, Pow): + exprs = expr.atoms(Pow) + else: + exprs = {expr} + + for arg in exprs: + if arg.has(x): + _, u = arg.as_independent(x, fx) + pow = self._degree((u.subs(fx, self.y), ), x)/self._degree((u.subs(fx, self.y), ), self.y) + pows.add(pow) + return pows + + def _verify(self, fx): + num, den = self.wilds_match() + x = self.ode_problem.sym + factor = simplify(x/fx*num/den) + # Try representing factor in terms of x^n*y + # where n is lowest power of x in factor; + # first remove terms like sqrt(2)*3 from factor.atoms(Mul) + num, dem = factor.as_numer_denom() + num = expand(num) + dem = expand(dem) + pows = self._powers(num) + pows.update(self._powers(dem)) + pows = list(pows) + if(len(pows)==1) and pows[0]!=zoo: + self.t = Dummy('t') + self.r2 = {'t': self.t} + num = num.subs(x**pows[0]*fx, self.t) + dem = dem.subs(x**pows[0]*fx, self.t) + test = num/dem + free = test.free_symbols + if len(free) == 1 and free.pop() == self.t: + self.r2.update({'power' : pows[0], 'u' : test}) + return True + return False + return False + + def _get_match_object(self): + fx = self.ode_problem.func + x = self.ode_problem.sym + u = self.r2['u'].subs(self.r2['t'], self.y) + ycoeff = 1/(self.y*(self.r2['power'] - u)) + m1 = {self.y: 1, x: -1/x, 'coeff': 1} + m2 = {self.y: ycoeff, x: 1, 'coeff': 1} + return m1, m2, x, x**self.r2['power']*fx + + # Avoid circular import: from .ode import dsolve diff --git a/sympy/solvers/ode/tests/test_single.py b/sympy/solvers/ode/tests/test_single.py index f67827d3a220..4ced7024f41f 100644 --- a/sympy/solvers/ode/tests/test_single.py +++ b/sympy/solvers/ode/tests/test_single.py @@ -916,7 +916,7 @@ def _get_examples_ode_sol_factorable(): 'fact_16': { 'eq': f(x).diff(x)**2 - f(x)**3, - 'sol': [Eq(f(x), 4/(C1**2 - 2*C1*x + x**2))] + 'sol': [Eq(f(x), 4/(C1**2 - 2*C1*x + x**2))], }, # kamke ode 1.1 diff --git a/sympy/tensor/indexed.py b/sympy/tensor/indexed.py index 2ef363af183e..4742d9b5284c 100644 --- a/sympy/tensor/indexed.py +++ b/sympy/tensor/indexed.py @@ -350,6 +350,10 @@ def free_symbols(self): @property def expr_free_symbols(self): + from sympy.utilities.exceptions import SymPyDeprecationWarning + SymPyDeprecationWarning(feature="expr_free_symbols method", + issue=21494, + deprecated_since_version="1.9").warn() return {self} diff --git a/sympy/utilities/lambdify.py b/sympy/utilities/lambdify.py index 52cc3d7e910e..cecfacc87e2f 100644 --- a/sympy/utilities/lambdify.py +++ b/sympy/utilities/lambdify.py @@ -83,7 +83,9 @@ "betainc_regularized": "betainc", } -NUMPY_TRANSLATIONS = {} # type: Dict[str, str] +NUMPY_TRANSLATIONS = { + "Heaviside": "heaviside", + } # type: Dict[str, str] SCIPY_TRANSLATIONS = {} # type: Dict[str, str] CUPY_TRANSLATIONS = {} # type: Dict[str, str] diff --git a/sympy/utilities/tests/test_lambdify.py b/sympy/utilities/tests/test_lambdify.py index b76e3ea0fb21..4752fa779630 100644 --- a/sympy/utilities/tests/test_lambdify.py +++ b/sympy/utilities/tests/test_lambdify.py @@ -1256,11 +1256,24 @@ def test_issue_17898(): f_ = lambdify([x], sympy.LambertW(x,-1), modules='scipy') assert f_(0.1) == mpmath.lambertw(0.1, -1) +def test_issue_13167_21411(): + if not numpy: + skip("numpy not installed") + f1 = lambdify(x, sympy.Heaviside(x)) + f2 = lambdify(x, sympy.Heaviside(x, 1)) + res1 = f1([-1, 0, 1]) + res2 = f2([-1, 0, 1]) + assert Abs(res1[0]).n() < 1e-15 # First functionality: only one argument passed + assert Abs(res1[1] - 1/2).n() < 1e-15 + assert Abs(res1[2] - 1).n() < 1e-15 + assert Abs(res2[0]).n() < 1e-15 # Second functionality: two arguments passed + assert Abs(res2[1] - 1).n() < 1e-15 + assert Abs(res2[2] - 1).n() < 1e-15 + def test_single_e(): f = lambdify(x, E) assert f(23) == exp(1.0) - def test_issue_16536(): if not scipy: skip("scipy not installed") diff --git a/sympy/utilities/tests/test_pickling.py b/sympy/utilities/tests/test_pickling.py index 73732f189a14..2f3bf75d4e01 100644 --- a/sympy/utilities/tests/test_pickling.py +++ b/sympy/utilities/tests/test_pickling.py @@ -34,6 +34,7 @@ '_assumptions', # This is a local cache that isn't automatically filled on creation '_mhash', # Cached after __hash__ is called but set to None after creation 'is_EmptySet', # Deprecated from SymPy 1.5. This can be removed when is_EmptySet is removed. + 'expr_free_symbols', # Deprecated from SymPy 1.9. This can be removed when exr_free_symbols is removed. }