diff --git a/.github/workflows/ci-meson.yml b/.github/workflows/ci-meson.yml index 98523d35951..730fbfd6d3c 100644 --- a/.github/workflows/ci-meson.yml +++ b/.github/workflows/ci-meson.yml @@ -142,20 +142,6 @@ jobs: # Use --no-deps and pip check below to verify that all necessary dependencies are installed via conda pip install --no-build-isolation --no-deps --config-settings=builddir=builddir ${{ matrix.editable && '--editable' || '' }} . -v - - name: Check update-meson - # this step must be after build, because meson.build creates a number of __init__.py files - # that is needed to make tools/update-meson.py run correctly - shell: bash -l {0} - if: matrix.tests == 'all' && matrix.python == '3.12' - run: | - python tools/update-meson.py - if ! ./tools/test-git-no-uncommitted-changes; then - git add --intent-to-add . # also show newly created files in git diff - git status - git diff - false - fi - - name: Verify dependencies shell: bash -l {0} run: pip check diff --git a/build/pkgs/sagelib/dependencies b/build/pkgs/sagelib/dependencies index 662fad365b5..f1db5e88139 100644 --- a/build/pkgs/sagelib/dependencies +++ b/build/pkgs/sagelib/dependencies @@ -1,4 +1,4 @@ -FORCE $(SCRIPTS) $(BLAS) brial cliquer cypari cysignals cython ecl eclib ecm flint libgd gap givaro glpk gmpy2 gsl iml importlib_metadata importlib_resources jupyter_core lcalc lrcalc_python libbraiding libhomfly libpng linbox m4ri m4rie memory_allocator mpc mpfi mpfr $(MP_LIBRARY) ntl numpy pari pip pkgconfig planarity ppl pplpy primesieve primecount primecountpy $(PYTHON) requests rw singular symmetrica typing_extensions $(PCFILES) | $(PYTHON_TOOLCHAIN) meson_python $(PYTHON) pythran +FORCE $(SCRIPTS) $(BLAS) brial cliquer cypari cysignals cython ecl eclib ecm flint libgd gap givaro glpk gmpy2 gsl iml importlib_metadata importlib_resources jupyter_core lcalc lrcalc_python libbraiding libhomfly libpng linbox m4ri m4rie memory_allocator mpc mpfi mpfr $(MP_LIBRARY) mpmath ntl numpy pari pip pkgconfig planarity ppl pplpy primesieve primecount primecountpy $(PYTHON) requests rw singular symmetrica typing_extensions $(PCFILES) | $(PYTHON_TOOLCHAIN) meson_python $(PYTHON) pythran ---------- All lines of this file are ignored except the first. diff --git a/src/sage/env.py b/src/sage/env.py index c2a1f7e562d..e8d5d37fa0c 100644 --- a/src/sage/env.py +++ b/src/sage/env.py @@ -231,8 +231,9 @@ def var(key: str, *fallbacks: Optional[str], force: bool = False) -> Optional[st OPENMP_CFLAGS = var("OPENMP_CFLAGS", "") OPENMP_CXXFLAGS = var("OPENMP_CXXFLAGS", "") -# Make sure mpmath uses Sage types -os.environ['MPMATH_SAGE'] = '1' +# Make sure that mpmath < 1.4 does not try to use Sage types +os.environ.pop('MPMATH_SAGE', None) +os.environ['MPMATH_NOSAGE'] = '1' # misc SAGE_BANNER = var("SAGE_BANNER", "") diff --git a/src/sage/features/mpmath13.py b/src/sage/features/mpmath13.py new file mode 100644 index 00000000000..3fc6066d2e9 --- /dev/null +++ b/src/sage/features/mpmath13.py @@ -0,0 +1,39 @@ +r""" +Feature for test for mpmath-1.3.x. Used to simultaneously +support v1.3 and v1.4 in the doctests. +""" +from . import Feature + + +class MpMath13(Feature): + r""" + A :class:`~sage.features.Feature`. + + EXAMPLES:: + + sage: from sage.features.mpmath13 import MpMath13 + sage: MpMath13() + Feature('mpmath13') + sage: MpMath13() is MpMath13() + True + """ + def __init__(self): + Feature.__init__(self, 'mpmath13') + + def _is_present(self): + r""" + Test the mpmath version. + + EXAMPLES:: + + sage: from sage.features.mpmath13 import MpMath13 + sage: MpMath()._is_present() # needs mpmath13 + FeatureTestResult('mpmath13', True) + + """ + from mpmath import __version__ as mpmath_version + return mpmath_version.startswith("1.3") + + +def all_features(): + return [MpMath13()] diff --git a/src/sage/libs/mpmath/ext_impl.pyx b/src/sage/libs/mpmath/ext_impl.pyx index 6c9d2e623a5..97396ab4542 100644 --- a/src/sage/libs/mpmath/ext_impl.pyx +++ b/src/sage/libs/mpmath/ext_impl.pyx @@ -1,3 +1,4 @@ +# sage.doctest: needs sage.libs.mpmath13 """ This module provides the core implementation of multiprecision floating-point arithmetic. Operations are done in-place. diff --git a/src/sage/libs/mpmath/ext_libmp.pyx b/src/sage/libs/mpmath/ext_libmp.pyx index 487a12d45bc..0f33de564ed 100644 --- a/src/sage/libs/mpmath/ext_libmp.pyx +++ b/src/sage/libs/mpmath/ext_libmp.pyx @@ -1,3 +1,4 @@ +# sage.doctest: needs sage.libs.mpmath13 """ Faster versions of some key functions in mpmath.libmp """ diff --git a/src/sage/libs/mpmath/ext_main.pyx b/src/sage/libs/mpmath/ext_main.pyx index f8c4420f2f4..45c848493a8 100644 --- a/src/sage/libs/mpmath/ext_main.pyx +++ b/src/sage/libs/mpmath/ext_main.pyx @@ -1,3 +1,4 @@ +# sage.doctest: needs sage.libs.mpmath13 """ mpmath floating-point numbers diff --git a/src/sage/libs/mpmath/meson.build b/src/sage/libs/mpmath/meson.build index a49cfaa64df..8ace5d259dd 100644 --- a/src/sage/libs/mpmath/meson.build +++ b/src/sage/libs/mpmath/meson.build @@ -1,31 +1,69 @@ -py.install_sources( - '__init__.py', - 'all.py', - 'ext_impl.pxd', - 'ext_impl.pyx', - 'ext_libmp.pyx', - 'ext_main.pxd', - 'ext_main.pyx', - 'utils.pxd', - 'utils.pyx', - subdir: 'sage/libs/mpmath', -) - -extension_data = { - 'ext_impl' : files('ext_impl.pyx'), - 'ext_libmp' : files('ext_libmp.pyx'), - 'ext_main' : files('ext_main.pyx'), - 'utils' : files('utils.pyx'), -} - -foreach name, pyx : extension_data - py.extension_module( - name, - sources: pyx, +# Determine the mpmath version at build time so that we can support +# both 1.3 and 1.4 simultaneously (for a while). +mpmath14 = run_command( + py, [ + '-c', + 'import mpmath; print(mpmath.__version__);' + ], + check: true, +).stdout().strip().startswith('1.4') + + +if mpmath14 + py.install_sources( + '__init__.py', + 'all.py', + 'utils.pxd', + 'utils.pyx', subdir: 'sage/libs/mpmath', - install: true, - include_directories: [inc_cpython, inc_ext, inc_rings], - dependencies: [py_dep, cysignals, gmp, mpfr], ) -endforeach + extension_data = { + 'utils' : files('utils.pyx'), + } + + foreach name, pyx : extension_data + py.extension_module( + name, + sources: pyx, + subdir: 'sage/libs/mpmath', + install: true, + include_directories: [inc_cpython, inc_ext, inc_rings], + dependencies: [py_dep, cysignals, gmp, gmpy2, mpfr], + ) + endforeach + +else + # mpmath-1.3.x + py.install_sources( + '__init__.py', + 'all.py', + 'ext_impl.pxd', + 'ext_impl.pyx', + 'ext_libmp.pyx', + 'ext_main.pxd', + 'ext_main.pyx', + 'utils.pxd', + 'utils13.pyx', + subdir: 'sage/libs/mpmath', + ) + + extension_data = { + 'ext_impl' : files('ext_impl.pyx'), + 'ext_libmp' : files('ext_libmp.pyx'), + 'ext_main' : files('ext_main.pyx'), + 'utils' : files('utils13.pyx'), + } + + foreach name, pyx : extension_data + py.extension_module( + name, + sources: pyx, + subdir: 'sage/libs/mpmath', + install: true, + include_directories: [inc_cpython, inc_ext, inc_rings], + dependencies: [py_dep, cysignals, gmp, mpfr], + ) + endforeach + +endif diff --git a/src/sage/libs/mpmath/utils.pyx b/src/sage/libs/mpmath/utils.pyx index 1acf14ea540..9209befc863 100644 --- a/src/sage/libs/mpmath/utils.pyx +++ b/src/sage/libs/mpmath/utils.pyx @@ -1,9 +1,9 @@ """ Utilities for Sage-mpmath interaction - -Also patches some mpmath functions for speed """ +cimport gmpy2 + from sage.ext.stdsage cimport PY_NEW from sage.rings.integer cimport Integer @@ -16,145 +16,7 @@ from sage.libs.gmp.all cimport * from sage.rings.real_mpfr cimport RealField -cpdef int bitcount(n) noexcept: - """ - Bitcount of a Sage Integer or Python int/long. - - EXAMPLES:: - - sage: from mpmath.libmp import bitcount - sage: bitcount(0) - 0 - sage: bitcount(1) - 1 - sage: bitcount(100) - 7 - sage: bitcount(-100) - 7 - sage: bitcount(2r) - 2 - sage: bitcount(2L) - 2 - """ - cdef Integer m - if isinstance(n, Integer): - m = n - else: - m = Integer(n) - if mpz_sgn(m.value) == 0: - return 0 - return mpz_sizeinbase(m.value, 2) - -cpdef isqrt(n): - """ - Square root (rounded to floor) of a Sage Integer or Python int/long. - The result is a Sage Integer. - - EXAMPLES:: - - sage: from mpmath.libmp import isqrt - sage: isqrt(0) - 0 - sage: isqrt(100) - 10 - sage: isqrt(10) - 3 - sage: isqrt(10r) - 3 - sage: isqrt(10L) - 3 - """ - cdef Integer m, y - if isinstance(n, Integer): - m = n - else: - m = Integer(n) - if mpz_sgn(m.value) < 0: - raise ValueError("square root of negative integer not defined.") - y = PY_NEW(Integer) - mpz_sqrt(y.value, m.value) - return y - -cpdef from_man_exp(man, exp, long prec=0, str rnd='d'): - """ - Create normalized mpf value tuple from mantissa and exponent. - - With prec > 0, rounds the result in the desired direction - if necessary. - - EXAMPLES:: - - sage: from mpmath.libmp import from_man_exp - sage: from_man_exp(-6, -1) - (1, 3, 0, 2) - sage: from_man_exp(-6, -1, 1, 'd') - (1, 1, 1, 1) - sage: from_man_exp(-6, -1, 1, 'u') - (1, 1, 2, 1) - """ - cdef Integer res - cdef long bc - res = Integer(man) - bc = mpz_sizeinbase(res.value, 2) - if not prec: - prec = bc - if mpz_sgn(res.value) < 0: - mpz_neg(res.value, res.value) - return normalize(1, res, exp, bc, prec, rnd) - else: - return normalize(0, res, exp, bc, prec, rnd) - -cpdef normalize(long sign, Integer man, exp, long bc, long prec, str rnd): - """ - Create normalized mpf value tuple from full list of components. - - EXAMPLES:: - - sage: from mpmath.libmp import normalize - sage: normalize(0, 4, 5, 3, 53, 'n') - (0, 1, 7, 1) - """ - cdef long shift - cdef Integer res - cdef unsigned long trail - if mpz_sgn(man.value) == 0: - from mpmath.libmp import fzero - return fzero - if bc <= prec and mpz_odd_p(man.value): - return (sign, man, exp, bc) - shift = bc - prec - res = PY_NEW(Integer) - if shift > 0: - if rnd == 'n': - if mpz_tstbit(man.value, shift-1) and (mpz_tstbit(man.value, shift) - or (mpz_scan1(man.value, 0) < (shift-1))): - mpz_cdiv_q_2exp(res.value, man.value, shift) - else: - mpz_fdiv_q_2exp(res.value, man.value, shift) - elif rnd == 'd': - mpz_fdiv_q_2exp(res.value, man.value, shift) - elif rnd == 'f': - if sign: - mpz_cdiv_q_2exp(res.value, man.value, shift) - else: - mpz_fdiv_q_2exp(res.value, man.value, shift) - elif rnd == 'c': - if sign: - mpz_fdiv_q_2exp(res.value, man.value, shift) - else: - mpz_cdiv_q_2exp(res.value, man.value, shift) - elif rnd == 'u': - mpz_cdiv_q_2exp(res.value, man.value, shift) - exp += shift - else: - mpz_set(res.value, man.value) - # Strip trailing bits - trail = mpz_scan1(res.value, 0) - if 0 < trail < bc: - mpz_tdiv_q_2exp(res.value, res.value, trail) - exp += trail - bc = mpz_sizeinbase(res.value, 2) - return (sign, res, int(exp), bc) +gmpy2.import_gmpy2() cdef mpfr_from_mpfval(mpfr_t res, tuple x): """ @@ -162,11 +24,11 @@ cdef mpfr_from_mpfval(mpfr_t res, tuple x): data tuple. """ cdef int sign - cdef Integer man + cdef gmpy2.mpz man cdef long exp sign, man, exp, _ = x if man: - mpfr_set_z(res, man.value, MPFR_RNDZ) + mpfr_set_z(res, man.z, MPFR_RNDZ) if sign: mpfr_neg(res, res, MPFR_RNDZ) mpfr_mul_2si(res, res, exp, MPFR_RNDZ) @@ -210,7 +72,7 @@ cdef mpfr_to_mpfval(mpfr_t value): mpz_tdiv_q_2exp(man.value, man.value, trailing) exp += trailing bc = mpz_sizeinbase(man.value, 2) - return (sign, man, int(exp), bc) + return (sign, man.__mpz__(), int(exp), bc) def mpmath_to_sage(x, prec): @@ -305,7 +167,7 @@ def sage_to_mpmath(x, prec): sage: print(a.sage_to_mpmath(1+pi, 53)) 4.14159265358979 sage: a.sage_to_mpmath(infinity, 53) - mpf('+inf') + mpf('inf') sage: a.sage_to_mpmath(-infinity, 53) mpf('-inf') sage: a.sage_to_mpmath(NaN, 53) @@ -413,7 +275,7 @@ def call(func, *args, **kwargs): Check that :issue:`11885` is fixed:: - sage: a.call(a.ei, 1.0r, parent=float) + sage: a.call(a.ei, 1.0r, parent=float) # abs tol 1e-15 1.8951178163559366 Check that :issue:`14984` is fixed:: diff --git a/src/sage/libs/mpmath/utils13.pyx b/src/sage/libs/mpmath/utils13.pyx new file mode 100644 index 00000000000..47c81cb2f97 --- /dev/null +++ b/src/sage/libs/mpmath/utils13.pyx @@ -0,0 +1,454 @@ +# sage.doctest: needs sage.libs.mpmath13 +""" +Utilities for Sage-mpmath interaction + +Also patches some mpmath functions for speed +""" + +from sage.ext.stdsage cimport PY_NEW + +from sage.rings.integer cimport Integer +from sage.rings.real_mpfr cimport RealNumber +from sage.rings.complex_mpfr cimport ComplexNumber +from sage.structure.element cimport Element + +from sage.libs.mpfr cimport * +from sage.libs.gmp.all cimport * + +from sage.rings.real_mpfr cimport RealField + +cpdef int bitcount(n) noexcept: + """ + Bitcount of a Sage Integer or Python int/long. + + EXAMPLES:: + + sage: from mpmath.libmp import bitcount + sage: bitcount(0) + 0 + sage: bitcount(1) + 1 + sage: bitcount(100) + 7 + sage: bitcount(-100) + 7 + sage: bitcount(2r) + 2 + sage: bitcount(2L) + 2 + """ + cdef Integer m + if isinstance(n, Integer): + m = n + else: + m = Integer(n) + if mpz_sgn(m.value) == 0: + return 0 + return mpz_sizeinbase(m.value, 2) + +cpdef isqrt(n): + """ + Square root (rounded to floor) of a Sage Integer or Python int/long. + The result is a Sage Integer. + + EXAMPLES:: + + sage: from mpmath.libmp import isqrt + sage: isqrt(0) + 0 + sage: isqrt(100) + 10 + sage: isqrt(10) + 3 + sage: isqrt(10r) + 3 + sage: isqrt(10L) + 3 + """ + cdef Integer m, y + if isinstance(n, Integer): + m = n + else: + m = Integer(n) + if mpz_sgn(m.value) < 0: + raise ValueError("square root of negative integer not defined.") + y = PY_NEW(Integer) + mpz_sqrt(y.value, m.value) + return y + +cpdef from_man_exp(man, exp, long prec=0, str rnd='d'): + """ + Create normalized mpf value tuple from mantissa and exponent. + + With prec > 0, rounds the result in the desired direction + if necessary. + + EXAMPLES:: + + sage: from mpmath.libmp import from_man_exp + sage: from_man_exp(-6, -1) + (1, 3, 0, 2) + sage: from_man_exp(-6, -1, 1, 'd') + (1, 1, 1, 1) + sage: from_man_exp(-6, -1, 1, 'u') + (1, 1, 2, 1) + """ + cdef Integer res + cdef long bc + res = Integer(man) + bc = mpz_sizeinbase(res.value, 2) + if not prec: + prec = bc + if mpz_sgn(res.value) < 0: + mpz_neg(res.value, res.value) + return normalize(1, res, exp, bc, prec, rnd) + else: + return normalize(0, res, exp, bc, prec, rnd) + +cpdef normalize(long sign, Integer man, exp, long bc, long prec, str rnd): + """ + Create normalized mpf value tuple from full list of components. + + EXAMPLES:: + + sage: from mpmath.libmp import normalize + sage: normalize(0, 4, 5, 3, 53, 'n') + (0, 1, 7, 1) + """ + cdef long shift + cdef Integer res + cdef unsigned long trail + if mpz_sgn(man.value) == 0: + from mpmath.libmp import fzero + return fzero + if bc <= prec and mpz_odd_p(man.value): + return (sign, man, exp, bc) + shift = bc - prec + res = PY_NEW(Integer) + if shift > 0: + if rnd == 'n': + if mpz_tstbit(man.value, shift-1) and (mpz_tstbit(man.value, shift) + or (mpz_scan1(man.value, 0) < (shift-1))): + mpz_cdiv_q_2exp(res.value, man.value, shift) + else: + mpz_fdiv_q_2exp(res.value, man.value, shift) + elif rnd == 'd': + mpz_fdiv_q_2exp(res.value, man.value, shift) + elif rnd == 'f': + if sign: + mpz_cdiv_q_2exp(res.value, man.value, shift) + else: + mpz_fdiv_q_2exp(res.value, man.value, shift) + elif rnd == 'c': + if sign: + mpz_fdiv_q_2exp(res.value, man.value, shift) + else: + mpz_cdiv_q_2exp(res.value, man.value, shift) + elif rnd == 'u': + mpz_cdiv_q_2exp(res.value, man.value, shift) + exp += shift + else: + mpz_set(res.value, man.value) + # Strip trailing bits + trail = mpz_scan1(res.value, 0) + if 0 < trail < bc: + mpz_tdiv_q_2exp(res.value, res.value, trail) + exp += trail + bc = mpz_sizeinbase(res.value, 2) + return (sign, res, int(exp), bc) + +cdef mpfr_from_mpfval(mpfr_t res, tuple x): + """ + Set value of an MPFR number (in place) to that of a given mpmath mpf + data tuple. + """ + cdef int sign + cdef Integer man + cdef long exp + sign, man, exp, _ = x + if man: + mpfr_set_z(res, man.value, MPFR_RNDZ) + if sign: + mpfr_neg(res, res, MPFR_RNDZ) + mpfr_mul_2si(res, res, exp, MPFR_RNDZ) + return + from mpmath.libmp import finf, fninf + if exp == 0: + mpfr_set_ui(res, 0, MPFR_RNDZ) + elif x == finf: + mpfr_set_inf(res, 1) + elif x == fninf: + mpfr_set_inf(res, -1) + else: + mpfr_set_nan(res) + +cdef mpfr_to_mpfval(mpfr_t value): + """ + Given an MPFR value, return an mpmath mpf data tuple representing + the same number. + """ + if mpfr_nan_p(value): + from mpmath.libmp import fnan + return fnan + if mpfr_inf_p(value): + from mpmath.libmp import finf, fninf + if mpfr_sgn(value) > 0: + return finf + else: + return fninf + if mpfr_sgn(value) == 0: + from mpmath.libmp import fzero + return fzero + sign = 0 + cdef Integer man = PY_NEW(Integer) + exp = mpfr_get_z_exp(man.value, value) + if mpz_sgn(man.value) < 0: + mpz_neg(man.value, man.value) + sign = 1 + cdef unsigned long trailing + trailing = mpz_scan1(man.value, 0) + if trailing: + mpz_tdiv_q_2exp(man.value, man.value, trailing) + exp += trailing + bc = mpz_sizeinbase(man.value, 2) + return (sign, man, int(exp), bc) + + +def mpmath_to_sage(x, prec): + """ + Convert any mpmath number (mpf or mpc) to a Sage RealNumber or + ComplexNumber of the given precision. + + EXAMPLES:: + + sage: import sage.libs.mpmath.all as a + sage: a.mpmath_to_sage(a.mpf('2.5'), 53) + 2.50000000000000 + sage: a.mpmath_to_sage(a.mpc('2.5','-3.5'), 53) + 2.50000000000000 - 3.50000000000000*I + sage: a.mpmath_to_sage(a.mpf('inf'), 53) + +infinity + sage: a.mpmath_to_sage(a.mpf('-inf'), 53) + -infinity + sage: a.mpmath_to_sage(a.mpf('nan'), 53) + NaN + sage: a.mpmath_to_sage(a.mpf('0'), 53) + 0.000000000000000 + + A real example:: + + sage: RealField(100)(pi) + 3.1415926535897932384626433833 + sage: t = RealField(100)(pi)._mpmath_(); t + mpf('3.1415926535897932') + sage: a.mpmath_to_sage(t, 100) + 3.1415926535897932384626433833 + + We can ask for more precision, but the result is undefined:: + + sage: a.mpmath_to_sage(t, 140) # random + 3.1415926535897932384626433832793333156440 + sage: ComplexField(140)(pi) + 3.1415926535897932384626433832795028841972 + + A complex example:: + + sage: ComplexField(100)([0, pi]) + 3.1415926535897932384626433833*I + sage: t = ComplexField(100)([0, pi])._mpmath_(); t + mpc(real='0.0', imag='3.1415926535897932') + sage: sage.libs.mpmath.all.mpmath_to_sage(t, 100) + 3.1415926535897932384626433833*I + + Again, we can ask for more precision, but the result is undefined:: + + sage: sage.libs.mpmath.all.mpmath_to_sage(t, 140) # random + 3.1415926535897932384626433832793333156440*I + sage: ComplexField(140)([0, pi]) + 3.1415926535897932384626433832795028841972*I + """ + cdef RealNumber y + cdef ComplexNumber z + if hasattr(x, "_mpf_"): + y = RealField(prec)() + mpfr_from_mpfval(y.value, x._mpf_) + return y + elif hasattr(x, "_mpc_"): + from sage.rings.complex_mpfr import ComplexField + z = ComplexField(prec)(0) + re, im = x._mpc_ + mpfr_from_mpfval(z.__re, re) + mpfr_from_mpfval(z.__im, im) + return z + else: + raise TypeError("cannot convert %r to Sage", x) + + +def sage_to_mpmath(x, prec): + """ + Convert any Sage number that can be coerced into a RealNumber + or ComplexNumber of the given precision into an mpmath mpf or mpc. + Integers are currently converted to int. + + Lists, tuples and dicts passed as input are converted + recursively. + + EXAMPLES:: + + sage: import sage.libs.mpmath.all as a + sage: a.mp.dps = 15 + sage: print(a.sage_to_mpmath(2/3, 53)) + 0.666666666666667 + sage: print(a.sage_to_mpmath(2./3, 53)) + 0.666666666666667 + sage: print(a.sage_to_mpmath(3+4*I, 53)) + (3.0 + 4.0j) + sage: print(a.sage_to_mpmath(1+pi, 53)) + 4.14159265358979 + sage: a.sage_to_mpmath(infinity, 53) + mpf('+inf') + sage: a.sage_to_mpmath(-infinity, 53) + mpf('-inf') + sage: a.sage_to_mpmath(NaN, 53) + mpf('nan') + sage: a.sage_to_mpmath(0, 53) + 0 + sage: a.sage_to_mpmath([0.5, 1.5], 53) + [mpf('0.5'), mpf('1.5')] + sage: a.sage_to_mpmath((0.5, 1.5), 53) + (mpf('0.5'), mpf('1.5')) + sage: a.sage_to_mpmath({'n':0.5}, 53) + {'n': mpf('0.5')} + """ + if isinstance(x, Element): + if isinstance(x, Integer): + return int(x) + try: + if isinstance(x, RealNumber): + return x._mpmath_() + else: + x = RealField(prec)(x) + return x._mpmath_() + except TypeError: + if isinstance(x, ComplexNumber): + return x._mpmath_() + else: + from sage.rings.complex_mpfr import ComplexField + x = ComplexField(prec)(x) + return x._mpmath_() + if isinstance(x, (tuple, list)): + return type(x)([sage_to_mpmath(v, prec) for v in x]) + if isinstance(x, dict): + return {k: sage_to_mpmath(v, prec) for k, v in x.items()} + return x + + +def call(func, *args, **kwargs): + """ + Call an mpmath function with Sage objects as inputs and + convert the result back to a Sage real or complex number. + + By default, a RealNumber or ComplexNumber with the current + working precision of mpmath (mpmath.mp.prec) will be returned. + + If prec=n is passed among the keyword arguments, the temporary + working precision will be set to n and the result will also + have this precision. + + If parent=P is passed, P.prec() will be used as working + precision and the result will be coerced to P (or the + corresponding complex field if necessary). + + Arguments should be Sage objects that can be coerced into RealField + or ComplexField elements. Arguments may also be tuples, lists or + dicts (which are converted recursively), or any type that mpmath + understands natively (e.g. Python floats, strings for options). + + EXAMPLES:: + + sage: import sage.libs.mpmath.all as a + sage: a.mp.prec = 53 + sage: a.call(a.erf, 3+4*I) + -120.186991395079 - 27.7503372936239*I + sage: a.call(a.polylog, 2, 1/3+4/5*I) + 0.153548951541433 + 0.875114412499637*I + sage: a.call(a.barnesg, 3+4*I) + -0.000676375932234244 - 0.0000442236140124728*I + sage: a.call(a.barnesg, -4) + 0.000000000000000 + sage: a.call(a.hyper, [2,3], [4,5], 1/3) + 1.10703578162508 + sage: a.call(a.hyper, [2,3], [4,(2,3)], 1/3) + 1.95762943509305 + sage: a.call(a.quad, a.erf, [0,1]) + 0.486064958112256 + sage: a.call(a.gammainc, 3+4*I, 2/3, 1-pi*I, prec=100) + -274.18871130777160922270612331 + 101.59521032382593402947725236*I + sage: x = (3+4*I).n(100) + sage: y = (2/3).n(100) + sage: z = (1-pi*I).n(100) + sage: a.call(a.gammainc, x, y, z, prec=100) + -274.18871130777160922270612331 + 101.59521032382593402947725236*I + sage: a.call(a.erf, infinity) + 1.00000000000000 + sage: a.call(a.erf, -infinity) + -1.00000000000000 + sage: a.call(a.gamma, infinity) + +infinity + sage: a.call(a.polylog, 2, 1/2, parent=RR) + 0.582240526465012 + sage: a.call(a.polylog, 2, 2, parent=RR) + 2.46740110027234 - 2.17758609030360*I + sage: a.call(a.polylog, 2, 1/2, parent=RealField(100)) + 0.58224052646501250590265632016 + sage: a.call(a.polylog, 2, 2, parent=RealField(100)) + 2.4674011002723396547086227500 - 2.1775860903036021305006888982*I + sage: a.call(a.polylog, 2, 1/2, parent=CC) + 0.582240526465012 + sage: type(_) + + sage: a.call(a.polylog, 2, 1/2, parent=RDF) + 0.5822405264650125 + sage: type(_) + + + Check that :issue:`11885` is fixed:: + + sage: a.call(a.ei, 1.0r, parent=float) + 1.8951178163559366 + + Check that :issue:`14984` is fixed:: + + sage: a.call(a.log, -1.0r, parent=float) + 3.141592653589793j + """ + from mpmath import mp + orig = mp.prec + prec = kwargs.pop('prec', orig) + parent = kwargs.pop('parent', None) + if parent is not None: + try: + prec = parent.prec() + except AttributeError: + pass + prec2 = prec + 20 + args = sage_to_mpmath(args, prec2) + kwargs = sage_to_mpmath(kwargs, prec2) + try: + mp.prec = prec + y = func(*args, **kwargs) + finally: + mp.prec = orig + y = mpmath_to_sage(y, prec) + if parent is None: + return y + try: + return parent(y) + except TypeError as error: + try: + return parent.complex_field()(y) + except AttributeError: + if parent is float: + return complex(y) + else: + raise TypeError(error) diff --git a/src/sage/rings/real_mpfr.pyx b/src/sage/rings/real_mpfr.pyx index 10ff6b8baae..931ece0dd92 100644 --- a/src/sage/rings/real_mpfr.pyx +++ b/src/sage/rings/real_mpfr.pyx @@ -150,8 +150,6 @@ from sage.structure.richcmp cimport rich_to_bool_sgn cdef bin_op from sage.structure.element import bin_op -from sage.libs.mpmath.utils cimport mpfr_to_mpfval - from sage.rings.integer cimport Integer from sage.rings.rational cimport Rational from sage.rings.real_double cimport RealDoubleElement diff --git a/src/sage/structure/coerce.pyx b/src/sage/structure/coerce.pyx index b00ac6487a0..541e0600fe1 100644 --- a/src/sage/structure/coerce.pyx +++ b/src/sage/structure/coerce.pyx @@ -100,6 +100,95 @@ import traceback from fractions import Fraction cdef type FractionType = Fraction +cpdef _py_scalar_parent_mpmath13(py_type): + r""" + Version of py_scalar_parent() suitable for use with mpmath-1.3.x. + The main py_scalar_parent() will invoke this one if needed. + + EXAMPLES:: + + sage: # needs mpmath13 + sage: from sage.structure.coerce import _py_scalar_parent_mpmath13 + sage: _py_scalar_parent_mpmath13(int) + Integer Ring + sage: _py_scalar_parent_mpmath13(float) + Real Double Field + sage: _py_scalar_parent_mpmath13(complex) # needs sage.rings.complex_double + Complex Double Field + sage: _py_scalar_parent_mpmath13(bool) + Integer Ring + sage: _py_scalar_parent_mpmath13(dict), + (None,) + + sage: # needs mpmath13 + sage: import fractions + sage: _py_scalar_parent_mpmath13(fractions.Fraction) + Rational Field + + sage: # needs numpy mpmath13 + sage: import numpy + sage: _py_scalar_parent_mpmath13(numpy.int16) + Integer Ring + sage: _py_scalar_parent_mpmath13(numpy.int32) + Integer Ring + sage: _py_scalar_parent_mpmath13(numpy.uint64) + Integer Ring + sage: _py_scalar_parent_mpmath13(numpy.double) + Real Double Field + + sage: # needs mpmath13 + sage: import gmpy2 + sage: _py_scalar_parent_mpmath13(gmpy2.mpz) + Integer Ring + sage: _py_scalar_parent_mpmath13(gmpy2.mpq) + Rational Field + sage: _py_scalar_parent_mpmath13(gmpy2.mpfr) + Real Double Field + sage: _py_scalar_parent_mpmath13(gmpy2.mpc) # needs sage.rings.complex_double + Complex Double Field + + """ + if issubclass(py_type, int): + import sage.rings.integer_ring + return sage.rings.integer_ring.ZZ + if py_type is FractionType: + import sage.rings.rational_field + return sage.rings.rational_field.QQ + elif issubclass(py_type, float): + import sage.rings.real_double + return sage.rings.real_double.RDF + elif issubclass(py_type, complex): + import sage.rings.complex_double + return sage.rings.complex_double.CDF + elif is_numpy_type(py_type): + import numpy + if issubclass(py_type, numpy.integer): + import sage.rings.integer_ring + return sage.rings.integer_ring.ZZ + elif issubclass(py_type, numpy.floating): + import sage.rings.real_double + return sage.rings.real_double.RDF + elif issubclass(py_type, numpy.complexfloating): + import sage.rings.complex_double + return sage.rings.complex_double.CDF + else: + return None + elif issubclass(py_type, gmpy2.mpz): + import sage.rings.integer_ring + return sage.rings.integer_ring.ZZ + elif issubclass(py_type, gmpy2.mpq): + import sage.rings.rational_field + return sage.rings.rational_field.QQ + elif issubclass(py_type, gmpy2.mpfr): + import sage.rings.real_double + return sage.rings.real_double.RDF + elif issubclass(py_type, gmpy2.mpc): + import sage.rings.complex_double + return sage.rings.complex_double.CDF + else: + return None + + cpdef py_scalar_parent(py_type): """ Return the Sage equivalent of the given python type, if one exists. @@ -143,46 +232,66 @@ cpdef py_scalar_parent(py_type): Real Double Field sage: py_scalar_parent(gmpy2.mpc) # needs sage.rings.complex_double Complex Double Field + + sage: # needs mpmath + sage: import mpmath + sage: py_scalar_parent(mpmath.mpf) + Real Double Field + sage: py_scalar_parent(mpmath.mpc) # needs sage.rings.complex_double + Complex Double Field """ + import mpmath + if not mpmath.__version__.startswith("1.4"): + # Support mpmath-1.3 for a while by copy/pasting the old + # version of the function. + return _py_scalar_parent_mpmath13(py_type) + if issubclass(py_type, int): import sage.rings.integer_ring return sage.rings.integer_ring.ZZ if py_type is FractionType: import sage.rings.rational_field return sage.rings.rational_field.QQ - elif issubclass(py_type, float): + if issubclass(py_type, float): import sage.rings.real_double return sage.rings.real_double.RDF - elif issubclass(py_type, complex): + if issubclass(py_type, complex): import sage.rings.complex_double return sage.rings.complex_double.CDF - elif is_numpy_type(py_type): + if is_numpy_type(py_type): import numpy if issubclass(py_type, numpy.integer): import sage.rings.integer_ring return sage.rings.integer_ring.ZZ - elif issubclass(py_type, numpy.floating): + if issubclass(py_type, numpy.floating): import sage.rings.real_double return sage.rings.real_double.RDF - elif issubclass(py_type, numpy.complexfloating): + if issubclass(py_type, numpy.complexfloating): import sage.rings.complex_double return sage.rings.complex_double.CDF - else: - return None - elif issubclass(py_type, gmpy2.mpz): + return None + if issubclass(py_type, gmpy2.mpz): import sage.rings.integer_ring return sage.rings.integer_ring.ZZ - elif issubclass(py_type, gmpy2.mpq): + if issubclass(py_type, gmpy2.mpq): import sage.rings.rational_field return sage.rings.rational_field.QQ - elif issubclass(py_type, gmpy2.mpfr): + if issubclass(py_type, gmpy2.mpfr): import sage.rings.real_double return sage.rings.real_double.RDF - elif issubclass(py_type, gmpy2.mpc): + if issubclass(py_type, gmpy2.mpc): import sage.rings.complex_double return sage.rings.complex_double.CDF - else: + if is_mpmath_type(py_type): + import mpmath + if issubclass(py_type, mpmath.mpf): + from sage.rings.real_double import RDF + return RDF + if issubclass(py_type, mpmath.mpc): + from sage.rings.complex_double import CDF + return CDF return None + return None cpdef py_scalar_to_element(x): """ @@ -468,10 +577,11 @@ cpdef bint is_numpy_type(t) noexcept: return True return False + cpdef bint is_mpmath_type(t) noexcept: r""" - Check whether the type ``t`` is a type whose name starts with either - ``mpmath.`` or ``sage.libs.mpmath.``. + Check whether the type ``t`` is a type whose name starts with + "mpmath." or "sage.libs.mpmath" EXAMPLES:: @@ -487,8 +597,11 @@ cpdef bint is_mpmath_type(t) noexcept: sage: is_mpmath_type(type(mpmath.mpf(2))) True """ - return isinstance(t, type) and \ - strncmp((t).tp_name, "sage.libs.mpmath.", 17) == 0 + return isinstance(t, type) and ( + t.__module__.startswith("mpmath.") + or + t.__module__.startswith("sage.libs.mpmath") + ) cdef class CoercionModel: diff --git a/src/sage/tests/books/computational-mathematics-with-sagemath/integration_doctest.py b/src/sage/tests/books/computational-mathematics-with-sagemath/integration_doctest.py index bf2701868ad..0603d419848 100644 --- a/src/sage/tests/books/computational-mathematics-with-sagemath/integration_doctest.py +++ b/src/sage/tests/books/computational-mathematics-with-sagemath/integration_doctest.py @@ -151,7 +151,7 @@ sage: mpmath.quad(f, [0, 1]) Traceback (most recent call last): ... - TypeError: no canonical coercion from to ... + TypeError: no canonical coercion from to ... Sage example in ./integration.tex, line 866::