Merge pull request #191 from SwayamInSync/spacing

Author: SwayamInSync

quaddtype/numpy_quaddtype/src/ops.hpp

@@ -1008,6 +1008,43 @@ quad_nextafter(const Sleef_quad *x, const Sleef_quad *y)
     return quad_set_words64(hx, lx);
 }
 
+static inline Sleef_quad
+quad_spacing(const Sleef_quad *x)
+{   
+    // spacing(x) returns the distance between x and the next representable value
+    // The result has the SAME SIGN as x (NumPy convention)
+    // For x >= 0: spacing(x) = nextafter(x, +inf) - x
+    // For x < 0:  spacing(x) = nextafter(x, -inf) - x (negative result)
+    
+    // Handle NaN
+    if (Sleef_iunordq1(*x, *x)) {
+        return *x;  // NaN
+    }
+    
+    // Handle infinity -> NaN (numpy convention)
+    if (quad_isinf(x)) {
+        return QUAD_NAN;
+    }
+    
+    // Determine direction based on sign of x
+    Sleef_quad direction;
+    if (Sleef_icmpltq1(*x, QUAD_ZERO)) {
+        // Negative: move toward -inf
+        direction = Sleef_negq1(QUAD_POS_INF);
+    } else {
+        // Positive or zero: move toward +inf
+        direction = QUAD_POS_INF;
+    }
+    
+    // Compute nextafter(x, direction)
+    Sleef_quad next = quad_nextafter(x, &direction);
+    
+    // spacing = next - x (preserves sign)
+    Sleef_quad result = Sleef_subq1_u05(next, *x);
+    
+    return result;
+}
+
 // Binary long double operations
 typedef long double (*binary_op_longdouble_def)(const long double *, const long double *);
 // Binary long double operations with 2 outputs (for divmod, modf, frexp)
@@ -1292,6 +1329,38 @@ ld_nextafter(const long double *x1, const long double *x2)
     return nextafterl(*x1, *x2);
 }
 
+static inline long double
+ld_spacing(const long double *x)
+{    
+    // Handle NaN
+    if (isnan(*x)) {
+        return *x;  // NaN
+    }
+    
+    // Handle infinity -> NaN (numpy convention)
+    if (isinf(*x)) {
+        return NAN;
+    }
+    
+    // Determine direction based on sign of x
+    long double direction;
+    if (*x < 0.0L) {
+        // Negative: move toward -inf
+        direction = -INFINITY;
+    } else {
+        // Positive or zero: move toward +inf
+        direction = INFINITY;
+    }
+    
+    // Compute nextafter(x, direction)
+    long double next = nextafterl(*x, direction);
+    
+    // spacing = next - x (preserves sign)
+    long double result = next - (*x);
+    
+    return result;
+}
+
 // comparison quad functions
 typedef npy_bool (*cmp_quad_def)(const Sleef_quad *, const Sleef_quad *);
 

quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp

@@ -385,7 +385,10 @@ init_quad_unary_ops(PyObject *numpy)
     if (create_quad_unary_ufunc<quad_radians, ld_radians>(numpy, "deg2rad") < 0) {
         return -1;
     }
-    
+    if (create_quad_unary_ufunc<quad_spacing, ld_spacing>(numpy, "spacing") < 0) {
+        return -1;
+    }
+
     // Logical operations (unary: not)
     if (create_quad_logical_not_ufunc<quad_logical_not, ld_logical_not>(numpy, "logical_not") < 0) {
         return -1;

quaddtype/release_tracker.md

@@ -78,7 +78,7 @@
 | signbit       | ✅    | ✅                                                                   |
 | copysign      | ✅    | ✅                                                                   |
 | nextafter     | ✅    | ✅                                                                   |
-| spacing       |       |                                                                      |
+| spacing       | ✅    | ✅                                                                   |
 | modf          |       |                                                                      |
 | ldexp         |       |                                                                      |
 | frexp         |       |                                                                      |

quaddtype/tests/test_quaddtype.py

@@ -2256,3 +2256,204 @@ def test_direction(self):
             elif expected_dir == "less":
                 assert result < q_x, f"nextafter({x}, {y}) should be < {x}, got {result}"
 
+class TestSpacing:
+    """Test cases for np.spacing function with QuadPrecision dtype"""
+    
+    @pytest.mark.parametrize("x", [
+        np.nan, -np.nan,
+        np.inf, -np.inf,
+    ])
+    def test_special_values_return_nan(self, x):
+        """Test spacing with NaN and infinity inputs returns NaN"""
+        q_x = QuadPrecision(x)
+        result = np.spacing(q_x)
+        
+        assert isinstance(result, QuadPrecision)
+        assert np.isnan(result), f"spacing({x}) should be NaN, got {result}"
+    
+    @pytest.mark.parametrize("x", [
+        1.0, -1.0,
+        10.0, -10.0,
+        100.0, -100.0,
+    ])
+    def test_sign_preservation(self, x):
+        """Test that spacing preserves the sign of the input"""
+        q_x = QuadPrecision(x)
+        result = np.spacing(q_x)
+        
+        q_zero = QuadPrecision(0)
+        # spacing should have the same sign as x
+        if x > 0:
+            assert result > q_zero, f"spacing({x}) should be positive, got {result}"
+        else:
+            assert result < q_zero, f"spacing({x}) should be negative, got {result}"
+        
+        # Compare with numpy behavior
+        np_result = np.spacing(np.float64(x))
+        assert np.signbit(result) == np.signbit(np_result), \
+            f"Sign mismatch for spacing({x}): quad signbit={np.signbit(result)}, numpy signbit={np.signbit(np_result)}"
+    
+    @pytest.mark.parametrize("x", [
+        0.0,
+        -0.0,
+    ])
+    def test_zero(self, x):
+        """Test spacing of zero returns smallest_subnormal"""
+        q_x = QuadPrecision(x)
+        result = np.spacing(q_x)
+        
+        finfo = np.finfo(QuadPrecDType())
+        q_zero = QuadPrecision(0)
+        
+        # spacing(±0.0) should return smallest_subnormal (positive)
+        assert result == finfo.smallest_subnormal, \
+            f"spacing({x}) should be smallest_subnormal, got {result}"
+        assert result > q_zero, f"spacing({x}) should be positive, got {result}"
+    
+    @pytest.mark.parametrize("x", [
+        1.0,
+        -1.0,
+    ])
+    def test_one_and_negative_one(self, x):
+        """Test spacing(±1.0) equals ±machine epsilon"""
+        q_x = QuadPrecision(x)
+        result = np.spacing(q_x)
+        
+        finfo = np.finfo(QuadPrecDType())
+        q_zero = QuadPrecision(0)
+        
+        # For binary floating point, spacing(±1.0) = ±eps
+        expected = finfo.eps if x > 0 else -finfo.eps
+        assert result == expected, \
+            f"spacing({x}) should equal {expected}, got {result}"
+        
+        if x > 0:
+            assert result > q_zero, "spacing(1.0) should be positive"
+        else:
+            assert result < q_zero, "spacing(-1.0) should be negative"
+    
+    @pytest.mark.parametrize("x", [
+        1.0, -1.0,
+        2.0, -2.0,
+        10.0, -10.0,
+        100.0, -100.0,
+        1e10, -1e10,
+        1e-10, -1e-10,
+        0.5, -0.5,
+        0.25, -0.25,
+    ])
+    def test_spacing_is_non_zero(self, x):
+        """Test that spacing always has positive magnitude"""
+        q_x = QuadPrecision(x)
+        result = np.spacing(q_x)
+        
+        q_zero = QuadPrecision(0)
+        # The absolute value should be positive
+        abs_result = np.abs(result)
+        assert abs_result > q_zero, f"|spacing({x})| should be positive, got {abs_result}"
+    
+    def test_smallest_subnormal(self):
+        """Test spacing at smallest subnormal value"""
+        finfo = np.finfo(QuadPrecDType())
+        smallest = finfo.smallest_subnormal
+        
+        result = np.spacing(smallest)
+        
+        q_zero = QuadPrecision(0)
+        # spacing(smallest_subnormal) should be smallest_subnormal itself
+        # (it's the minimum spacing in the subnormal range)
+        assert result == smallest, \
+            f"spacing(smallest_subnormal) should be smallest_subnormal, got {result}"
+        assert result > q_zero, "Result should be positive"
+    
+    @pytest.mark.parametrize("x", [
+        1.5, -1.5,
+        3.7, -3.7,
+        42.0, -42.0,
+        1e8, -1e8,
+    ])
+    def test_finite_values(self, x):
+        """Test spacing on various finite values"""
+        q_x = QuadPrecision(x)
+        result = np.spacing(q_x)
+        
+        q_zero = QuadPrecision(0)
+        # Result should be finite
+        assert np.isfinite(result), \
+            f"spacing({x}) should be finite, got {result}"
+        
+        # Result should be non-zero
+        assert result != q_zero, \
+            f"spacing({x}) should be non-zero, got {result}"
+        
+        # Result should have same sign as input
+        assert np.signbit(result) == np.signbit(q_x), \
+            f"spacing({x}) should have same sign as {x}"
+    
+    def test_array_spacing(self):
+        """Test spacing on an array of QuadPrecision values"""
+        values = [1.0, -1.0, 2.0, -2.0, 0.0, 10.0, -10.0]
+        q_array = np.array([QuadPrecision(v) for v in values])
+        
+        result = np.spacing(q_array)
+        
+        q_zero = QuadPrecision(0)
+        # Check each result
+        for i, val in enumerate(values):
+            q_val = QuadPrecision(val)
+            if val != 0:
+                assert np.signbit(result[i]) == np.signbit(q_val), \
+                    f"Sign mismatch for spacing({val})"
+            else:
+                assert result[i] > q_zero, \
+                    f"spacing(0) should be positive, got {result[i]}"
+    
+    @pytest.mark.parametrize("x", [
+        1e-100, -1e-100,
+        1e-200, -1e-200,
+    ])
+    def test_subnormal_range(self, x):
+        """Test spacing in subnormal range"""
+        q_x = QuadPrecision(x)
+        result = np.spacing(q_x)
+        
+        finfo = np.finfo(QuadPrecDType())
+        
+        # In subnormal range, spacing should be smallest_subnormal
+        # or at least very small
+        assert np.abs(result) >= finfo.smallest_subnormal, \
+            f"spacing({x}) should be >= smallest_subnormal"
+        
+        q_zero = QuadPrecision(0)
+        # Sign should match (but for very small subnormals, spacing might underflow to zero)
+        if result != q_zero:
+            assert np.signbit(result) == np.signbit(q_x), \
+                f"spacing({x}) should have same sign as {x}"
+
+    def test_smallest_normal_spacing(self):
+        """Test spacing for smallest normal value and 2*smallest normal"""
+        finfo = np.finfo(QuadPrecDType())
+        smallest_normal = finfo.smallest_normal
+        
+        # Test spacing at smallest normal value
+        result1 = np.spacing(smallest_normal)
+        
+        # Test spacing at 2 * smallest normal value
+        two_smallest_normal = 2 * smallest_normal
+        result2 = np.spacing(two_smallest_normal)
+        
+        q_zero = QuadPrecision("0")
+        
+        # spacing(smallest_normal) should be smallest_subnormal
+        # (the spacing at the boundary between subnormal and normal range)
+        expected1 = finfo.smallest_subnormal
+        assert result1 == expected1, \
+            f"spacing(smallest_normal) should be {expected1}, got {result1}"
+        assert result1 > q_zero, "Result should be positive"
+        
+        # The scaling relationship: spacing(2*x) = 2*spacing(x) for normal numbers
+        expected2 = 2 * finfo.smallest_subnormal
+        assert result2 == expected2, \
+            f"spacing(2*smallest_normal) should be {expected2}, got {result2}"
+        assert result2 > q_zero, "Result should be positive"
+