Fix failing unit test in Python 3.13 (#288)

* Fix failing unit test and code in `length` (math_util.c)  in Python 3.13

* add a code comment for PyUFunc_clearfperr

CHANGES.rst

@@ -4,6 +4,14 @@
 Release Notes
 =============
 
+1.3.3 (30-January-2025)
+=======================
+
+- Fixed a bug in the C implementation of ``dot_qd`` that may result in
+  "invalid value encountered in length" runtime warning when input vectors
+  contain NaN values. [#288]
+
+
 1.3.2 (12-June-2024)
 ====================
 
@@ -16,6 +24,7 @@ Release Notes
 
 - Build wheel with Numpy 2.0 release candidate ahead of Numpy 2.0 release [#279]
 
+
 1.3.1 (5-November-2023)
 =======================
 

spherical_geometry/tests/test_basic.py

@@ -529,8 +529,27 @@ def test_math_util_angle_domain():
 
 
 def test_math_util_length_domain():
+    a = [
+        [np.nan, 0, 0],
+        [np.nan, 0, 0],
+        [np.nan, np.nan, np.nan],
+        [0, 0, 0],
+        [0, 0, 0],
+        [0, 0, np.nan],
+        [0, 0, 0],
+    ]
+    b = [
+        [0, 0, np.inf],
+        [0, 0, np.inf],
+        [0, 0, 0],
+        [np.inf, np.inf, np.inf],
+        [0, 0, np.inf],
+        [0, 0, 0],
+        [0, 0, 0],
+    ]
+
     with pytest.raises(ValueError):
-        great_circle_arc.length([[np.nan, 0, 0]], [[0, 0, np.inf]])
+        great_circle_arc.length(a, b)
 
 
 def test_math_util_angle_nearly_coplanar_vec():

src/math_util.c

@@ -181,14 +181,23 @@ dot_qd(const qd *A, const qd *B, qd *C) {
 
     c_qd_add(tmp[0], tmp[1], tmp[3]);
     c_qd_add(tmp[3], tmp[2], C->x);
+
+    // In Python 3.13 it seems that the code above sets a floating point error
+    // flag (when input vectors contain nan/inf values) and this raises a
+    // warning/error "RuntimeWarning: invalid value encountered in length"
+    // and which results in a SystemError once
+    // PyErr_SetString(PyExc_ValueError, "Out of domain for acos") is called
+    // in length_qd. This clears FP error flags before raising the above
+    // exception.
+    // Also See https://github.com/spacetelescope/spherical_geometry/pull/288
+    PyUFunc_clearfperr();
 }
 
 /*
     normalized_dot_qd returns dot product of normalized input vectors.
 */
 static NPY_INLINE int
 normalized_dot_qd(const qd *A, const qd *B, qd *dot_val) {
-    int flag;
     qd aa, bb, ab;
     double aabb[4];
     double norm[4];
@@ -206,7 +215,7 @@ normalized_dot_qd(const qd *A, const qd *B, qd *dot_val) {
         return 1;
     }
 
-    flag = c_qd_sqrt(aabb, norm);
+    c_qd_sqrt(aabb, norm);
 
     if (norm[0] == 0.0) {
         /* return non-normalized value: */
@@ -646,11 +655,15 @@ DOUBLE_length(char **args, const intp *dimensions, const intp *steps, void *NPY_
         load_point_qd(ip1, is1, A);
         load_point_qd(ip2, is2, B);
 
-        if (normalize_qd(A, A)) return;
-        if (normalize_qd(B, B)) return;
-
-        if (length_qd(A, B, &s)) return;
-
+        if (normalize_qd(A, A)) {
+            return;
+        }
+        if (normalize_qd(B, B)) {
+            return;
+        }
+        if (length_qd(A, B, &s)) {
+            return;
+        }
         *((double *)op) = s.x[0];
     END_OUTER_LOOP