Repository: mpmath/mpmath

Issue:
Interval arithmetic: incorrect results for min and max
```
Consider the following operation:

a = iv.mpf([0,1])
b = 0.5

min(a,b)

Mpmath returns the range [0,1] though there is no posible value of a that can give a result greater than 0.5

The correct response is [0,0.5]
```

Original issue for #258: http://code.google.com/p/mpmath/issues/detail?id=218
Original author: https://code.google.com/u/117464931809270512527/

Hints from issue comments:
@fredrik-johansson: ```
This behavior is impossible to change since min is a Python function implemented in terms of comparison operators.

We could perhaps add new interval min/max functions to mpmath that behave this way?
```

Original comment: http://code.google.com/p/mpmath/issues/detail?id=218#c1
Original author: https://code.google.com/u/111502149103757882156/
@warnerjon12: Solution would be to implement comparison operators in ivmpf that could work with other scalar types@skirpichev: @warnerjon12, the problem is that rich-comparison helpers for iv sometimes return ``None``.  In a context, where boolean value is required (e.g. in min/max) - Python calls ``bool(None)``, that's ``False``.@skirpichev: One option might be raise a ``ValueError`` instead of return ``None``.  This will prevent silent wrong answer from min/max builtins (and perhaps, from more, like ``sorted``).

@fredrik-johansson [mentioned](https://github.com/mpmath/mpmath/issues/191#issuecomment-41250056), that library code might be adapted for triple-valued logic.  But I'm not sure that it's a good idea to try handle unordered "numbers" in generic code.  Better - fail quickly.

BTW, python-flint just returns False:
```pycon
>>> flint.arb('0.0 +/- 0.1') > flint.arb("0.0 +/- 0.1")
False
>>> flint.arb('0.0 +/- 0.1') < flint.arb("0.0 +/- 0.1")
False
>>> flint.arb('0.0 +/- 0.1') <= flint.arb("0.0 +/- 0.1")
False
>>> flint.arb('0.0 +/- 0.1') >= flint.arb("0.0 +/- 0.1")
False
```@oscarbenjamin: > BTW, python-flint just returns False:

We might want to change that. Flint's generics use three-valued logic and I expect that we would want to expose that somehow although I'm not sure if that might be by adding something like a new `T_UNKNOWN` object rather than using `None`. The idea would be that `T_UNKNOWN.__bool__` would raise an exception if so.@skirpichev: > Flint's generics use three-valued logic

Hmm, how this works?  I'm not sure what helpful could be done if both ``bool(x > y)`` and ``bool(x <= y)`` are ``False``.

> The idea would be that T_UNKNOWN.__bool__ would raise an exception if so.

That looks more like my proposal.@oscarbenjamin: > Hmm, how this works?

It is explained here:
https://flintlib.org/doc/gr.html#predicates

I'm just not sure how to translate that into Python yet. In the PR https://github.com/flintlib/python-flint/pull/218 I have mapped `T_UNKNOWN` to `None` so far although ordered comparisons with `gr_` types are not implemented there yet:
```python
In [11]: from flint.types import _gr

In [12]: print(ctx('0.0 +/- 0.1') == ctx('0.0 +/- 0.1'))
None

In [13]: print(ctx('0.0 +/- 0.1') != ctx('0.0 +/- 0.1'))
None

In [20]: print(ctx('0.0 +/- 0.1') == ctx('2.0 +/- 0.1'))
False

In [21]: print(ctx('0.0 +/- 0.1') != ctx('2.0 +/- 0.1'))
True
```
That is just the logic implemented by the generic rings code in FLINT except with `T_UNKNOWN` being returned as `None` in Python.

Also worth noting that the PR adds versions of `arf` and `acf` that use non-global context for precision which is what mpmath could use as faster versions of `mpf` and `mpc` internally:
```python
In [15]: ctx = _gr.gr_real_float_arf_ctx.new(53)

In [16]: ctx(1) / 3
Out[16]: 0.3333333333333333

In [17]: ctx.real_prec = 100

In [18]: ctx(1) / 3
Out[18]: 0.3333333333333333333333333333335

In [10]: from mpmath import mp as mpmath_ctx

In [11]: flint_ctx = _gr.gr_real_float_arf_ctx.new(53)

In [12]: a_flint = flint_ctx(20)

In [13]: a_mpmath = mpmath_ctx.mpf(20)

In [14]: %timeit 2*a_mpmath
4.2 μs ± 14.9 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [15]: %timeit 2*a_flint
684 ns ± 19.7 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

In [16]: %timeit 1/a_mpmath
6.78 μs ± 243 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [17]: %timeit 1/a_flint
1.13 μs ± 6.93 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
```
I expect that any Python wrapper class around `arf/acf` would end up having practically as much overhead as `mpf/mpc` already do in comparison but if some internal code could be made to be able to work with one or the other transparently then it could bring significant speedup.

The `_gr` interface should be considered unstable but could be suitable for testing python-flint/mpmath integration if merged and put in a prerelease.@skirpichev: > It is explained here: https://flintlib.org/doc/gr.html#predicates

I meant, how it's used in some generic code (e.g. to implement special functions)?  I would expect that numeric code on some T_UNKNOWN condition - should just fail.@oscarbenjamin: > I meant, how it's used in some generic code (e.g. to implement special functions)?

It isn't used in that way in arf/arb. I don't think they use `T_UNKNOWN` directly at all so none of the internal code for computing special functions uses this. Just when arb is wrapped in the `gr_` interface comparison operators like `==`, `<` etc can return `T_UNKNOWN` and potentially it makes sense to expose that at the Python level for these types. I think the internal code for arb would do something like `if midpoint_as_arf(some_arb) < some_value` which would always return True or False e.g.:
```python
In [2]: a = arb('1 +/- 0.1')

In [3]: a
Out[3]: [1e+0 +/- 0.101]

In [4]: a.mid()
Out[4]: 1.00000000000000

In [5]: a.mid() < 1
Out[5]: False

In [6]: a.mid() == 1
Out[6]: True

In [7]: type(a.mid())
Out[7]: flint.types.arb.arb
```@skirpichev: > think the internal code for arb would do something like if midpoint_as_arf(some_arb) < some_value

In this way we loose most benefits of the interval arithmetic, isn't?@oscarbenjamin: Oh, actually it is just that the arb functions don't use three values for logical predicates so you have `arb_le(a, b)` but it returns True if it is definitely true but False otherwise i.e. it does not distinguish between False and None: you would have to call `arb_gt(a, b)` to distinguish those. It looks like currently python-flint just exposes this sort of logic directly for the `arb` type whereas the `gr` wrapper of arb uses `truth_t` with `T_UNKNOWN` instead.

Probably in python-flint it makes more sense to use three-valued logic somehow in a way that is consistent across all types.@skirpichev: Well, I think that raising an exception seems to be a simplest option than using three-value logic (that will require special implementation of all functions for the iv context).

Proposed patch breaks only few test.

Solution patch:
diff --git a/docs/contexts.rst b/docs/contexts.rst
index 77ab898..606def7 100644
--- a/docs/contexts.rst
+++ b/docs/contexts.rst
@@ -144,7 +144,11 @@ Intervals may be infinite or half-infinite::
     >>> print(1 / iv.mpf([2, 'inf']))
     [0.0, 0.5]
 
-The equality testing operators ``==`` and ``!=`` check whether their operands are identical as intervals; that is, have the same endpoints. The ordering operators ``< <= > >=`` permit inequality testing using triple-valued logic: a guaranteed inequality returns ``True`` or ``False`` while an indeterminate inequality returns ``None``::
+The equality testing operators ``==`` and ``!=`` check whether their operands
+are identical as intervals; that is, have the same endpoints. The ordering
+operators ``< <= > >=`` permit inequality testing using triple-valued logic: a
+guaranteed inequality returns ``True`` or ``False`` while an indeterminate
+inequality raises :exc:`ValueError`::
 
     >>> iv.mpf([1,2]) == iv.mpf([1,2])
     True
@@ -156,12 +160,18 @@ The equality testing operators ``==`` and ``!=`` check whether their operands ar
     True
     >>> iv.mpf([1,2]) < 1
     False
-    >>> iv.mpf([1,2]) < 2    # returns None
+    >>> iv.mpf([1,2]) < 2
+    Traceback (most recent call last):
+      ...
+    ValueError
     >>> iv.mpf([2,2]) < 2
     False
     >>> iv.mpf([1,2]) <= iv.mpf([2,3])
     True
-    >>> iv.mpf([1,2]) < iv.mpf([2,3])  # returns None
+    >>> iv.mpf([1,2]) < iv.mpf([2,3])
+    Traceback (most recent call last):
+      ...
+    ValueError
     >>> iv.mpf([1,2]) < iv.mpf([-1,0])
     False
 
@@ -236,13 +246,19 @@ seen by increasing the precision::
     >>> print(mp.exp(mp.pi*mp.sqrt(163)))
     262537412640768743.99999999999925007259719818568888
 
-With interval arithmetic, the comparison returns ``None`` until the precision
-is large enough for `x-y` to have a definite sign::
+With interval arithmetic, the comparison raises :exc:`ValueError` until the
+precision is large enough for `x-y` to have a definite sign::
 
     >>> iv.dps = 15
     >>> iv.exp(iv.pi*iv.sqrt(163)) > (640320**3+744)
+    Traceback (most recent call last):
+      ...
+    ValueError
     >>> iv.dps = 30
     >>> iv.exp(iv.pi*iv.sqrt(163)) > (640320**3+744)
+    Traceback (most recent call last):
+      ...
+    ValueError
     >>> iv.dps = 60
     >>> iv.exp(iv.pi*iv.sqrt(163)) > (640320**3+744)
     False
diff --git a/mpmath/ctx_fp.py b/mpmath/ctx_fp.py
index 354fe83..e6bd5fa 100644
--- a/mpmath/ctx_fp.py
+++ b/mpmath/ctx_fp.py
@@ -202,16 +202,40 @@ def _is_real_type(ctx, z):
     def _is_complex_type(ctx, z):
         return isinstance(z, complex)
 
-    def hypsum(ctx, p, q, types, coeffs, z, maxterms=6000, **kwargs):
-        coeffs = list(coeffs)
+    def hypsum(ctx, p, q, flags, coeffs, z, maxterms=6000, **kwargs):
+        for i, c in enumerate(coeffs[p:], start=p):
+            if flags[i] == 'Z':
+                if c <= 0:
+                    ok = False
+                    for ii, cc in enumerate(coeffs[:p]):
+                        # Note: c <= cc or c < cc, depending on convention
+                        if flags[ii] == 'Z' and cc <= 0 and c <= cc:
+                            ok = True
+                    if not ok:
+                        raise ZeroDivisionError("pole in hypergeometric series")
         num = range(p)
         den = range(p,p+q)
+        if ctx.isinf(z):
+            n = max(((n, c) for n, c in enumerate(coeffs[:p])
+                     if flags[n] == 'Z' and c < 0), default=(-1, 0),
+                    key=lambda x: x[1])[0]
+            if n >= 0:
+                n = -coeffs[n]
+                t = z**n
+                for k in range(n):
+                    for i in num: t *= (coeffs[i]+k)
+                    for i in den: t /= (coeffs[i]+k)
+                    t /= (k+1)
+                return t
         tol = ctx.eps
         s = t = 1.0
         k = 0
         while 1:
             for i in num: t *= (coeffs[i]+k)
-            for i in den: t /= (coeffs[i]+k)
+            try:
+                for i in den: t /= (coeffs[i]+k)
+            except ZeroDivisionError:
+                raise NotImplementedError
             k += 1; t /= k; t *= z; s += t
             if abs(t) < tol:
                 return s
diff --git a/mpmath/ctx_mp.py b/mpmath/ctx_mp.py
index d815dde..d2fde03 100644
--- a/mpmath/ctx_mp.py
+++ b/mpmath/ctx_mp.py
@@ -645,6 +645,30 @@ def hypsum(ctx, p, q, flags, coeffs, z, accurate_small=True, **kwargs):
         elif hasattr(z, "_mpc_"):
             key = p, q, flags, 'C'
             v = z._mpc_
+        for i, c in enumerate(coeffs[p:], start=p):
+            if flags[i] == 'Z':
+                if c <= 0:
+                    ok = False
+                    for ii, cc in enumerate(coeffs[:p]):
+                        # Note: c <= cc or c < cc, depending on convention
+                        if flags[ii] == 'Z' and cc <= 0 and c <= cc:
+                            ok = True
+                    if not ok:
+                        raise ZeroDivisionError("pole in hypergeometric series")
+        num = range(p)
+        den = range(p,p+q)
+        if ctx.isinf(z):
+            n = max(((n, c) for n, c in enumerate(coeffs[:p])
+                     if flags[n] == 'Z' and c < 0), default=(-1, 0),
+                    key=lambda x: x[1])[0]
+            if n >= 0:
+                n = -coeffs[n]
+                t = z**n
+                for k in range(n):
+                    for i in num: t *= (coeffs[i]+k)
+                    for i in den: t /= (coeffs[i]+k)
+                    t /= (k+1)
+                return t
         if key not in ctx.hyp_summators:
             ctx.hyp_summators[key] = libmp.make_hyp_summator(key)[1]
         summator = ctx.hyp_summators[key]
@@ -659,14 +683,6 @@ def hypsum(ctx, p, q, flags, coeffs, z, accurate_small=True, **kwargs):
         max_total_jump = 0
         for i, c in enumerate(coeffs):
             if flags[i] == 'Z':
-                if i >= p and c <= 0:
-                    ok = False
-                    for ii, cc in enumerate(coeffs[:p]):
-                        # Note: c <= cc or c < cc, depending on convention
-                        if flags[ii] == 'Z' and cc <= 0 and c <= cc:
-                            ok = True
-                    if not ok:
-                        raise ZeroDivisionError("pole in hypergeometric series")
                 continue
             n, d = ctx.nint_distance(c)
             n = -int(n)
diff --git a/mpmath/libmp/libmpi.py b/mpmath/libmp/libmpi.py
index e3d8e70..8aec5f5 100644
--- a/mpmath/libmp/libmpi.py
+++ b/mpmath/libmp/libmpi.py
@@ -37,14 +37,14 @@ def mpi_lt(s, t):
     ta, tb = t
     if mpf_lt(sb, ta): return True
     if mpf_ge(sa, tb): return False
-    return None
+    raise ValueError
 
 def mpi_le(s, t):
     sa, sb = s
     ta, tb = t
     if mpf_le(sb, ta): return True
     if mpf_gt(sa, tb): return False
-    return None
+    raise ValueError
 
 def mpi_gt(s, t): return mpi_lt(t, s)
 def mpi_ge(s, t): return mpi_le(t, s)


Test patch:
diff --git a/mpmath/tests/test_functions2.py b/mpmath/tests/test_functions2.py
index c57ec95..1176fb4 100644
--- a/mpmath/tests/test_functions2.py
+++ b/mpmath/tests/test_functions2.py
@@ -165,6 +165,8 @@ def test_hyper_misc():
     mp.dps = 25
     v = mpc('1.2282306665029814734863026', '-0.1225033830118305184672133')
     assert hyper([(3,4),2+j,1],[1,5,j/3],mpf(1)/5+j/8).ae(v)
+    pytest.raises(ZeroDivisionError, lambda: mp.hyper([1, 2, -2], [-1, 3], 1.1))
+    pytest.raises(ZeroDivisionError, lambda: fp.hyper([1, 2, -2], [-1, 3], 1.1))
 
 def test_elliptic_integrals():
     assert ellipk(0).ae(pi/2)
@@ -528,6 +530,7 @@ def test_hyper_2f1():
 def test_hyper_2f1_hard():
     # Singular cases
     assert hyp2f1(2,-1,-1,3).ae(7)
+    pytest.raises(NotImplementedError, lambda: fp.hyp2f1(2,-1,-1,3))
     assert hyp2f1(2,-1,-1,3,eliminate_all=True).ae(0.25)
     assert hyp2f1(2,-2,-2,3).ae(34)
     assert hyp2f1(2,-2,-2,3,eliminate_all=True).ae(0.25)
@@ -2412,3 +2415,8 @@ def test_issue_505():
 
 def test_issue_653():
     pytest.raises(ValueError, lambda: zeta(2, -2))
+
+def test_issue_511():
+    assert mp.laguerre(1, 2, mp.inf) == -mp.inf
+    assert mp.laguerre(1, 7.2, mp.inf) == -mp.inf
+    assert fp.laguerre(1, 7.2, fp.inf) == -fp.inf
diff --git a/mpmath/tests/test_interval.py b/mpmath/tests/test_interval.py
index 93e7f8b..f39a270 100644
--- a/mpmath/tests/test_interval.py
+++ b/mpmath/tests/test_interval.py
@@ -1,4 +1,6 @@
-from mpmath import inf, iv, mp, mpf, mpi, pi, sqrt
+import pytest
+
+from mpmath import inf, iv, mp, mpf, mpi, pi, sqrt, workprec
 
 
 def test_interval_identity():
@@ -303,8 +305,9 @@ def test_interval_complex():
     assert iv.mpc(2,2) ** (-2) == (2+2j) ** (-2)
     assert iv.cos(2).ae(mp.cos(2))
     assert iv.sin(2).ae(mp.sin(2))
-    assert iv.cos(2+3j).ae(mp.cos(2+3j))
-    assert iv.sin(2+3j).ae(mp.sin(2+3j))
+    with workprec(54):
+        assert iv.cos(2+3j).ae(mp.cos(2+3j))
+        assert iv.sin(2+3j).ae(mp.sin(2+3j))
 
 def test_interval_complex_arg():
     assert iv.arg(3) == 0
@@ -351,9 +354,9 @@ def test_interval_complex_arg():
 
 def test_interval_ae():
     x = iv.mpf([1,2])
-    assert x.ae(1) is None
-    assert x.ae(1.5) is None
-    assert x.ae(2) is None
+    pytest.raises(ValueError, lambda: x.ae(1))
+    pytest.raises(ValueError, lambda: x.ae(1.5))
+    pytest.raises(ValueError, lambda: x.ae(2))
     assert x.ae(2.01) is False
     assert x.ae(0.99) is False
     x = iv.mpf(3.5)
@@ -362,8 +365,8 @@ def test_interval_ae():
     assert x.ae(3.5-1e-15) is True
     assert x.ae(3.501) is False
     assert x.ae(3.499) is False
-    assert x.ae(iv.mpf([3.5,3.501])) is None
-    assert x.ae(iv.mpf([3.5,4.5+1e-15])) is None
+    pytest.raises(ValueError, lambda: x.ae(iv.mpf([3.5,3.501])))
+    pytest.raises(ValueError, lambda: x.ae(iv.mpf([3.5,4.5+1e-15])))
 
 def test_interval_nstr():
     iv.dps = n = 30
@@ -439,3 +442,9 @@ def test_interval_conversions():
         assert float(r.b) == float(b)
         assert complex(r.a) == complex(a)
         assert complex(r.b) == complex(b)
+
+def test_issue_258():
+    a = iv.mpf([0, 1])
+    b = 0.5
+    pytest.raises(ValueError, lambda: min(a, b))
+    pytest.raises(ValueError, lambda: max(a, b))