From d1f9db0c4a0c8322ff2cfdb39ca9910e6c277a34 Mon Sep 17 00:00:00 2001 From: SwayamInSync Date: Sun, 19 Oct 2025 19:37:30 +0000 Subject: [PATCH 1/2] works on my machine :) --- quaddtype/numpy_quaddtype/src/ops.hpp | 131 ++++++++++++ .../numpy_quaddtype/src/umath/binary_ops.cpp | 3 + quaddtype/release_tracker.md | 2 +- quaddtype/tests/test_quaddtype.py | 191 ++++++++++++++++++ 4 files changed, 326 insertions(+), 1 deletion(-) diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index c76904d..a8c4ef3 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/quaddtype/numpy_quaddtype/src/ops.hpp @@ -883,6 +883,131 @@ quad_hypot(const Sleef_quad *x1, const Sleef_quad *x2) return Sleef_hypotq1_u05(*x1, *x2); } +// todo: we definitely need to refactor this file, getting too clumsy everything here + +static inline void quad_get_words64(int64_t *hx, uint64_t *lx, Sleef_quad x) +{ + union { + Sleef_quad q; + struct { +#if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_BIG_ENDIAN__) + uint64_t hi; + uint64_t lo; +#else + uint64_t lo; + uint64_t hi; +#endif + } i; + } u; + u.q = x; + *hx = (int64_t)u.i.hi; + *lx = u.i.lo; +} + +static inline Sleef_quad quad_set_words64(int64_t hx, uint64_t lx) +{ + union { + Sleef_quad q; + struct { +#if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_BIG_ENDIAN__) + uint64_t hi; + uint64_t lo; +#else + uint64_t lo; + uint64_t hi; +#endif + } i; + } u; + u.i.hi = (uint64_t)hx; + u.i.lo = lx; + return u.q; +} + + +static inline Sleef_quad +quad_nextafter(const Sleef_quad *x, const Sleef_quad *y) +{ + int64_t hx, hy, ix, iy; + uint64_t lx, ly; + + quad_get_words64(&hx, &lx, *x); + quad_get_words64(&hy, &ly, *y); + + // extracting absolute value + ix = hx & 0x7fffffffffffffffLL; + iy = hy & 0x7fffffffffffffffLL; + + // NaN if either is NaN + if (Sleef_iunordq1(*x, *y)) { + return Sleef_addq1_u05(*x, *y); // still NaN + } + + // x == y then return y + if (Sleef_icmpeqq1(*x, *y)) { + return *y; + } + + // both input 0 then extract sign from y and return correspondingly signed smallest subnormal + if ((ix | lx) == 0) { + Sleef_quad result = quad_set_words64(hy & 0x8000000000000000LL, 1); // quad_set_words64(sign_y, 1) + return result; + } + + if (hx >= 0) + { + if (hx > hy || ((hx == hy) && (lx > ly))) + { + // Moving toward smaller y (x > y) + // low word is 0 then decrement high word first (borrowing) + if (lx == 0) + hx--; + lx--; + } + else + { + lx++; + if (lx == 0) + // carry to high words + hx++; + } + } + else + { + // Moving toward larger y + // similar to above case just direction will be swapped + if (hy >= 0 || hx > hy || ((hx == hy) && (lx > ly))) + { + if (lx == 0) + hx--; + lx--; + } + else + { + lx++; + if (lx == 0) + hx++; + } + } + + // check if reached infinity + // this can be NaN XOR inf but NaN are already checked at start + hy = hx & 0x7fff000000000000LL; + if (hy == 0x7fff000000000000LL) { + Sleef_quad result = quad_set_words64(hx, lx); + return result; + } + // check whether entered into subnormal range + // 0 exponent i.e. either (0 or subnormal) + if (hy == 0) { + Sleef_quad result = quad_set_words64(hx, lx); + return result; + } + + // well I did not read to have those above checks + // but they can be important when settng FPE flag manually + return quad_set_words64(hx, lx); +} + // 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) @@ -1161,6 +1286,12 @@ ld_hypot(const long double *x1, const long double *x2) return hypotl(*x1, *x2); } +static inline long double +ld_nextafter(const long double *x1, const long double *x2) +{ + return nextafterl(*x1, *x2); +} + // comparison quad functions typedef npy_bool (*cmp_quad_def)(const Sleef_quad *, const Sleef_quad *); diff --git a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp index 9cc19a9..562f277 100644 --- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp @@ -448,6 +448,9 @@ init_quad_binary_ops(PyObject *numpy) if (create_quad_binary_ufunc(numpy, "copysign") < 0) { return -1; } + if (create_quad_binary_ufunc(numpy, "nextafter") < 0) { + return -1; + } if (create_quad_binary_ufunc(numpy, "logaddexp") < 0) { return -1; } diff --git a/quaddtype/release_tracker.md b/quaddtype/release_tracker.md index 267fee9..6b374ec 100644 --- a/quaddtype/release_tracker.md +++ b/quaddtype/release_tracker.md @@ -77,7 +77,7 @@ | isnan | ✅ | ✅ | | signbit | ✅ | ✅ | | copysign | ✅ | ✅ | -| nextafter | | | +| nextafter | ✅ | ✅ | | spacing | | | | modf | | | | ldexp | | | diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index ec22dd5..b919ea3 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -2065,3 +2065,194 @@ def test_radians(op, degrees, expected_radians): else: np.testing.assert_allclose(float(result), expected_radians, rtol=1e-13) + +class TestNextAfter: + """Test cases for np.nextafter function with QuadPrecision dtype""" + + @pytest.mark.parametrize("x1,x2", [ + # NaN tests + (np.nan, 1.0), + (1.0, np.nan), + (np.nan, np.nan), + ]) + def test_nan(self, x1, x2): + """Test nextafter with NaN inputs returns NaN""" + q_x1 = QuadPrecision(x1) + q_x2 = QuadPrecision(x2) + + result = np.nextafter(q_x1, q_x2) + + assert isinstance(result, QuadPrecision) + assert np.isnan(float(result)) + + def test_precision(self): + """Test that nextafter provides the exact next representable value""" + # Start with 1.0 and move towards 2.0 + x1 = QuadPrecision(1.0) + x2 = QuadPrecision(2.0) + + result = np.nextafter(x1, x2) + + # Get machine epsilon from finfo + finfo = np.finfo(QuadPrecDType()) + expected = x1 + finfo.eps + + # result should be exactly 1.0 + eps + assert result == expected + + # Moving the other direction should give us back 1.0 + result_back = np.nextafter(result, x1) + assert result_back == x1 + + def test_smallest_subnormal(self): + """Test that nextafter(0.0, 1.0) returns the smallest positive subnormal (TRUE_MIN)""" + zero = QuadPrecision(0.0) + one = QuadPrecision(1.0) + + result = np.nextafter(zero, one) # smallest_subnormal + finfo = np.finfo(QuadPrecDType()) + + assert result == finfo.smallest_subnormal, \ + f"nextafter(0.0, 1.0) should equal smallest_subnormal, got {result} vs {finfo.smallest_subnormal}" + + # Verify it's positive and very small + assert result > zero, "nextafter(0.0, 1.0) should be positive" + + # Moving back towards zero should give us zero + result_back = np.nextafter(result, zero) + assert result_back == zero, f"nextafter(smallest_subnormal, 0.0) should be 0.0, got {result_back}" + + def test_negative_zero(self): + """Test nextafter with negative zero""" + neg_zero = QuadPrecision(-0.0) + pos_zero = QuadPrecision(0.0) + one = QuadPrecision(1.0) + neg_one = QuadPrecision(-1.0) + + finfo = np.finfo(QuadPrecDType()) + + # nextafter(-0.0, 1.0) should return smallest positive subnormal + result = np.nextafter(neg_zero, one) + assert result == finfo.smallest_subnormal, \ + f"nextafter(-0.0, 1.0) should be smallest_subnormal, got {result}" + assert result > pos_zero, "Result should be positive" + + # nextafter(+0.0, -1.0) should return smallest negative subnormal + result_neg = np.nextafter(pos_zero, neg_one) + expected_neg_subnormal = -finfo.smallest_subnormal + assert result_neg == expected_neg_subnormal, \ + f"nextafter(+0.0, -1.0) should be -smallest_subnormal, got {result_neg}" + assert result_neg < pos_zero, "Result should be negative" + + def test_infinity_cases(self): + """Test nextafter with infinity edge cases""" + pos_inf = QuadPrecision(np.inf) + neg_inf = QuadPrecision(-np.inf) + one = QuadPrecision(1.0) + neg_one = QuadPrecision(-1.0) + zero = QuadPrecision(0.0) + + finfo = np.finfo(QuadPrecDType()) + + # nextafter(+inf, finite) should return max finite value + result = np.nextafter(pos_inf, zero) + assert not np.isinf(result), "nextafter(+inf, 0) should be finite" + assert result < pos_inf, "Result should be less than +inf" + assert result == finfo.max, f"nextafter(+inf, 0) should be max, got {result} vs {finfo.max}" + + # nextafter(-inf, finite) should return -max (most negative finite) + result_neg = np.nextafter(neg_inf, zero) + assert not np.isinf(result_neg), "nextafter(-inf, 0) should be finite" + assert result_neg > neg_inf, "Result should be greater than -inf" + assert result_neg == -finfo.max, f"nextafter(-inf, 0) should be -max, got {result_neg}" + + # Verify symmetry: nextafter(result, +inf) should give us +inf back + back_to_inf = np.nextafter(result, pos_inf) + assert back_to_inf == pos_inf, "nextafter(max_finite, +inf) should be +inf" + + # nextafter(+inf, +inf) should return +inf + result_inf = np.nextafter(pos_inf, pos_inf) + assert result_inf == pos_inf, "nextafter(+inf, +inf) should be +inf" + + # nextafter(-inf, -inf) should return -inf + result_neg_inf = np.nextafter(neg_inf, neg_inf) + assert result_neg_inf == neg_inf, "nextafter(-inf, -inf) should be -inf" + + def test_max_to_infinity(self): + """Test nextafter from max finite value to infinity""" + finfo = np.finfo(QuadPrecDType()) + max_val = finfo.max + pos_inf = QuadPrecision(np.inf) + neg_inf = QuadPrecision(-np.inf) + + # nextafter(max_finite, +inf) should return +inf + result = np.nextafter(max_val, pos_inf) + assert np.isinf(result), f"nextafter(max, +inf) should be inf, got {result}" + assert result > max_val, "Result should be greater than max" + assert result == pos_inf, "Result should be +inf" + + # nextafter(-max_finite, -inf) should return -inf + neg_max_val = -max_val + result_neg = np.nextafter(neg_max_val, neg_inf) + assert np.isinf(result_neg), f"nextafter(-max, -inf) should be -inf, got {result_neg}" + assert result_neg < neg_max_val, "Result should be less than -max" + assert result_neg == neg_inf, "Result should be -inf" + + def test_near_max(self): + """Test nextafter near maximum finite value""" + finfo = np.finfo(QuadPrecDType()) + max_val = finfo.max + zero = QuadPrecision(0.0) + pos_inf = QuadPrecision(np.inf) + + # nextafter(max, 0) should return a value less than max + result = np.nextafter(max_val, zero) + assert result < max_val, "nextafter(max, 0) should be less than max" + assert not np.isinf(result), "Result should be finite" + + # The difference should be one ULP at that scale + # Moving back should give us max again + result_back = np.nextafter(result, pos_inf) + assert result_back == max_val, f"Moving back should return max, got {result_back}" + + def test_symmetry(self): + """Test symmetry properties of nextafter""" + values = [0.0, 1.0, -1.0, 1e10, -1e10, 1e-10, -1e-10] + + for val in values: + q_val = QuadPrecision(val) + + # nextafter(x, +direction) then nextafter(result, x) should return x + if not np.isinf(val): + result_up = np.nextafter(q_val, QuadPrecision(np.inf)) + result_back = np.nextafter(result_up, q_val) + assert result_back == q_val, \ + f"Symmetry failed for {val}: nextafter then back should return original" + + # Same for down direction + result_down = np.nextafter(q_val, QuadPrecision(-np.inf)) + result_back_down = np.nextafter(result_down, q_val) + assert result_back_down == q_val, \ + f"Symmetry failed for {val}: nextafter down then back should return original" + + def test_direction(self): + """Test that nextafter moves in the correct direction""" + test_cases = [ + (1.0, 2.0, "greater"), # towards larger + (2.0, 1.0, "less"), # towards smaller + (-1.0, -2.0, "less"), # towards more negative + (-2.0, -1.0, "greater"), # towards less negative + (1.0, np.inf, "greater"), # towards +inf + (1.0, -np.inf, "less"), # towards -inf + ] + + for x, y, expected_dir in test_cases: + q_x = QuadPrecision(x) + q_y = QuadPrecision(y) + result = np.nextafter(q_x, q_y) + + if expected_dir == "greater": + assert result > q_x, f"nextafter({x}, {y}) should be > {x}, got {result}" + elif expected_dir == "less": + assert result < q_x, f"nextafter({x}, {y}) should be < {x}, got {result}" + From e32505328b3f2ce893198d5d73b9e07e60fc1369 Mon Sep 17 00:00:00 2001 From: swayaminsync Date: Mon, 20 Oct 2025 03:00:23 +0530 Subject: [PATCH 2/2] spelling fix --- quaddtype/numpy_quaddtype/src/ops.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index a8c4ef3..2eac92e 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/quaddtype/numpy_quaddtype/src/ops.hpp @@ -1003,8 +1003,8 @@ quad_nextafter(const Sleef_quad *x, const Sleef_quad *y) return result; } - // well I did not read to have those above checks - // but they can be important when settng FPE flag manually + // well I did not need to have those above checks + // but they can be important when setting FPE flag manually return quad_set_words64(hx, lx); }