Repository: sympy/sympy

Issue:
Simplify of small imaginary number yields 0
Using the most recent release (0.7.6) the following expression yields 0:  

    from sympy import *
    (I*Number(1.)*Number(10)**Number(-14)).simplify()

Using the absolute value this does not happen.

     (Number(1.)*Number(10)**Number(-14)).simplify()


Discussion:

Patch:
diff --git a/sympy/core/mod.py b/sympy/core/mod.py
index 23bdd33399..cefbf8c913 100644
--- a/sympy/core/mod.py
+++ b/sympy/core/mod.py
@@ -84,7 +84,7 @@ def doit(p, q):
 
         # extract gcd; any further simplification should be done by the user
         G = gcd(p, q)
-        if G is not S.One:
+        if G != 1:
             p, q = [
                 gcd_terms(i/G, clear=False, fraction=False) for i in (p, q)]
         pwas, qwas = p, q
diff --git a/sympy/polys/constructor.py b/sympy/polys/constructor.py
index bb86da60fa..409e0c2b8f 100644
--- a/sympy/polys/constructor.py
+++ b/sympy/polys/constructor.py
@@ -6,6 +6,7 @@
 from sympy.polys.polyoptions import build_options
 from sympy.polys.polyerrors import GeneratorsNeeded
 from sympy.polys.domains import ZZ, QQ, RR, EX
+from sympy.polys.domains.realfield import RealField
 from sympy.utilities import public
 from sympy.core import sympify
 
@@ -44,7 +45,10 @@ def _construct_simple(coeffs, opt):
         domain, result = _construct_algebraic(coeffs, opt)
     else:
         if reals:
-            domain = RR
+            # Use the maximum precision of all coefficients for the RR's
+            # precision
+            max_prec = max([c._prec for c in coeffs])
+            domain = RealField(prec=max_prec)
         else:
             if opt.field or rationals:
                 domain = QQ
diff --git a/sympy/polys/domains/mpelements.py b/sympy/polys/domains/mpelements.py
index d44a6f38c7..e5a89220e2 100644
--- a/sympy/polys/domains/mpelements.py
+++ b/sympy/polys/domains/mpelements.py
@@ -1,4 +1,4 @@
-"""Real and complex elements with built-in truncation. """
+"""Real and complex elements. """
 
 from __future__ import print_function, division
 
@@ -19,13 +19,7 @@ class RealElement(_mpf, DomainElement):
     __slots__ = ['__mpf__']
 
     def _set_mpf(self, val):
-        prec, rounding = self.context._prec_rounding
-        tol = self.context.tol
-
-        if mpf_lt(mpf_abs(val, prec, rounding), tol):
-            self.__mpf__ = fzero
-        else:
-            self.__mpf__ = val
+        self.__mpf__ = val
 
     _mpf_ = property(lambda self: self.__mpf__, _set_mpf)
 
@@ -39,20 +33,7 @@ class ComplexElement(_mpc, DomainElement):
     __slots__ = ['__mpc__']
 
     def _set_mpc(self, val):
-        prec, rounding = self.context._prec_rounding
-        tol = self.context.tol
-
-        # norm = mpc_abs(val, prec, rounding)
-        # tol = mpf_max(tol, mpf_mul(norm, tol))
-
-        re, im = val
-
-        if mpf_lt(mpf_abs(re, prec, rounding), tol):
-            re = fzero
-        if mpf_lt(mpf_abs(im, prec, rounding), tol):
-            im = fzero
-
-        self.__mpc__ = (re, im)
+        self.__mpc__ = val
 
     _mpc_ = property(lambda self: self.__mpc__, _set_mpc)
 


Test patch:
diff --git a/sympy/core/tests/test_arit.py b/sympy/core/tests/test_arit.py
index eaff343528..222c732c24 100644
--- a/sympy/core/tests/test_arit.py
+++ b/sympy/core/tests/test_arit.py
@@ -1553,7 +1553,14 @@ def test_Mod():
     assert (x + 3) % 1 == Mod(x, 1)
     assert (x + 3.0) % 1 == Mod(1.*x, 1)
     assert (x - S(33)/10) % 1 == Mod(x + S(7)/10, 1)
-    assert str(Mod(.6*x + y, .3*y)) == str(Mod(0.1*y + 0.6*x, 0.3*y))
+
+    a = Mod(.6*x + y, .3*y)
+    b = Mod(0.1*y + 0.6*x, 0.3*y)
+    # Test that a, b are equal, with 1e-14 accuracy in coefficients
+    eps = 1e-14
+    assert abs((a.args[0] - b.args[0]).subs({x: 1, y: 1})) < eps
+    assert abs((a.args[1] - b.args[1]).subs({x: 1, y: 1})) < eps
+
     assert (x + 1) % x == 1 % x
     assert (x + y) % x == y % x
     assert (x + y + 2) % x == (y + 2) % x
@@ -1564,7 +1571,7 @@ def test_Mod():
     assert (-3*x) % (-2*y) == -Mod(3*x, 2*y)
     assert (.6*pi) % (.3*x*pi) == 0.3*pi*Mod(2, x)
     assert (.6*pi) % (.31*x*pi) == pi*Mod(0.6, 0.31*x)
-    assert (6*pi) % (.3*x*pi) == pi*Mod(6, 0.3*x)
+    assert (6*pi) % (.3*x*pi) == 0.3*pi*Mod(20, x)
     assert (6*pi) % (.31*x*pi) == pi*Mod(6, 0.31*x)
     assert (6*pi) % (.42*x*pi) == pi*Mod(6, 0.42*x)
     assert (12*x) % (2*y) == 2*Mod(6*x, y)
diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py
index 650c24b306..9b0575bdab 100644
--- a/sympy/integrals/tests/test_integrals.py
+++ b/sympy/integrals/tests/test_integrals.py
@@ -11,7 +11,7 @@
 from sympy.physics import units
 from sympy.core.compatibility import range
 from sympy.utilities.pytest import XFAIL, raises, slow
-
+from sympy.utilities.randtest import verify_numerically
 
 
 x, y, a, t, x_1, x_2, z, s = symbols('x y a t x_1 x_2 z s')
@@ -1026,12 +1026,9 @@ def test_risch_option():
     # TODO: How to test risch=False?
 
 def test_issue_6828():
-    # TODO: Currently `h' is the result (all three are equivalent). Improve
-    # simplify() to find the form with simplest real coefficients.
     f = 1/(1.08*x**2 - 4.3)
-    g = 300.0/(324.0*x**2 - 1290.0)
-    h = 0.925925925925926/(1.0*x**2 - 3.98148148148148)
-    assert integrate(f, x).diff(x).simplify().equals(f) is True
+    g = integrate(f, x).diff(x)
+    assert verify_numerically(f, g, tol=1e-12)
 
 @XFAIL
 def test_integrate_Piecewise_rational_over_reals():
diff --git a/sympy/physics/optics/tests/test_medium.py b/sympy/physics/optics/tests/test_medium.py
index 66ec8664cf..e574e719b5 100644
--- a/sympy/physics/optics/tests/test_medium.py
+++ b/sympy/physics/optics/tests/test_medium.py
@@ -29,8 +29,11 @@ def test_medium():
     assert m4.refractive_index < m1.refractive_index
     assert m4 < m1
     m5 = Medium('m5', permittivity=710*10**(-12)*s**4*A**2/(m**3*kg), n=1.33)
-    assert simplify(m5.intrinsic_impedance - 6.24845417765552*kg*m**2/(A**2*s**3)) == 0
+    assert abs(m5.intrinsic_impedance - 6.24845417765552*kg*m**2/(A**2*s**3)) \
+                < 1e-12*kg*m**2/(A**2*s**3)
     assert abs(m5.speed - 225407863.157895*m/s) < 1e-6*m/s
-    assert simplify(m5.refractive_index - 1.33000000000000) == 0
-    assert simplify(m5.permittivity - 7.1e-10*A**2*s**4/(kg*m**3)) == 0
-    assert simplify(m5.permeability - 2.77206575232851e-8*kg*m/(A**2*s**2)) == 0
+    assert abs(m5.refractive_index - 1.33000000000000) < 1e-12
+    assert abs(m5.permittivity - 7.1e-10*A**2*s**4/(kg*m**3)) \
+                < 1e-20*A**2*s**4/(kg*m**3)
+    assert abs(m5.permeability - 2.77206575232851e-8*kg*m/(A**2*s**2)) \
+                < 1e-20*kg*m/(A**2*s**2)
diff --git a/sympy/polys/domains/tests/test_domains.py b/sympy/polys/domains/tests/test_domains.py
index c71ee686d4..6efd1575cc 100644
--- a/sympy/polys/domains/tests/test_domains.py
+++ b/sympy/polys/domains/tests/test_domains.py
@@ -1,9 +1,10 @@
 """Tests for classes defining properties of ground domains, e.g. ZZ, QQ, ZZ[x] ... """
 
-from sympy import S, sqrt, sin, oo, Poly
+from sympy import S, sqrt, sin, oo, Poly, Float
 from sympy.abc import x, y, z
 
 from sympy.polys.domains import ZZ, QQ, RR, CC, FF, GF, EX
+from sympy.polys.domains.realfield import RealField
 
 from sympy.polys.rings import ring
 from sympy.polys.fields import field
@@ -732,3 +733,38 @@ def test_ModularInteger():
 def test_QQ_int():
     assert int(QQ(2**2000, 3**1250)) == 455431
     assert int(QQ(2**100, 3)) == 422550200076076467165567735125
+
+def test_RR_double():
+    assert RR(3.14) > 1e-50
+    assert RR(1e-13) > 1e-50
+    assert RR(1e-14) > 1e-50
+    assert RR(1e-15) > 1e-50
+    assert RR(1e-20) > 1e-50
+    assert RR(1e-40) > 1e-50
+
+def test_RR_Float():
+    f1 = Float("1.01")
+    f2 = Float("1.0000000000000000000001")
+    assert f1._prec == 53
+    assert f2._prec == 80
+    assert RR(f1)-1 > 1e-50
+    assert RR(f2)-1 < 1e-50 # RR's precision is lower than f2's
+
+    RR2 = RealField(prec=f2._prec)
+    assert RR2(f1)-1 > 1e-50
+    assert RR2(f2)-1 > 1e-50 # RR's precision is equal to f2's
+
+def test_CC_double():
+    assert CC(3.14).real > 1e-50
+    assert CC(1e-13).real > 1e-50
+    assert CC(1e-14).real > 1e-50
+    assert CC(1e-15).real > 1e-50
+    assert CC(1e-20).real > 1e-50
+    assert CC(1e-40).real > 1e-50
+
+    assert CC(3.14j).imag > 1e-50
+    assert CC(1e-13j).imag > 1e-50
+    assert CC(1e-14j).imag > 1e-50
+    assert CC(1e-15j).imag > 1e-50
+    assert CC(1e-20j).imag > 1e-50
+    assert CC(1e-40j).imag > 1e-50
diff --git a/sympy/polys/tests/test_constructor.py b/sympy/polys/tests/test_constructor.py
index 62466e27a5..5238ec52af 100644
--- a/sympy/polys/tests/test_constructor.py
+++ b/sympy/polys/tests/test_constructor.py
@@ -2,8 +2,9 @@
 
 from sympy.polys.constructor import construct_domain
 from sympy.polys.domains import ZZ, QQ, RR, EX
+from sympy.polys.domains.realfield import RealField
 
-from sympy import S, sqrt, sin
+from sympy import S, sqrt, sin, Float
 from sympy.abc import x, y
 
 def test_construct_domain():
@@ -14,7 +15,9 @@ def test_construct_domain():
     assert construct_domain([S(1), S(2), S(3)], field=True) == (QQ, [QQ(1), QQ(2), QQ(3)])
 
     assert construct_domain([S(1)/2, S(2)]) == (QQ, [QQ(1, 2), QQ(2)])
-    assert construct_domain([3.14, 1, S(1)/2]) == (RR, [RR(3.14), RR(1.0), RR(0.5)])
+    result = construct_domain([3.14, 1, S(1)/2])
+    assert isinstance(result[0], RealField)
+    assert result[1] == [RR(3.14), RR(1.0), RR(0.5)]
 
     assert construct_domain([3.14, sqrt(2)], extension=None) == (EX, [EX(3.14), EX(sqrt(2))])
     assert construct_domain([3.14, sqrt(2)], extension=True) == (EX, [EX(3.14), EX(sqrt(2))])
@@ -101,3 +104,20 @@ def test_composite_option():
 
     assert construct_domain({(1, 0): y}, composite=False) == \
         (EX, {(1, 0): EX(y)})
+
+def test_precision():
+    f1 = Float("1.01")
+    f2 = Float("1.0000000000000000000001")
+    for x in [1, 1e-2, 1e-6, 1e-13, 1e-14, 1e-16, 1e-20, 1e-100, 1e-300,
+            f1, f2]:
+        result = construct_domain([x])
+        y = float(result[1][0])
+        assert abs(x - y) / x < 1e-14  # Test relative accuracy
+
+    result = construct_domain([f1])
+    y = result[1][0]
+    assert y-1 > 1e-50
+
+    result = construct_domain([f2])
+    y = result[1][0]
+    assert y-1 > 1e-50
diff --git a/sympy/polys/tests/test_polytools.py b/sympy/polys/tests/test_polytools.py
index 56ccb0c726..38deab20d5 100644
--- a/sympy/polys/tests/test_polytools.py
+++ b/sympy/polys/tests/test_polytools.py
@@ -47,6 +47,7 @@
 
 from sympy.polys.fields import field
 from sympy.polys.domains import FF, ZZ, QQ, RR, EX
+from sympy.polys.domains.realfield import RealField
 from sympy.polys.orderings import lex, grlex, grevlex
 
 from sympy import (
@@ -346,7 +347,7 @@ def test_Poly__new__():
         Poly(3*x**5 + 65536*x**4 + x**3 + 65536*x** 2 + 1, x,
              modulus=65537, symmetric=False)
 
-    assert Poly(x**2 + x + 1.0).get_domain() == RR
+    assert isinstance(Poly(x**2 + x + 1.0).get_domain(), RealField)
 
 
 def test_Poly__args():
@@ -560,7 +561,7 @@ def test_Poly_get_domain():
     raises(CoercionFailed, lambda: Poly(x/2, domain='ZZ'))
     assert Poly(x/2, domain='QQ').get_domain() == QQ
 
-    assert Poly(0.2*x).get_domain() == RR
+    assert isinstance(Poly(0.2*x).get_domain(), RealField)
 
 
 def test_Poly_set_domain():
diff --git a/sympy/simplify/tests/test_simplify.py b/sympy/simplify/tests/test_simplify.py
index 4705bdd359..22deb4f78f 100644
--- a/sympy/simplify/tests/test_simplify.py
+++ b/sympy/simplify/tests/test_simplify.py
@@ -229,7 +229,7 @@ def test_issue_3210():
 
 def test_issue_7263():
     assert abs((simplify(30.8**2 - 82.5**2 * sin(rad(11.6))**2)).evalf() - \
-            673.447451402970) < 1e-15
+            673.447451402970) < 1e-12
 
 
 def test_trigsimp_issues():
@@ -1914,6 +1914,34 @@ def test_inequality_no_auto_simplify():
     assert e == Lt(lhs, rhs, evaluate=False)
     assert simplify(e)
 
+def test_issue_9398():
+    from sympy import Number, cancel
+    assert cancel(1e-14) != 0
+    assert cancel(1e-14*I) != 0
+
+    assert simplify(1e-14) != 0
+    assert simplify(1e-14*I) != 0
+
+    assert (I*Number(1.)*Number(10)**Number(-14)).simplify() != 0
+
+    assert cancel(1e-20) != 0
+    assert cancel(1e-20*I) != 0
+
+    assert simplify(1e-20) != 0
+    assert simplify(1e-20*I) != 0
+
+    assert cancel(1e-100) != 0
+    assert cancel(1e-100*I) != 0
+
+    assert simplify(1e-100) != 0
+    assert simplify(1e-100*I) != 0
+
+    f = Float("1e-1000")
+    assert cancel(f) != 0
+    assert cancel(f*I) != 0
+
+    assert simplify(f) != 0
+    assert simplify(f*I) != 0
 
 def test_issue_9324_simplify():
     M = MatrixSymbol('M', 10, 10)
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index 09fa1537ae..cfd5ca0c74 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -16,7 +16,7 @@
 from sympy.physics.units import cm
 from sympy.polys.rootoftools import RootOf
 
-from sympy.utilities.pytest import slow, XFAIL, raises, skip, ON_TRAVIS
+from sympy.utilities.pytest import slow, XFAIL, SKIP, raises, skip, ON_TRAVIS
 from sympy.utilities.randtest import verify_numerically as tn
 
 from sympy.abc import a, b, c, d, k, h, p, x, y, z, t, q, m
@@ -743,8 +743,8 @@ def test_issue_5335():
     assert len(solve(eqs, sym, manual=True, minimal=True, simplify=False)) == 2
 
 
-@XFAIL
-def test_issue_5335_float():
+@SKIP("Hangs")
+def _test_issue_5335_float():
     # gives ZeroDivisionError: polynomial division
     lam, a0, conc = symbols('lam a0 conc')
     eqs = [lam + 2*y - a0*(1 - x/2)*x - 0.005*x/2*x,