Repository: sympy/sympy

Issue:
Question about sympy.solve / sympy.nonlinsolve
I am trying to solve the nonlinear system of equations below. (These are the three 'law of cosine' equations, I found them in the internet.)
```
a1, b1, c1, ca, cb, cg = sm.symbols('a1, b1, c1, ca, cb, cg')
eq1 = a1*a1 + b1*b1 - 2.*a1*b1*cg - c1*c1
eq2 = a1*a1 + c1*c1 - 2.*a1*c1*cb - b1*b1
eq3 = b1*b1 + c1*c1 - 2.*b1*c1*ca - a1*a1

loesung = sm.nonlinsolve([eq1, eq2, eq3], [c1, cb, cg])
```
I first tried sympy.solve.
It gave a solution, actually two solutions, but in the documentation I could not find, that sympy.solve will solve systems of nonlinear equations.
Also, the solutions seemed to be wrong, e.g. the cos(angle) > 1

Then I tried sympy.nonlinsolve.
After about 45 min of running, I aborted it.

My question:
Is it normal, that it takes so long to find a solution, or is this an indication, that none exists?

Many thanks for any explanations!

Discussion:

Patch:
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index f3839e22d1..aa3411864e 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -4248,6 +4248,53 @@ def equal_valued(x, y):
         return (1 << neg_exp) == q
 
 
+def all_close(expr1, expr2, rtol=1e-5, atol=1e-8):
+    """Return True if expr1 and expr2 are numerically close.
+
+    The expressions must have the same structure, but any Rational, Integer, or
+    Float numbers they contain are compared approximately using rtol and atol.
+    Any other parts of expressions are compared exactly.
+
+    Relative tolerance is measured with respect to expr2 so when used in
+    testing expr2 should be the expected correct answer.
+
+    Examples
+    ========
+
+    >>> from sympy import exp
+    >>> from sympy.abc import x, y
+    >>> from sympy.core.numbers import all_close
+    >>> expr1 = 0.1*exp(x - y)
+    >>> expr2 = exp(x - y)/10
+    >>> expr1
+    0.1*exp(x - y)
+    >>> expr2
+    exp(x - y)/10
+    >>> expr1 == expr2
+    False
+    >>> all_close(expr1, expr2)
+    True
+    """
+    NUM_TYPES = (Rational, Float)
+
+    def _all_close(expr1, expr2, rtol, atol):
+        num1 = isinstance(expr1, NUM_TYPES)
+        num2 = isinstance(expr2, NUM_TYPES)
+        if num1 != num2:
+            return False
+        elif num1:
+            return bool(abs(expr1 - expr2) <= atol + rtol*abs(expr2))
+        elif expr1.is_Atom:
+            return expr1 == expr2
+        elif expr1.func != expr2.func or len(expr1.args) != len(expr2.args):
+            return False
+        else:
+            args = zip(expr1.args, expr2.args)
+            return all(_all_close(a1, a2, rtol, atol) for a1, a2 in args)
+
+    return _all_close(_sympify(expr1), _sympify(expr2), rtol, atol)
+
+
 @dispatch(Tuple, Number) # type:ignore
 def _eval_is_eq(self, other): # noqa: F811
     return False
diff --git a/sympy/polys/domains/complexfield.py b/sympy/polys/domains/complexfield.py
index a36e94ebdf..4642b20249 100644
--- a/sympy/polys/domains/complexfield.py
+++ b/sympy/polys/domains/complexfield.py
@@ -5,6 +5,7 @@
 from sympy.core.numbers import Float, I
 from sympy.polys.domains.characteristiczero import CharacteristicZero
 from sympy.polys.domains.field import Field
+from sympy.polys.domains.gaussiandomains import QQ_I
 from sympy.polys.domains.mpelements import MPContext
 from sympy.polys.domains.simpledomain import SimpleDomain
 from sympy.polys.polyerrors import DomainError, CoercionFailed
@@ -136,7 +137,7 @@ def get_ring(self):
 
     def get_exact(self):
         """Returns an exact domain associated with ``self``. """
-        raise DomainError("there is no exact domain associated with %s" % self)
+        return QQ_I
 
     def is_negative(self, element):
         """Returns ``False`` for any ``ComplexElement``. """
diff --git a/sympy/polys/domains/compositedomain.py b/sympy/polys/domains/compositedomain.py
index 560720a024..a8f63ba7bb 100644
--- a/sympy/polys/domains/compositedomain.py
+++ b/sympy/polys/domains/compositedomain.py
@@ -30,3 +30,23 @@ def drop(self, *symbols):
             return domain
         else:
             return self.__class__(domain, newsyms, self.order)
+
+    def set_domain(self, domain):
+        """Set the ground domain of this domain. """
+        return self.__class__(domain, self.symbols, self.order)
+
+    @property
+    def is_Exact(self):
+        """Returns ``True`` if this domain is exact. """
+        return self.domain.is_Exact
+
+    def get_exact(self):
+        """Returns an exact version of this domain. """
+        return self.set_domain(self.domain.get_exact())
+
+    @property
+    def has_CharacteristicZero(self):
+        return self.domain.has_CharacteristicZero
+
+    def characteristic(self):
+        return self.domain.characteristic()
diff --git a/sympy/polys/domains/fractionfield.py b/sympy/polys/domains/fractionfield.py
index 6dbfc904a4..47bc25436b 100644
--- a/sympy/polys/domains/fractionfield.py
+++ b/sympy/polys/domains/fractionfield.py
@@ -49,13 +49,6 @@ def one(self):
     def order(self):
         return self.field.order
 
-    @property
-    def is_Exact(self):
-        return self.domain.is_Exact
-
-    def get_exact(self):
-        return FractionField(self.domain.get_exact(), self.symbols)
-
     def __str__(self):
         return str(self.domain) + '(' + ','.join(map(str, self.symbols)) + ')'
 
@@ -68,13 +61,6 @@ def __eq__(self, other):
             (self.dtype.field, self.domain, self.symbols) ==\
             (other.dtype.field, other.domain, other.symbols)
 
-    @property
-    def has_CharacteristicZero(self):
-        return self.domain.has_CharacteristicZero
-
-    def characteristic(self):
-        return self.domain.characteristic()
-
     def to_sympy(self, a):
         """Convert ``a`` to a SymPy object. """
         return a.as_expr()
diff --git a/sympy/polys/domains/gaussiandomains.py b/sympy/polys/domains/gaussiandomains.py
index e8b9701d91..bf3df50d5d 100644
--- a/sympy/polys/domains/gaussiandomains.py
+++ b/sympy/polys/domains/gaussiandomains.py
@@ -678,4 +678,9 @@ def from_GaussianRationalField(K1, a, K0):
         """Convert a QQ_I element to QQ_I."""
         return a
 
+    def from_ComplexField(K1, a, K0):
+        """Convert a ComplexField element to QQ_I."""
+        return K1.new(QQ.convert(a.real), QQ.convert(a.imag))
+
+
 QQ_I = GaussianRational._parent = GaussianRationalField()
diff --git a/sympy/polys/domains/old_fractionfield.py b/sympy/polys/domains/old_fractionfield.py
index c05e6f99c1..25d849c39e 100644
--- a/sympy/polys/domains/old_fractionfield.py
+++ b/sympy/polys/domains/old_fractionfield.py
@@ -3,14 +3,13 @@
 
 from sympy.polys.domains.field import Field
 from sympy.polys.domains.compositedomain import CompositeDomain
-from sympy.polys.domains.characteristiczero import CharacteristicZero
 from sympy.polys.polyclasses import DMF
 from sympy.polys.polyerrors import GeneratorsNeeded
 from sympy.polys.polyutils import dict_from_basic, basic_from_dict, _dict_reorder
 from sympy.utilities import public
 
 @public
-class FractionField(Field, CharacteristicZero, CompositeDomain):
+class FractionField(Field, CompositeDomain):
     """A class for representing rational function fields. """
 
     dtype = DMF
@@ -32,6 +31,10 @@ def __init__(self, dom, *gens):
         self.domain = self.dom = dom
         self.symbols = self.gens = gens
 
+    def set_domain(self, dom):
+        """Make a new fraction field with given domain. """
+        return self.__class__(dom, *self.gens)
+
     def new(self, element):
         return self.dtype(element, self.dom, len(self.gens) - 1)
 
@@ -46,13 +49,6 @@ def __eq__(self, other):
         return isinstance(other, FractionField) and \
             self.dtype == other.dtype and self.dom == other.dom and self.gens == other.gens
 
-    @property
-    def has_CharacteristicZero(self):
-        return self.dom.has_CharacteristicZero
-
-    def characteristic(self):
-        return self.dom.characteristic()
-
     def to_sympy(self, a):
         """Convert ``a`` to a SymPy object. """
         return (basic_from_dict(a.numer().to_sympy_dict(), *self.gens) /
diff --git a/sympy/polys/domains/old_polynomialring.py b/sympy/polys/domains/old_polynomialring.py
index bd637b257d..c29a4529aa 100644
--- a/sympy/polys/domains/old_polynomialring.py
+++ b/sympy/polys/domains/old_polynomialring.py
@@ -2,7 +2,6 @@
 
 
 from sympy.polys.agca.modules import FreeModulePolyRing
-from sympy.polys.domains.characteristiczero import CharacteristicZero
 from sympy.polys.domains.compositedomain import CompositeDomain
 from sympy.polys.domains.old_fractionfield import FractionField
 from sympy.polys.domains.ring import Ring
@@ -14,10 +13,9 @@
 from sympy.utilities import public
 from sympy.utilities.iterables import iterable
 
-# XXX why does this derive from CharacteristicZero???
 
 @public
-class PolynomialRingBase(Ring, CharacteristicZero, CompositeDomain):
+class PolynomialRingBase(Ring, CompositeDomain):
     """
     Base class for generalized polynomial rings.
 
@@ -47,6 +45,10 @@ def __init__(self, dom, *gens, **opts):
         # NOTE 'order' may not be set if inject was called through CompositeDomain
         self.order = opts.get('order', monomial_key(self.default_order))
 
+    def set_domain(self, dom):
+        """Return a new polynomial ring with given domain. """
+        return self.__class__(dom, *self.gens, order=self.order)
+
     def new(self, element):
         return self.dtype(element, self.dom, len(self.gens) - 1)
 
@@ -72,13 +74,6 @@ def __eq__(self, other):
             self.dtype == other.dtype and self.dom == other.dom and \
             self.gens == other.gens and self.order == other.order
 
-    @property
-    def has_CharacteristicZero(self):
-        return self.dom.has_CharacteristicZero
-
-    def characteristic(self):
-        return self.dom.characteristic()
-
     def from_ZZ(K1, a, K0):
         """Convert a Python ``int`` object to ``dtype``. """
         return K1._ground_new(K1.dom.convert(a, K0))
diff --git a/sympy/polys/domains/polynomialring.py b/sympy/polys/domains/polynomialring.py
index 0063b6f34b..bad73208f8 100644
--- a/sympy/polys/domains/polynomialring.py
+++ b/sympy/polys/domains/polynomialring.py
@@ -67,13 +67,6 @@ def __eq__(self, other):
             (self.dtype.ring, self.domain, self.symbols) == \
             (other.dtype.ring, other.domain, other.symbols)
 
-    @property
-    def has_CharacteristicZero(self):
-        return self.domain.has_CharacteristicZero
-
-    def characteristic(self):
-        return self.domain.characteristic()
-
     def is_unit(self, a):
         """Returns ``True`` if ``a`` is a unit of ``self``"""
         if not a.is_ground:
diff --git a/sympy/polys/euclidtools.py b/sympy/polys/euclidtools.py
index 1a919e1f10..768a44a949 100644
--- a/sympy/polys/euclidtools.py
+++ b/sympy/polys/euclidtools.py
@@ -1489,7 +1489,19 @@ def dup_inner_gcd(f, g, K):
     (x - 1, x + 1, x - 2)
 
     """
-    if not K.is_Exact:
+    # XXX: This used to check for K.is_Exact but leads to awkward results when
+    # the domain is something like RR[z] e.g.:
+    #
+    # >>> g, p, q = Poly(1, x).cancel(Poly(51.05*x*y - 1.0, x))
+    # >>> g
+    # 1.0
+    # >>> p
+    # Poly(17592186044421.0, x, domain='RR[y]')
+    # >>> q
+    # Poly(898081097567692.0*y*x - 17592186044421.0, x, domain='RR[y]'))
+    #
+    # Maybe it would be better to flatten into multivariate polynomials first.
+    if K.is_RR or K.is_CC:
         try:
             exact = K.get_exact()
         except DomainError:


Test patch:
diff --git a/sympy/core/tests/test_numbers.py b/sympy/core/tests/test_numbers.py
index 82a221cf44..fb4c6a86e2 100644
--- a/sympy/core/tests/test_numbers.py
+++ b/sympy/core/tests/test_numbers.py
@@ -10,7 +10,7 @@
 from sympy.core.numbers import (mpf_norm, seterr,
     Integer, I, pi, comp, Rational, E, nan,
     oo, AlgebraicNumber, Number, Float, zoo, equal_valued,
-    int_valued)
+    int_valued, all_close)
 from sympy.core.intfunc import (igcd, igcdex, igcd2, igcd_lehmer,
     ilcm, integer_nthroot, isqrt, integer_log, mod_inverse)
 from sympy.core.power import Pow
@@ -2280,3 +2280,18 @@ def test_equal_valued():
                     continue
                 for value_j in values_n:
                     assert equal_valued(value_i, value_j) is False
+
+
+def test_all_close():
+    x = Symbol('x')
+    assert all_close(2, 2) is True
+    assert all_close(2, 2.0000) is True
+    assert all_close(2, 2.0001) is False
+    assert all_close(1/3, 1/3.0001) is False
+    assert all_close(1/3, 1/3.0001, 1e-3, 1e-3) is True
+    assert all_close(1/3, Rational(1, 3)) is True
+    assert all_close(0.1*exp(0.2*x), exp(x/5)/10) is True
+    # The expressions should be structurally the same:
+    assert all_close(1.4142135623730951, sqrt(2)) is False
+    assert all_close(1.4142135623730951, sqrt(2).evalf()) is True
+    assert all_close(x + 1e-20, x) is False
diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py
index 994c7e5b2c..2b739d184f 100644
--- a/sympy/integrals/tests/test_integrals.py
+++ b/sympy/integrals/tests/test_integrals.py
@@ -5,7 +5,7 @@
 from sympy.core.expr import Expr
 from sympy.core.function import (Derivative, Function, Lambda, diff)
 from sympy.core import EulerGamma
-from sympy.core.numbers import (E, Float, I, Rational, nan, oo, pi, zoo)
+from sympy.core.numbers import (E, I, Rational, nan, oo, pi, zoo, all_close)
 from sympy.core.relational import (Eq, Ne)
 from sympy.core.singleton import S
 from sympy.core.symbol import (Symbol, symbols)
@@ -438,13 +438,12 @@ def test_issue_18133():
 
 
 def test_issue_21741():
-    a = Float('3999999.9999999995', precision=53)
-    b = Float('2.5000000000000004e-7', precision=53)
-    r = Piecewise((b*I*exp(-a*I*pi*t*y)*exp(-a*I*pi*x*z)/(pi*x),
-                   Ne(1.0*pi*x*exp(a*I*pi*t*y), 0)),
+    a = 4e6
+    b = 2.5e-7
+    r = Piecewise((b*I*exp(-a*I*pi*t*y)*exp(-a*I*pi*x*z)/(pi*x), Ne(x, 0)),
                   (z*exp(-a*I*pi*t*y), True))
     fun = E**((-2*I*pi*(z*x+t*y))/(500*10**(-9)))
-    assert integrate(fun, z) == r
+    assert all_close(integrate(fun, z), r)
 
 
 def test_matrices():
diff --git a/sympy/polys/domains/tests/test_domains.py b/sympy/polys/domains/tests/test_domains.py
index 4eb97d12c6..3fe40cdb8f 100644
--- a/sympy/polys/domains/tests/test_domains.py
+++ b/sympy/polys/domains/tests/test_domains.py
@@ -562,22 +562,70 @@ def test_Domain_get_field():
     assert QQ[x, y].get_field() == QQ.frac_field(x, y)
 
 
+def test_Domain_set_domain():
+    doms = [GF(5), ZZ, QQ, ALG, RR, CC, EX, ZZ[z], QQ[z], RR[z], CC[z], EX[z]]
+    for D1 in doms:
+        for D2 in doms:
+            assert D1[x].set_domain(D2) == D2[x]
+            assert D1[x, y].set_domain(D2) == D2[x, y]
+            assert D1.frac_field(x).set_domain(D2) == D2.frac_field(x)
+            assert D1.frac_field(x, y).set_domain(D2) == D2.frac_field(x, y)
+            assert D1.old_poly_ring(x).set_domain(D2) == D2.old_poly_ring(x)
+            assert D1.old_poly_ring(x, y).set_domain(D2) == D2.old_poly_ring(x, y)
+            assert D1.old_frac_field(x).set_domain(D2) == D2.old_frac_field(x)
+            assert D1.old_frac_field(x, y).set_domain(D2) == D2.old_frac_field(x, y)
+
+
+def test_Domain_is_Exact():
+    exact = [GF(5), ZZ, QQ, ALG, EX]
+    inexact = [RR, CC]
+    for D in exact + inexact:
+        for R in D, D[x], D.frac_field(x), D.old_poly_ring(x), D.old_frac_field(x):
+            if D in exact:
+                assert R.is_Exact is True
+            else:
+                assert R.is_Exact is False
+
+
 def test_Domain_get_exact():
     assert EX.get_exact() == EX
     assert ZZ.get_exact() == ZZ
     assert QQ.get_exact() == QQ
     assert RR.get_exact() == QQ
-    # XXX: This should also be like RR:
-    # assert CC.get_exact() == QQ_I
+    assert CC.get_exact() == QQ_I
     assert ALG.get_exact() == ALG
     assert ZZ[x].get_exact() == ZZ[x]
     assert QQ[x].get_exact() == QQ[x]
+    assert RR[x].get_exact() == QQ[x]
+    assert CC[x].get_exact() == QQ_I[x]
     assert ZZ[x, y].get_exact() == ZZ[x, y]
     assert QQ[x, y].get_exact() == QQ[x, y]
+    assert RR[x, y].get_exact() == QQ[x, y]
+    assert CC[x, y].get_exact() == QQ_I[x, y]
     assert ZZ.frac_field(x).get_exact() == ZZ.frac_field(x)
     assert QQ.frac_field(x).get_exact() == QQ.frac_field(x)
+    assert RR.frac_field(x).get_exact() == QQ.frac_field(x)
+    assert CC.frac_field(x).get_exact() == QQ_I.frac_field(x)
     assert ZZ.frac_field(x, y).get_exact() == ZZ.frac_field(x, y)
     assert QQ.frac_field(x, y).get_exact() == QQ.frac_field(x, y)
+    assert RR.frac_field(x, y).get_exact() == QQ.frac_field(x, y)
+    assert CC.frac_field(x, y).get_exact() == QQ_I.frac_field(x, y)
+    assert ZZ.old_poly_ring(x).get_exact() == ZZ.old_poly_ring(x)
+    assert QQ.old_poly_ring(x).get_exact() == QQ.old_poly_ring(x)
+    assert RR.old_poly_ring(x).get_exact() == QQ.old_poly_ring(x)
+    assert CC.old_poly_ring(x).get_exact() == QQ_I.old_poly_ring(x)
+    assert ZZ.old_poly_ring(x, y).get_exact() == ZZ.old_poly_ring(x, y)
+    assert QQ.old_poly_ring(x, y).get_exact() == QQ.old_poly_ring(x, y)
+    assert RR.old_poly_ring(x, y).get_exact() == QQ.old_poly_ring(x, y)
+    assert CC.old_poly_ring(x, y).get_exact() == QQ_I.old_poly_ring(x, y)
+    assert ZZ.old_frac_field(x).get_exact() == ZZ.old_frac_field(x)
+    assert QQ.old_frac_field(x).get_exact() == QQ.old_frac_field(x)
+    assert RR.old_frac_field(x).get_exact() == QQ.old_frac_field(x)
+    assert CC.old_frac_field(x).get_exact() == QQ_I.old_frac_field(x)
+    assert ZZ.old_frac_field(x, y).get_exact() == ZZ.old_frac_field(x, y)
+    assert QQ.old_frac_field(x, y).get_exact() == QQ.old_frac_field(x, y)
+    assert RR.old_frac_field(x, y).get_exact() == QQ.old_frac_field(x, y)
+    assert CC.old_frac_field(x, y).get_exact() == QQ_I.old_frac_field(x, y)
 
 
 def test_Domain_characteristic():
@@ -614,8 +662,8 @@ def check_element(e1, e2, K1, K2, K3):
 
     def check_domains(K1, K2):
         K3 = K1.unify(K2)
-        check_element(K3.convert_from( K1.one, K1),  K3.one, K1, K2, K3)
-        check_element(K3.convert_from( K2.one, K2),  K3.one, K1, K2, K3)
+        check_element(K3.convert_from(K1.one, K1),  K3.one,  K1, K2, K3)
+        check_element(K3.convert_from(K2.one, K2),  K3.one,  K1, K2, K3)
         check_element(K3.convert_from(K1.zero, K1), K3.zero, K1, K2, K3)
         check_element(K3.convert_from(K2.zero, K2), K3.zero, K1, K2, K3)
 
@@ -648,6 +696,11 @@ def composite_domains(K):
     assert CC.convert(ZZ_I(1, 2)) == CC(1, 2)
     assert CC.convert(QQ_I(1, 2)) == CC(1, 2)
 
+    assert QQ.convert_from(RR(0.5), RR) == QQ(1, 2)
+    assert RR.convert_from(QQ(1, 2), QQ) == RR(0.5)
+    assert QQ_I.convert_from(CC(0.5, 0.75), CC) == QQ_I(QQ(1, 2), QQ(3, 4))
+    assert CC.convert_from(QQ_I(QQ(1, 2), QQ(3, 4)), QQ_I) == CC(0.5, 0.75)
+
     K1 = QQ.frac_field(x)
     K2 = ZZ.frac_field(x)
     K3 = QQ[x]
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index 87dd4cb424..1c4dd02dc3 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -1899,6 +1899,20 @@ def test_solve_nonlinear_trans():
     assert nonlinsolve([x**2 - y**2/exp(x)], [x, y]) == soln4
 
 
+def test_nonlinsolve_issue_25182():
+    a1, b1, c1, ca, cb, cg = symbols('a1, b1, c1, ca, cb, cg')
+    eq1 = a1*a1 + b1*b1 - 2.*a1*b1*cg - c1*c1
+    eq2 = a1*a1 + c1*c1 - 2.*a1*c1*cb - b1*b1
+    eq3 = b1*b1 + c1*c1 - 2.*b1*c1*ca - a1*a1
+    assert nonlinsolve([eq1, eq2, eq3], [c1, cb, cg]) == FiniteSet(
+        (1.0*b1*ca - 1.0*sqrt(a1**2 + b1**2*ca**2 - b1**2),
+        -1.0*sqrt(a1**2 + b1**2*ca**2 - b1**2)/a1,
+        -1.0*b1*(ca - 1)*(ca + 1)/a1 + 1.0*ca*sqrt(a1**2 + b1**2*ca**2 - b1**2)/a1),
+        (1.0*b1*ca + 1.0*sqrt(a1**2 + b1**2*ca**2 - b1**2),
+        1.0*sqrt(a1**2 + b1**2*ca**2 - b1**2)/a1,
+        -1.0*b1*(ca - 1)*(ca + 1)/a1 - 1.0*ca*sqrt(a1**2 + b1**2*ca**2 - b1**2)/a1))
+
+
 def test_issue_14642():
     x = Symbol('x')
     n1 = 0.5*x**3+x**2+0.5+I #add I in the Polynomials