Repository: numpy/numpy
Pull request: BUG: Loss of precision in longdouble min
Generic reuse in the latest changes around `min` works unless the macro is redefined for SIMD.

This change avoids `scalar_min_f` for generic comparisons (as it can be redefined) and defines is separately as `scalar_min`

This should resolve https://github.com/numpy/numpy/issues/20863

Problem statement:
BUG: loss of precision in longdouble min (nightly build)
### Describe the issue:

My pre-release tests have begun breaking, and I've tracked it down to a small issue. See the minimal reproduction below.

The use case is rescaling an array into a smaller dtype, which requires finding the bounds of the array and target dtype. The specific failing use case is when squeezing float128 into int32, but I provide below a more thorough regression test that demonstrates that this occurs in a number of cases, including float16 -> int16, float32 -> int32, and float128 -> int64.

### Reproduce the code example:

```python
import numpy as np
np.min(np.float128([np.inf, 2147483647]))
# Result: 2147483648.0

# More thorough test
for ft in np.sctypes["float"]:
    for it in np.sctypes["int"]:
        fmax = np.finfo(ft).max
        imax = np.iinfo(it).max
        shared_max = np.min(ft([fmax, imax]))
        check = shared_max <= imax
        print(f"{it} {ft} {check} {fmax} {imax} {shared_max}")
```


### Error message:

```shell
<class 'numpy.int8'> <class 'numpy.float16'> True 65504.0 127 127.0
<class 'numpy.int16'> <class 'numpy.float16'> False 65504.0 32767 32768.0
<class 'numpy.int32'> <class 'numpy.float16'> True 65504.0 2147483647 65504.0
<class 'numpy.int64'> <class 'numpy.float16'> True 65504.0 9223372036854775807 65504.0
<class 'numpy.int8'> <class 'numpy.float32'> True 3.4028234663852886e+38 127 127.0
<class 'numpy.int16'> <class 'numpy.float32'> True 3.4028234663852886e+38 32767 32767.0
<class 'numpy.int32'> <class 'numpy.float32'> False 3.4028234663852886e+38 2147483647 2147483648.0
<class 'numpy.int64'> <class 'numpy.float32'> True 3.4028234663852886e+38 9223372036854775807 9.223372036854776e+18
<class 'numpy.int8'> <class 'numpy.float64'> True 1.7976931348623157e+308 127 127.0
<class 'numpy.int16'> <class 'numpy.float64'> True 1.7976931348623157e+308 32767 32767.0
<class 'numpy.int32'> <class 'numpy.float64'> True 1.7976931348623157e+308 2147483647 2147483647.0
<class 'numpy.int64'> <class 'numpy.float64'> True 1.7976931348623157e+308 9223372036854775807 9.223372036854776e+18
<class 'numpy.int8'> <class 'numpy.float128'> True inf 127 127.0
<class 'numpy.int16'> <class 'numpy.float128'> True inf 32767 32767.0
<class 'numpy.int32'> <class 'numpy.float128'> False inf 2147483647 2147483648.0
<class 'numpy.int64'> <class 'numpy.float128'> False inf 9223372036854775807 9.223372036854776e+18
```

### NumPy/Python version information:

1.23.0.dev0+437.g4b985e2b8 3.8.12 (default, Oct 12 2021, 13:49:34) 
[GCC 7.5.0]

Hints from issue discussion:
Thanks!  Linked/mentioned in the above PR which is likely the cause.  It sounds a bit like a temporary or so may just be incorrectly typed, because your magic "off by one" is gotten also if you do the following:
```
np.float64(np.float32(2147483647))
```
I am not sure the table is very helpful.  The problem is that `np.float32(2147483647)` will lose precision, and that is completely fine (the float needs to have more bits than the integer, or precision is lost earlier).

However:
```
np.min(np.float128([np.inf, 2147483647]))
```
should not lose any precision.

However, I can only preproduce this for float128, so it seems like the issue is limited to `longdouble` only.
> I am not sure the table is very helpful.

<details>
<summary>You're right. Only the last two change when reverting to stable. (Click to expand table.)</summary>

```
<class 'numpy.int8'> <class 'numpy.float16'> True 65504.0 127 127.0
<class 'numpy.int16'> <class 'numpy.float16'> False 65504.0 32767 32768.0
<class 'numpy.int32'> <class 'numpy.float16'> True 65504.0 2147483647 65504.0
<class 'numpy.int64'> <class 'numpy.float16'> True 65504.0 9223372036854775807 65504.0
<class 'numpy.int8'> <class 'numpy.float32'> True 3.4028234663852886e+38 127 127.0
<class 'numpy.int16'> <class 'numpy.float32'> True 3.4028234663852886e+38 32767 32767.0
<class 'numpy.int32'> <class 'numpy.float32'> False 3.4028234663852886e+38 2147483647 2147483648.0
<class 'numpy.int64'> <class 'numpy.float32'> True 3.4028234663852886e+38 9223372036854775807 9.223372036854776e+18
<class 'numpy.int8'> <class 'numpy.float64'> True 1.7976931348623157e+308 127 127.0
<class 'numpy.int16'> <class 'numpy.float64'> True 1.7976931348623157e+308 32767 32767.0
<class 'numpy.int32'> <class 'numpy.float64'> True 1.7976931348623157e+308 2147483647 2147483647.0
<class 'numpy.int64'> <class 'numpy.float64'> True 1.7976931348623157e+308 9223372036854775807 9.223372036854776e+18
<class 'numpy.int8'> <class 'numpy.float128'> True inf 127 127.0
<class 'numpy.int16'> <class 'numpy.float128'> True inf 32767 32767.0
<class 'numpy.int32'> <class 'numpy.float128'> True inf 2147483647 2147483647.0
<class 'numpy.int64'> <class 'numpy.float128'> True inf 9223372036854775807 9.223372036854776e+18
```

</details>

I suppose the check should be:

```Python
for ft in np.sctypes["float"]:
    for it in np.sctypes["int"]:
        fmax = np.finfo(ft).max
        imax = np.iinfo(it).max
        cast_imax = ft(imax)
        shared_max = np.min([fmax, cast_imax])
        check = shared_max <= cast_imax
        print(f"{it} {ft} {check} {fmax} {imax}({cast_imax}) {shared_max}")
```
https://github.com/numpy/numpy/pull/20872 intended to resolve this issue
Thinking more about trying to capture the precision loss problem suggests a refined test case. Using `fmax = np.info(dt).max`, any loss of precision must produce a different value. For reducing functions, which will (almost?) always have identity values, we could write a check:

```Python
def check(dt, func, identity):
    fmax = np.finfo(dt).max
    assert func([fmax, identity]) == fmax
```

Some obvious results:

```Python
check(np.float128, np.min, np.inf)
check(np.float128, np.max, -np.inf)
check(np.float128, np.sum, 0.)
check(np.float128, np.prod, 1.)
```

Not sure if there are other broad classes that could be checked as easily.
I would suggest using something like `np.finfo(...).max`, a subnormal (any will do), and maybe an additional "typical" but full precision value, e.g. `eps`, `0.1`...).
EDIT: OK, I guess `0.1` doesn't quite do the trick :), `np.nextafter` could be used also for simplicit, though.

Commit message: Add tests as requested

Adding tests like suggested in https://github.com/numpy/numpy/issues/20863

Patch to review:
diff --git a/numpy/core/src/umath/loops_minmax.dispatch.c.src b/numpy/core/src/umath/loops_minmax.dispatch.c.src
index 708270016d63..ba2288f0b9dc 100644
--- a/numpy/core/src/umath/loops_minmax.dispatch.c.src
+++ b/numpy/core/src/umath/loops_minmax.dispatch.c.src
@@ -22,12 +22,14 @@
 #define scalar_max_i(A, B) ((A > B) ? A : B)
 #define scalar_min_i(A, B) ((A < B) ? A : B)
 // fp, propagates NaNs
-#define scalar_max_f(A, B) ((A >= B || npy_isnan(A)) ? A : B)
-#define scalar_max_d scalar_max_f
-#define scalar_max_l scalar_max_f
-#define scalar_min_f(A, B) ((A <= B || npy_isnan(A)) ? A : B)
-#define scalar_min_d scalar_min_f
-#define scalar_min_l scalar_min_f
+#define scalar_max(A, B) ((A >= B || npy_isnan(A)) ? A : B)
+#define scalar_max_f scalar_max
+#define scalar_max_d scalar_max
+#define scalar_max_l scalar_max
+#define scalar_min(A, B) ((A <= B || npy_isnan(A)) ? A : B)
+#define scalar_min_f scalar_min
+#define scalar_min_d scalar_min
+#define scalar_min_l scalar_min
 // fp, ignores NaNs
 #define scalar_maxp_f fmaxf
 #define scalar_maxp_d fmax
diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py
index e7fee46b7b79..734a2fcf85e3 100644
--- a/numpy/core/tests/test_umath.py
+++ b/numpy/core/tests/test_umath.py
@@ -1728,6 +1728,29 @@ def test_strided_array(self):
         assert_equal(np.maximum(arr1[:6:2], arr2[::3], out=out[::3]), np.array([-2.0, 10., np.nan]))
         assert_equal(out, out_maxtrue)
 
+    def test_precision(self):
+        dtypes = [np.float16, np.float32, np.float64]
+        if(hasattr(np, 'float128')):
+            dtypes.append(np.float128)
+
+        for dt in dtypes:
+            dtmin = np.finfo(dt).min
+            dtmax = np.finfo(dt).max
+            d1 = dt(0.1)
+            d1_next = np.nextafter(d1, np.inf)
+
+            test_cases = [
+                # v1    v2          expected
+                (dtmin, -np.inf,    dtmin),
+                (dtmax, -np.inf,    dtmax),
+                (d1,    d1_next,    d1_next),
+                (dtmax, np.nan,     np.nan),
+            ]
+
+            for v1, v2, expected in test_cases:
+                assert_equal(np.maximum([v1], [v2]), [expected])
+                assert_equal(np.maximum.reduce([v1, v2]), expected)
+
 
 class TestMinimum(_FilterInvalids):
     def test_reduce(self):
@@ -1799,6 +1822,30 @@ def test_strided_array(self):
         assert_equal(np.minimum(arr1[:6:2], arr2[::3], out=out[::3]), np.array([-4.0, 1.0, np.nan]))
         assert_equal(out, out_mintrue)
 
+    def test_precision(self):
+        dtypes = [np.float16, np.float32, np.float64]
+        if(hasattr(np, 'float128')):
+            dtypes.append(np.float128)
+
+        for dt in dtypes:
+            dtmin = np.finfo(dt).min
+            dtmax = np.finfo(dt).max
+            d1 = dt(0.1)
+            d1_next = np.nextafter(d1, np.inf)
+
+            test_cases = [
+                # v1    v2          expected
+                (dtmin, np.inf,     dtmin),
+                (dtmax, np.inf,     dtmax),
+                (d1,    d1_next,    d1),
+                (dtmin, np.nan,     np.nan),
+            ]
+
+            for v1, v2, expected in test_cases:
+                assert_equal(np.minimum([v1], [v2]), [expected])
+                assert_equal(np.minimum.reduce([v1, v2]), expected)
+
+
 class TestFmax(_FilterInvalids):
     def test_reduce(self):
         dflt = np.typecodes['AllFloat']
@@ -1840,6 +1887,29 @@ def test_complex_nans(self):
             out = np.array([0,    0, nan], dtype=complex)
             assert_equal(np.fmax(arg1, arg2), out)
 
+    def test_precision(self):
+        dtypes = [np.float16, np.float32, np.float64]
+        if(hasattr(np, 'float128')):
+            dtypes.append(np.float128)
+
+        for dt in dtypes:
+            dtmin = np.finfo(dt).min
+            dtmax = np.finfo(dt).max
+            d1 = dt(0.1)
+            d1_next = np.nextafter(d1, np.inf)
+
+            test_cases = [
+                # v1    v2          expected
+                (dtmin, -np.inf,    dtmin),
+                (dtmax, -np.inf,    dtmax),
+                (d1,    d1_next,    d1_next),
+                (dtmax, np.nan,     dtmax),
+            ]
+
+            for v1, v2, expected in test_cases:
+                assert_equal(np.fmax([v1], [v2]), [expected])
+                assert_equal(np.fmax.reduce([v1, v2]), expected)
+
 
 class TestFmin(_FilterInvalids):
     def test_reduce(self):
@@ -1882,6 +1952,30 @@ def test_complex_nans(self):
             out = np.array([0,    0, nan], dtype=complex)
             assert_equal(np.fmin(arg1, arg2), out)
 
+    def test_precision(self):
+        dtypes = [np.float16, np.float32, np.float64]
+        if(hasattr(np, 'float128')):
+            dtypes.append(np.float128)
+
+        for dt in dtypes:
+            dtmin = np.finfo(dt).min
+            dtmax = np.finfo(dt).max
+            d1 = dt(0.1)
+            d1_next = np.nextafter(d1, np.inf)
+
+            test_cases = [
+                # v1    v2          expected
+                (dtmin, np.inf,     dtmin),
+                (dtmax, np.inf,     dtmax),
+                (d1,    d1_next,    d1),
+                (dtmin, np.nan,     dtmin),
+            ]
+
+            for v1, v2, expected in test_cases:
+                assert_equal(np.fmin([v1], [v2]), [expected])
+                assert_equal(np.fmin.reduce([v1, v2]), expected)
+
+
 
 class TestBool:
     def test_exceptions(self):


Code review comments:
diff_hunk: @@ -1840,6 +1887,29 @@ def test_complex_nans(self):
             out = np.array([0,    0, nan], dtype=complex)
             assert_equal(np.fmax(arg1, arg2), out)
 
+    def test_precision(self):
+        dtypes = [np.float16, np.float32, np.float64]
+        if(hasattr(np, 'float128')):
+            dtypes.append(np.float128)
path: numpy/core/tests/test_umath.py
text: @user1:
```suggestion
        dtypes = [np.float16, np.float32, np.float64, np.longdouble]
```

diff_hunk: @@ -1728,6 +1728,29 @@ def test_strided_array(self):
         assert_equal(np.maximum(arr1[:6:2], arr2[::3], out=out[::3]), np.array([-2.0, 10., np.nan]))
         assert_equal(out, out_maxtrue)
 
+    def test_precision(self):
+        dtypes = [np.float16, np.float32, np.float64]
+        if(hasattr(np, 'float128')):
+            dtypes.append(np.float128)
path: numpy/core/tests/test_umath.py
text: @user1:
```suggestion
        dtypes = [np.float16, np.float32, np.float64, np.longdouble]
```

diff_hunk: @@ -1799,6 +1822,30 @@ def test_strided_array(self):
         assert_equal(np.minimum(arr1[:6:2], arr2[::3], out=out[::3]), np.array([-4.0, 1.0, np.nan]))
         assert_equal(out, out_mintrue)
 
+    def test_precision(self):
+        dtypes = [np.float16, np.float32, np.float64]
+        if(hasattr(np, 'float128')):
+            dtypes.append(np.float128)
path: numpy/core/tests/test_umath.py
text: @user1:
```suggestion
        dtypes = [np.float16, np.float32, np.float64, np.longdouble]
```

diff_hunk: @@ -1882,6 +1952,30 @@ def test_complex_nans(self):
             out = np.array([0,    0, nan], dtype=complex)
             assert_equal(np.fmin(arg1, arg2), out)
 
+    def test_precision(self):
+        dtypes = [np.float16, np.float32, np.float64]
+        if(hasattr(np, 'float128')):
+            dtypes.append(np.float128)
path: numpy/core/tests/test_umath.py
text: @user1:
```suggestion
        dtypes = [np.float16, np.float32, np.float64, np.longdouble]
```