`gcd()` hangs in infinite loop

[aerobio]

Describe the issue:

When receiving infinity as a (theoretically illegal) argument, gcd() hangs instead of throwing an exception. In this case, infinity is an element of an array of object dtype.

Reproduce the code example:

import numpy as np

inf = np.array([np.inf], dtype=np.object_)
gcd = np.gcd(inf, 33)

Error message:

No error. Function `gcd()` remains in an infinite loop.

Python and NumPy Versions:

Numpy: 1.26.4 Python: 3.10.12

Runtime Environment:

[{'numpy_version': '1.26.4', 'python': '3.10.12 (main, Nov 20 2023, 15:14:05) [GCC 11.4.0]', 'uname': uname_result(system='Linux', node='VA-HNbs0rMjnJlA', release='5.15.146.1-microsoft-standard-WSL2', version='#1 SMP Thu Jan 11 04:09:03 UTC 2024', machine='x86_64')}, {'simd_extensions': {'baseline': ['SSE', 'SSE2', 'SSE3'], 'found': ['SSSE3', 'SSE41', 'POPCNT', 'SSE42', 'AVX', 'F16C', 'FMA3', 'AVX2'], 'not_found': ['AVX512F', 'AVX512CD', 'AVX512_KNL', 'AVX512_KNM', 'AVX512_SKX', 'AVX512_CLX', 'AVX512_CNL', 'AVX512_ICL']}}, {'architecture': 'Haswell', 'filepath': '[REDACTED]/.venv/lib/python3.10/site-packages/numpy.libs/libopenblas64_p-r0-0cf96a72.3.23.dev.so', 'internal_api': 'openblas', 'num_threads': 12, 'prefix': 'libopenblas', 'threading_layer': 'pthreads', 'user_api': 'blas', 'version': '0.3.23.dev'}]

Context for the issue:

It should raise an exception (as np.gcd(np.inf, 33) does) instead of hanging in an infinite loop. Infinite loops are bad for your cloud bill.

[ngoldbaum]

I can reproduce this on main. Here's a C traceback from lldb if I do ctrl-C during the hang:

* thread #1, queue = 'com.apple.main-thread', stop reason = signal SIGSTOP
  * frame #0: 0x0000000100c654e8 libpython3.13t.dylib`fmod + 4
    frame #1: 0x0000000100a48e14 libpython3.13t.dylib`float_rem + 368
    frame #2: 0x0000000100a07c58 libpython3.13t.dylib`binary_op1 + 264
    frame #3: 0x0000000100a081dc libpython3.13t.dylib`PyNumber_Remainder + 32
    frame #4: 0x0000000100b5e3b8 libpython3.13t.dylib`_PyEval_EvalFrameDefault + 1928
    frame #5: 0x0000000100a2582c libpython3.13t.dylib`_PyObject_CallFunctionVa + 152
    frame #6: 0x0000000100a25788 libpython3.13t.dylib`PyObject_CallFunction + 56
    frame #7: 0x0000000100696c94 _multiarray_umath.cpython-313t-darwin.so`npy_ObjectGCD + 120
    frame #8: 0x0000000100544558 _multiarray_umath.cpython-313t-darwin.so`PyUFunc_OO_O + 140
    frame #9: 0x000000010067b614 _multiarray_umath.cpython-313t-darwin.so`generic_wrapped_legacy_loop + 40
    frame #10: 0x0000000100681e4c _multiarray_umath.cpython-313t-darwin.so`PyUFunc_GenericFunctionInternal + 2440
    frame #11: 0x000000010068012c _multiarray_umath.cpython-313t-darwin.so`ufunc_generic_fastcall + 3508
    frame #12: 0x0000000100a24e3c libpython3.13t.dylib`PyObject_Vectorcall + 88
    frame #13: 0x0000000100b6022c libpython3.13t.dylib`_PyEval_EvalFrameDefault + 9724
    frame #14: 0x0000000100b5d928 libpython3.13t.dylib`PyEval_EvalCode + 360
... further frames elided ...

And it actually turns out this calls into a pure-python function:

https://github.com/numpy/numpy/blob/a4cddb60489f821a1a4dffc16cd5c69755d43bdb/numpy/_core/_internal.py#L861-L865

Not surprisingly, it causes an infinite loop if you pass it infinity and zero. Should probably catch any edge cases that would cause the loop to never terminate. Or maybe deprecate float inputs while we're at it? Python doesn't allow it:

>>> math.gcd(float('inf'), 0)
Traceback (most recent call last):
  File "<python-input-2>", line 1, in <module>
    math.gcd(float('inf'), 0)
    ~~~~~~~~^^^^^^^^^^^^^^^^^
TypeError: 'float' object cannot be interpreted as an integer

[JuliaPoo]

Is this the kind of change you're going for when you mentioned depreciation? This will depreciate arbitrary objects that numpy can't cast to int, and are not accepted by math.gcd.

diff --git a/numpy/_core/src/umath/funcs.inc.src b/numpy/_core/src/umath/funcs.inc.src
index 3825bd8694..383c42e4de 100644
--- a/numpy/_core/src/umath/funcs.inc.src
+++ b/numpy/_core/src/umath/funcs.inc.src
@@ -177,35 +177,8 @@ npy_ObjectTrunc(PyObject *obj) {
 static PyObject *
 npy_ObjectGCD(PyObject *i1, PyObject *i2)
 {
-    PyObject *gcd = NULL;
-
-    /* use math.gcd if valid on the provided types */
-    {
-        gcd = PyObject_CallFunction(npy_static_pydata.math_gcd_func,
-                                    "OO", i1, i2);
-        if (gcd != NULL) {
-            return gcd;
-        }
-        /* silence errors, and fall back on pure-python gcd */
-        PyErr_Clear();
-    }
-
-    /* otherwise, use our internal one, written in python */
-    {
-        npy_cache_import("numpy._core._internal", "_gcd",
-                         &npy_thread_unsafe_state.internal_gcd_func);
-        if (npy_thread_unsafe_state.internal_gcd_func == NULL) {
-            return NULL;
-        }
-        gcd = PyObject_CallFunction(npy_thread_unsafe_state.internal_gcd_func,
+    return PyObject_CallFunction(npy_static_pydata.math_gcd_func,
                                     "OO", i1, i2);
-        if (gcd == NULL) {
-            return NULL;
-        }
-        /* _gcd has some unusual behaviour regarding sign */
-        Py_SETREF(gcd, PyNumber_Absolute(gcd));
-        return gcd;
-    }
 }

Object arrays/arrays containing objects are already handled by PyUFunc_SimpleUniformOperationTypeResolver which uses safe casting.

Update:

The above change fails with Decimal in the test suite, which was what the above code was intending to support

[ngoldbaum]

Hey, sorry for taking a couple days to get back to you. This week has been hectic! I was more thinking of adding some type-checking to the pure-python GCD, and explicitly rejecting float inputs there. We'd still need to keep the C code with the fallback to pure python to support decimals apparently. I had no idea that test was there when I was debugging this issue and proposed my possible fix.

[JuliaPoo]

Oh my bad, I missed this. Currently it's a little unclear the kind of type checking that would be needed. I think the core problem here is that np.gcd is written to handle object arrays containing objects that can't be casted to an integer type via safe casting. Will it be okay to just not fix this and the onus for a sensible output lies on the user?

[ngoldbaum]

Well, one thing I think we can do is ban infinity and NaN at least, to avoid the infinite loop:

diff --git a/numpy/_core/_internal.py b/numpy/_core/_internal.py
index 058e93644d..fcb6c37f0d 100644
--- a/numpy/_core/_internal.py
+++ b/numpy/_core/_internal.py
@@ -5,6 +5,7 @@

 """
 import ast
+import math
 import re
 import sys
 import warnings
@@ -860,6 +861,9 @@ def _prod(a):

 def _gcd(a, b):
     """Calculate the greatest common divisor of a and b"""
+    if not (math.isfinite(a) and math.isfinite(b)):
+        raise ValueError('Can only find greatest common demoninator of a '
+                         f'finite argument, found "{a}" and "{b}"')
     while b:
         a, b = b, a % b
     return a

You could add isinstance checks to deprecate floats I guess but I think banning infinity and nan is probably sufficient to call this issue fixed. If people really want to find the GCD of floats I guess they can continue to use dtype=object.

[guan404ming]

Hi @ngoldbaum~ I had opened a PR based on ur solution. Could you please help me to review it? Thank you~

[WarrenWeckesser]

Fixed by https://github.com/numpy/numpy/pull/27014.