Incompatibilities in Sage mpmath cause downstream issues in SymPy

[asmeurer]

The Sage mpmath has some incompatibilities with non-sage mpmath, which causes some issues in SymPy, since mpmath always uses the Sage backend if Sage is installed.

You can get a full idea of the issues by running the SymPy tests in Sage (import sympy;sympy.test()). A few failures are also due to some gmpy2 issues (gmpy2 types fail on Sage integers).

The issues I've noticed so far are:

• The sage mpmath types do not subclass from mpmath.ctx_mp.mpnumeric, which is the base class of mpf and mpc according to its own docstring.

• The sage mpf type does not have a context attribute.

Another issue, reported in sympy/sympy#19690:

import sympy.physics
from sympy.physics.units import meter
import copy
copy.deepcopy(3.7*meter)
[...]
File "sage/rings/integer.pyx", line 7143, in sage.rings.integer.mpz_set_str_python (build/cythonized/sage/rings/inte
ger.c:43554)
TypeError: unable to convert '0xecccccccccccd' to an integer

(fixed in sympy/sympy#21996, merged for inclusion in Sympy 1.9)

References:

• Install Sage on Travis using conda-forge sympy/sympy#14703
• doctesting sage-9.3.beta failures cschwan/sage-on-gentoo#628
• Enable sage backend by default only if SAGE_ROOT is set mpmath/mpmath#466
• sage-devel, 2021-03, pointy python question related to sage
• sage-devel, 2021-03, Running MathicsSession inside Sage
• Explicitly set MPMATH_SAGE env variable #32799: Explicitly set MPMATH_SAGE env variable

CC: @antonio-rojas @kiwifb @fredrik-johansson @isuruf @nbruin @saraedum @slel @soehms @strogdon @videlec

Component: distribution

Keywords: mpmath

Issue created by migration from https://trac.sagemath.org/ticket/25445

[slel]

Changed keywords from none to mpmath

[kiwifb]

comment:2

The content of the description echoes some of what I said on sage-devel, we are probably in this situation because sage misuses mpmath, or uses some parts that should be private.

[dimpase]

comment:7

mpmath update on #31564

[mkoeppe]

comment:8

Moving to 9.4, as 9.3 has been released.

[slel]

comment:13

See also #32799: Explicitly set MPMATH_SAGE env variable.

[sheerluck]

Replying to @asmeurer:

(fixed in sympy/sympy#21996, merged for inclusion in Sympy 1.9)

Yes, commit 9737261 is in master branch, but not in 1.9 branch and not in released sympy-1.9.tar.gz so deepcopy(3.7*meter) not gonna work :(

[sheerluck]

comment:16

Speaking about import sympy;sympy.test(), there are 75 failed tests, most of them
failed due to an incomplete int_types in mpmath:

$ ipython

In [1]: import mpmath.libmp.backend as mb

In [2]: mb.int_types
Out[2]: (int, mpz)

and

$ sage -q
sage: import mpmath.libmp.backend as mb
sage: mb.int_types
(<class 'int'>, <class 'sage.rings.integer.Integer'>)

But when mpmath is patched this way:

--- a/mpmath/libmp/backend.py
+++ b/mpmath/libmp/backend.py
@@ -76,6 +76,7 @@
         if gmpy.version() >= '1.03':
             BACKEND = 'gmpy'
             MPZ = gmpy.mpz
+            GMPY_MPZ = gmpy.mpz
     except:
         pass
 
@@ -111,5 +112,7 @@
 except NameError:
     if BACKEND == 'python':
         int_types = (int,)
+    elif BACKEND == 'sage':
+        int_types = (int, GMPY_MPZ, MPZ_TYPE)
     else:
         int_types = (int, MPZ_TYPE)

we have

$ sage -q
sage: import mpmath.libmp.backend as mb
sage: mb.int_types
(<class 'int'>, <class 'mpz'>, <class 'sage.rings.integer.Integer'>)

and number of failed tests went from 75 to 30.

When sympy is patched this way:

--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -4469,6 +4469,14 @@
 _sympy_converter[mpnumeric] = sympify_mpmath


+def sympify_sage(x):
+    from mpmath import mp
+    from sage.libs.mpmath.utils import sage_to_mpmath
+    return Expr._from_mpmath(sage_to_mpmath(x, mp.prec), mp.prec)
+
+_sympy_converter["sage"] = sympify_sage
+
+
 def sympify_complex(a):
     real, imag = list(map(sympify, (a.real, a.imag)))
     return real + S.ImaginaryUnit*imag
--- a/sympy/core/sympify.py
+++ b/sympy/core/sympify.py
@@ -382,6 +382,11 @@
         if conv is not None:
             return conv(a)

+    if hasattr(a, "_mpf_") or hasattr(a, "_mpc_"):
+        conv = _sympy_converter.get("sage")
+        if conv is not None:
+            return conv(a)
+
     if cls is type(None):
         if strict:
             raise SympifyError(a)
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -2919,7 +2919,8 @@
                 f[i] = fi.lhs - fi.rhs
         f = Matrix(f).T
     if iterable(x0):
-        x0 = list(x0)
+        if not hasattr(x0, "imag"):
+            x0 = list(x0)
     if not isinstance(f, Matrix):
         # assume it's a SymPy expression
         if isinstance(f, Eq):

number of failed tests went from 30 to 13.

We can fix "TypeError: '>=' not supported between instances of 'mpz' and
'sage.libs.mpmath.ext_main.mpf'" with something like this

--- a/sage/libs/mpmath/ext_main.pyx
+++ b/sage/libs/mpmath/ext_main.pyx
@@ -46,6 +46,8 @@
 import mpmath.function_docs as function_docs
 from mpmath.libmp import to_str
 from mpmath.libmp import repr_dps, prec_to_dps, dps_to_prec
+from mpmath.libmp import int_types, mpf_lt, mpf_gt, mpf_le, mpf_ge, mpf_eq
+from mpmath.libmp.libmpf import from_int

 DEF OP_ADD = 0
 DEF OP_SUB = 1
@@ -2101,6 +2103,13 @@
         MPF_normalize(&r.value, global_opts)
         return r

+    def __lt__(s, t): return s._cmp(t, mpf_lt, OP_LT)
+    def __gt__(s, t): return s._cmp(t, mpf_gt, OP_GT)
+    def __le__(s, t): return s._cmp(t, mpf_le, OP_LE)
+    def __ge__(s, t): return s._cmp(t, mpf_ge, OP_GE)
+    def __ne__(s, t): return not s.__eq__(t)
+    def __eq__(s, t): return s._cmp(t, mpf_eq, OP_EQ)
+
     def __abs__(s):
         """
         Computes the absolute value, rounded to the current
@@ -2128,7 +2137,7 @@
         MPF_sqrt(&r.value, &s.value, global_opts)
         return r

-    def __richcmp__(self, other, int op):
+    def _cmp(s, t, func, int op):
         """
         Compares numbers ::

@@ -2139,8 +2148,15 @@
             True
             sage: mpf(3) == 4
             False
+            sage: from gmpy2 import fac
+            sage: mpf(120) == fac(5)
+            True
+            sage: fac(5) >= mpf(3)
+            True
         """
-        return binop(OP_RICHCMP+op, self, other, global_opts)
+        if isinstance(t, int_types):
+            return func(s._mpf_, from_int(t))
+        return binop(op, s, t, global_opts)

Sadly that huge patch fixes only one failed test :)

So we left with 12 test, two of them (in test_holonomic.py) can be fixed using
mpmath.almosteq instead of comparing strings, one test in test_pickling.py
is about SymPyDeprecationWarning, one test in test_lambdify.py is just
"TypeError vs NotImplementedError", and all other failures are just
AssertionErrors and they actually work if you just try it:

sage: import mpmath
sage: from mpmath import mpc
sage: assert mpmath.besselj(2, 1 + 1j).ae(mpc("0.04157988694396212", "0.24739764151330632"))
sage: from sympy.functions.elementary.exponential import exp
sage: assert exp(pi * sqrt(163)).evalf(50).num.ae(262537412640768744)
sage: import sympy
....: import mpmath
....: from sympy.utilities.lambdify import lambdify
....: 
....: f_ = lambdify([x], sympy.LambertW(x, -1), modules="scipy")
sage: assert f_(0.1) == mpmath.lambertw(0.1, -1)
sage: from sympy.functions.special.gamma_functions import uppergamma, lowergamma
sage: expr1 = lambdify((x, y), uppergamma(x, y), "mpmath")(1, 2)
sage: assert expr1 == uppergamma(1, 2).evalf()
sage: from sympy.core.symbol import symbols
sage: from sympy.solvers import nsolve
sage: x, y = symbols("x y")
sage: nsolve(x**2 + 2, 1j)
-5.12588859584049e-25 + 1.4142135623731*I

The only test I can do nothing about is test_P25 in test_wester.py
Oh, and also that 3.7*meter too:

sage: import sympy.physics
....: from sympy.physics.units import meter
....: import copy 
....:
....: copy.deepcopy(3.7 * meter)
---------------------------------------------------------------------------
NotImplementedError: conversion to SageMath is not implemented

All References in this ticket are already fixed even without those three patches earlier in this comment. And this ticket would have been closed as solved if it wasn't for the test_P25 (and 3.7*meter)

[mkoeppe]

Can we solve this problem by vendoring a copy of mpmath for the internal uses within Sage?

[mkoeppe]

... done in #35095

• https://github.com/mkoeppe/sage/tree/t/32432/modularization_of_sagelib__break_out_a_separate_package_sagemath_polyhedra/pkgs/sagemath-mpmath
• https://github.com/mkoeppe/sage/tree/t/32432/modularization_of_sagelib__break_out_a_separate_package_sagemath_polyhedra/src/sage/libs/mpmath

[jhpalmieri]

Could this be the cause of the problem reported at https://ask.sagemath.org/question/73706/sages-python-not-giving-the-expected-result/ ?

[skirpichev]

Also, cython implementation of mpmath types in sage seems to be outdated (there are test failures in the mpmath test suite for 1.3.0): #36447 (comment)