Fix np.mod(0, inf)

src/include/ops.hpp

@@ -629,6 +629,11 @@ quad_mod(const Sleef_quad *a, const Sleef_quad *b)
 
     // finite % inf
     if (quad_isfinite(a) && quad_isinf(b)) {
+        // 0 % inf
+        if (Sleef_icmpeqq1(*a, QUAD_PRECISION_ZERO)) {
+            return Sleef_copysignq1(*a, *b);
+        }
+
         int sign_a = quad_signbit(a);
         int sign_b = quad_signbit(b);
 
@@ -1241,6 +1246,10 @@ ld_mod(const long double *a, const long double *b)
         return NAN;
 
     if (isfinite(*a) && isinf(*b)) {
+        // 0 % inf
+        if (*a == 0.0L) {
+            return copysignl(*a, *b);
+        }
         int sign_a = signbit(*a);
         int sign_b = signbit(*b);
         return (sign_a == sign_b) ? *a : *b;

tests/test_quaddtype.py

@@ -2629,7 +2629,8 @@ def test_array_operations():
     # Finite % infinity cases
     (5.0, float('inf')), (-5.0, float('inf')),
     (5.0, float('-inf')), (-5.0, float('-inf')),
-    (0.0, float('inf')), (-0.0, float('-inf')),
+    (0.0, float('inf')), (-0.0, float('inf')),
+    (0.0, float('-inf')), (-0.0, float('-inf')),
 
     # NaN cases (should return NaN)
     (float('nan'), 3.0), (3.0, float('nan')), (float('nan'), float('nan')),
@@ -2674,6 +2675,7 @@ def test_mod(a, b, backend, op):
     if numpy_result == 0.0:
         numpy_sign = np.signbit(numpy_result)
         quad_sign = np.signbit(quad_result)
+        assert quad_result == 0, f"Zero mismatch for {a} % {b}: numpy={numpy_result}, quad={quad_result}"
         assert numpy_sign == quad_sign, f"Zero sign mismatch for {a} % {b}: numpy={numpy_sign}, quad={quad_sign}"
 
     # Check that non-zero results have correct sign relative to divisor