diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index 2eac92ee..40ee029b 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/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 *); diff --git a/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp index 09c11a07..c2710747 100644 --- a/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp @@ -385,7 +385,10 @@ init_quad_unary_ops(PyObject *numpy) if (create_quad_unary_ufunc(numpy, "deg2rad") < 0) { return -1; } - + if (create_quad_unary_ufunc(numpy, "spacing") < 0) { + return -1; + } + // Logical operations (unary: not) if (create_quad_logical_not_ufunc(numpy, "logical_not") < 0) { return -1; diff --git a/quaddtype/release_tracker.md b/quaddtype/release_tracker.md index 6b374ecb..c5e1ad2b 100644 --- a/quaddtype/release_tracker.md +++ b/quaddtype/release_tracker.md @@ -78,7 +78,7 @@ | signbit | ✅ | ✅ | | copysign | ✅ | ✅ | | nextafter | ✅ | ✅ | -| spacing | | | +| spacing | ✅ | ✅ | | modf | | | | ldexp | | | | frexp | | | diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index b919ea3b..693b9090 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/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" +