From ee1e16fbb103766bb69779eb90f611671053af2f Mon Sep 17 00:00:00 2001 From: swayaminsync Date: Wed, 22 Oct 2025 20:41:59 +0530 Subject: [PATCH 1/2] frexp impl --- quaddtype/numpy_quaddtype/src/ops.hpp | 65 ++++++ .../numpy_quaddtype/src/umath/binary_ops.cpp | 2 +- .../numpy_quaddtype/src/umath/unary_ops.cpp | 164 ++++++++++++++ quaddtype/release_tracker.md | 2 +- quaddtype/tests/test_quaddtype.py | 202 +++++++++++++++++- 5 files changed, 432 insertions(+), 3 deletions(-) diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index 718b68a4..fd54e99b 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/quaddtype/numpy_quaddtype/src/ops.hpp @@ -1051,6 +1051,10 @@ quad_spacing(const Sleef_quad *x) typedef Sleef_quad (*ldexp_op_quad_def)(const Sleef_quad *, const int *); typedef long double (*ldexp_op_longdouble_def)(const long double *, const int *); +// Frexp operations: quad -> (quad mantissa, int exponent) +typedef Sleef_quad (*frexp_op_quad_def)(const Sleef_quad *, int *); +typedef long double (*frexp_op_longdouble_def)(const long double *, int *); + static inline Sleef_quad quad_ldexp(const Sleef_quad *x, const int *exp) { @@ -1101,6 +1105,67 @@ ld_ldexp(const long double *x, const int *exp) return result; } +static inline Sleef_quad +quad_frexp(const Sleef_quad *x, int *exp) +{ + // frexp(x) returns mantissa m and exponent e such that x = m * 2^e + // where 0.5 <= |m| < 1.0 + // NumPy's documentation says "between -1 and 1" but actual behavior is: + // - Positive x: mantissa in [0.5, 1.0) + // - Negative x: mantissa in (-1.0, -0.5] + // This matches SLEEF's Sleef_frexpq1 behavior exactly. + + // NaN input -> NaN output with exponent 0 + if (Sleef_iunordq1(*x, *x)) { + *exp = 0; + return *x; + } + + // ±0 -> mantissa ±0 with exponent 0 (preserves sign of zero) + if (Sleef_icmpeqq1(*x, QUAD_ZERO)) { + *exp = 0; + return *x; + } + + // ±inf -> mantissa ±inf with exponent 0 (preserves sign of infinity) + if (quad_isinf(x)) { + *exp = 0; + return *x; + } + + Sleef_quad mantissa = Sleef_frexpq1(*x, exp); + + return mantissa; +} + +static inline long double +ld_frexp(const long double *x, int *exp) +{ + // frexp(x) returns mantissa m and exponent e such that x = m * 2^e + + // NaN input -> NaN output with exponent 0 + if (isnan(*x)) { + *exp = 0; + return *x; + } + + // ±0 -> mantissa ±0 with exponent 0 (preserves sign of zero) + if (*x == 0.0L) { + *exp = 0; + return *x; + } + + // ±inf -> mantissa ±inf with exponent 0 (preserves sign of infinity) + if (isinf(*x)) { + *exp = 0; + return *x; + } + + long double mantissa = frexpl(*x, exp); + + return mantissa; +} + // 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) diff --git a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp index 4c9dba9c..97cc42ac 100644 --- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp @@ -451,7 +451,7 @@ create_quad_ldexp_ufunc(PyObject *numpy, const char *ufunc_name) return 0; } -// Create binary ufunc with 2 outputs (generic for divmod, modf, frexp, etc.) +// Create binary ufunc with 2 outputs template int create_quad_binary_2out_ufunc(PyObject *numpy, const char *ufunc_name) diff --git a/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp index 32cae501..d30483f3 100644 --- a/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp @@ -408,6 +408,165 @@ create_quad_unary_2out_ufunc(PyObject *numpy, const char *ufunc_name) return 0; } +// Frexp-specific resolver: QuadPrecDType -> (QuadPrecDType, int32) +static NPY_CASTING +quad_frexp_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtypes[], + PyArray_Descr *const given_descrs[], PyArray_Descr *loop_descrs[], + npy_intp *NPY_UNUSED(view_offset)) +{ + // Input descriptor + Py_INCREF(given_descrs[0]); + loop_descrs[0] = given_descrs[0]; + + // Output 1: mantissa (same dtype as input) + if (given_descrs[1] == NULL) { + Py_INCREF(given_descrs[0]); + loop_descrs[1] = given_descrs[0]; + } + else { + Py_INCREF(given_descrs[1]); + loop_descrs[1] = given_descrs[1]; + } + + // Output 2: exponent (int32) + if (given_descrs[2] == NULL) { + loop_descrs[2] = PyArray_DescrFromType(NPY_INT32); + } + else { + Py_INCREF(given_descrs[2]); + loop_descrs[2] = given_descrs[2]; + } + + QuadPrecDTypeObject *descr_in = (QuadPrecDTypeObject *)given_descrs[0]; + QuadPrecDTypeObject *descr_out1 = (QuadPrecDTypeObject *)loop_descrs[1]; + + if (descr_in->backend != descr_out1->backend) { + return NPY_UNSAFE_CASTING; + } + + return NPY_NO_CASTING; +} + +// Strided loop for frexp (unaligned) +template +int +quad_frexp_strided_loop_unaligned(PyArrayMethod_Context *context, char *const data[], + npy_intp const dimensions[], npy_intp const strides[], + NpyAuxData *auxdata) +{ + npy_intp N = dimensions[0]; + char *in_ptr = data[0]; + char *out_mantissa_ptr = data[1]; + char *out_exp_ptr = data[2]; + npy_intp in_stride = strides[0]; + npy_intp out_mantissa_stride = strides[1]; + npy_intp out_exp_stride = strides[2]; + + QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0]; + QuadBackendType backend = descr->backend; + size_t elem_size = (backend == BACKEND_SLEEF) ? sizeof(Sleef_quad) : sizeof(long double); + + quad_value in, out_mantissa; + int out_exp; + + while (N--) { + memcpy(&in, in_ptr, elem_size); + if (backend == BACKEND_SLEEF) { + out_mantissa.sleef_value = sleef_op(&in.sleef_value, &out_exp); + } + else { + out_mantissa.longdouble_value = longdouble_op(&in.longdouble_value, &out_exp); + } + memcpy(out_mantissa_ptr, &out_mantissa, elem_size); + memcpy(out_exp_ptr, &out_exp, sizeof(int)); + + in_ptr += in_stride; + out_mantissa_ptr += out_mantissa_stride; + out_exp_ptr += out_exp_stride; + } + return 0; +} + +// Strided loop for frexp (aligned) +template +int +quad_frexp_strided_loop_aligned(PyArrayMethod_Context *context, char *const data[], + npy_intp const dimensions[], npy_intp const strides[], + NpyAuxData *auxdata) +{ + npy_intp N = dimensions[0]; + char *in_ptr = data[0]; + char *out_mantissa_ptr = data[1]; + char *out_exp_ptr = data[2]; + npy_intp in_stride = strides[0]; + npy_intp out_mantissa_stride = strides[1]; + npy_intp out_exp_stride = strides[2]; + + QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0]; + QuadBackendType backend = descr->backend; + + int out_exp; + + while (N--) { + if (backend == BACKEND_SLEEF) { + Sleef_quad mantissa = sleef_op((Sleef_quad *)in_ptr, &out_exp); + memcpy(out_mantissa_ptr, &mantissa, sizeof(Sleef_quad)); + } + else { + long double mantissa = longdouble_op((long double *)in_ptr, &out_exp); + memcpy(out_mantissa_ptr, &mantissa, sizeof(long double)); + } + memcpy(out_exp_ptr, &out_exp, sizeof(int)); + + in_ptr += in_stride; + out_mantissa_ptr += out_mantissa_stride; + out_exp_ptr += out_exp_stride; + } + return 0; +} + + +template +int +create_quad_frexp_ufunc(PyObject *numpy, const char *ufunc_name) +{ + PyObject *ufunc = PyObject_GetAttrString(numpy, ufunc_name); + if (ufunc == NULL) { + return -1; + } + + // 1 input (QuadPrecDType), 2 outputs (QuadPrecDType, int32) + PyArray_DTypeMeta *dtypes[3] = { + &QuadPrecDType, + &QuadPrecDType, + &PyArray_Int32DType + }; + + PyType_Slot slots[] = { + {NPY_METH_resolve_descriptors, (void *)&quad_frexp_resolve_descriptors}, + {NPY_METH_strided_loop, + (void *)&quad_frexp_strided_loop_aligned}, + {NPY_METH_unaligned_strided_loop, + (void *)&quad_frexp_strided_loop_unaligned}, + {0, NULL}}; + + PyArrayMethod_Spec Spec = { + .name = "quad_frexp", + .nin = 1, + .nout = 2, + .casting = NPY_NO_CASTING, + .flags = NPY_METH_SUPPORTS_UNALIGNED, + .dtypes = dtypes, + .slots = slots, + }; + + if (PyUFunc_AddLoopFromSpec(ufunc, &Spec) < 0) { + return -1; + } + + return 0; +} + int init_quad_unary_ops(PyObject *numpy) { @@ -540,5 +699,10 @@ init_quad_unary_ops(PyObject *numpy) return -1; } + // Frexp: special case with (QuadPrecDType, int32) outputs + if (create_quad_frexp_ufunc(numpy, "frexp") < 0) { + return -1; + } + return 0; } \ No newline at end of file diff --git a/quaddtype/release_tracker.md b/quaddtype/release_tracker.md index 5fd413db..441f370c 100644 --- a/quaddtype/release_tracker.md +++ b/quaddtype/release_tracker.md @@ -81,7 +81,7 @@ | spacing | ✅ | ✅ | | modf | ✅ | ✅ | | ldexp | ✅ | ✅ | -| frexp | | | +| frexp | ✅ | ✅ | | floor | ✅ | ✅ | | ceil | ✅ | ✅ | | trunc | ✅ | ✅ | diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index a00a8dd5..c6737a62 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -2669,4 +2669,204 @@ def test_ldexp_underflow(self, x_val, exp_val): assert float_result == 0.0, f"numpy ldexp({x_val}, {exp_val}) should underflow to zero" # Sign should be preserved assert np.signbit(quad_result) == np.signbit(float(x_val)), \ - f"Sign should be preserved for underflow ldexp({x_val}, {exp_val})" \ No newline at end of file + f"Sign should be preserved for underflow ldexp({x_val}, {exp_val})" + + +class TestFrexp: + """Tests for frexp function (decompose x into mantissa and exponent)""" + + @pytest.mark.parametrize("x_val", [ + "1.0", + "2.0", + "3.0", + "4.0", + "0.5", + "0.25", + "8.0", + "16.0", + "-1.0", + "-2.0", + "-3.0", + "-4.0", + "-0.5", + "-0.25", + "-8.0", + "-16.0", + "1.5", + "2.5", + "3.5", + "-1.5", + "-2.5", + "0.1", + "0.9", + "1000.0", + "-1000.0", + ]) + def test_frexp_basic(self, x_val): + """Test frexp with basic values - work directly with QuadPrecision""" + quad_x = QuadPrecision(x_val) + + # Get results + quad_m, quad_e = np.frexp(quad_x) + + # Check types + assert isinstance(quad_m, QuadPrecision), f"Mantissa should be QuadPrecision for frexp({x_val})" + assert isinstance(quad_e, (int, np.integer)), f"Exponent should be integer for frexp({x_val})" + + # Verify mantissa is in correct range: 0.5 <= |mantissa| < 1.0 + abs_m = abs(quad_m) + half = QuadPrecision("0.5") + one = QuadPrecision("1.0") + assert abs_m >= half and abs_m < one, \ + f"Mantissa {quad_m} for frexp({x_val}) not in [0.5, 1.0) range" + + # Verify reconstruction: x = mantissa * 2^exponent using ldexp + reconstructed = np.ldexp(quad_m, int(quad_e)) + # Compare directly without float conversion + assert reconstructed == quad_x, \ + f"Reconstruction failed for frexp({x_val}): {reconstructed} != {quad_x}" + + # Compare with NumPy float64 to ensure results are close + float_x = np.float64(x_val) + float_m, float_e = np.frexp(float_x) + + # Mantissa should be close to float64 result (within float64 precision) + np.testing.assert_allclose( + quad_m, float_m, rtol=1e-15, atol=1e-15, + err_msg=f"Mantissa differs from NumPy float64 for frexp({x_val})" + ) + + # Exponent should match exactly for values in float64 range + assert quad_e == float_e, \ + f"Exponent mismatch with NumPy for frexp({x_val}): {quad_e} != {float_e}" + + @pytest.mark.parametrize("x_val", [ + "0.0", + "-0.0", + ]) + def test_frexp_zero(self, x_val): + """Test frexp with zero values (should return ±0 mantissa, exponent 0)""" + quad_x = QuadPrecision(x_val) + + quad_m, quad_e = np.frexp(quad_x) + + # Mantissa should be zero with same sign + zero = QuadPrecision("0.0") + assert quad_m == zero, f"Mantissa should be zero for frexp({x_val})" + assert np.signbit(quad_m) == np.signbit(quad_x), \ + f"Sign mismatch for frexp({x_val}) mantissa" + + # Exponent should be 0 + assert quad_e == 0, f"Exponent should be 0 for frexp({x_val})" + + # Compare with NumPy float64 + float_x = np.float64(x_val) + float_m, float_e = np.frexp(float_x) + + # Mantissa should match (both zero with same sign) + np.testing.assert_allclose( + quad_m, float_m, rtol=0, atol=0, + err_msg=f"Mantissa differs from NumPy float64 for frexp({x_val})" + ) + assert np.signbit(quad_m) == np.signbit(float_m), \ + f"Sign mismatch with NumPy for frexp({x_val})" + + # Exponent should match + assert quad_e == float_e, \ + f"Exponent mismatch with NumPy for frexp({x_val}): {quad_e} != {float_e}" + + @pytest.mark.parametrize("x_val", [ + "inf", + "-inf", + ]) + def test_frexp_inf(self, x_val): + """Test frexp with infinity (should return ±inf mantissa, exponent 0)""" + quad_x = QuadPrecision(x_val) + + quad_m, quad_e = np.frexp(quad_x) + + # Mantissa should be infinity with same sign + assert np.isinf(quad_m), f"Mantissa should be infinity for frexp({x_val})" + assert np.signbit(quad_m) == np.signbit(quad_x), \ + f"Sign mismatch for frexp({x_val}) mantissa" + + # Exponent should be 0 + assert quad_e == 0, f"Exponent should be 0 for frexp({x_val})" + + # Compare with NumPy float64 + float_x = np.float64(x_val) + float_m, float_e = np.frexp(float_x) + + # Both should be infinity with same sign + assert np.isinf(float_m), f"NumPy mantissa should also be infinity for frexp({x_val})" + assert np.signbit(quad_m) == np.signbit(float_m), \ + f"Sign mismatch with NumPy for frexp({x_val})" + + # Exponent should match + assert quad_e == float_e, \ + f"Exponent mismatch with NumPy for frexp({x_val}): {quad_e} != {float_e}" + + @pytest.mark.parametrize("x_val", [ + "nan", + "-nan", + ]) + def test_frexp_nan(self, x_val): + """Test frexp with NaN (should return NaN mantissa, exponent 0)""" + quad_x = QuadPrecision(x_val) + + quad_m, quad_e = np.frexp(quad_x) + + # Mantissa should be NaN + assert np.isnan(quad_m), f"Mantissa should be NaN for frexp({x_val})" + + # Exponent should be 0 + assert quad_e == 0, f"Exponent should be 0 for frexp({x_val})" + + # Compare with NumPy float64 + float_x = np.float64(x_val) + float_m, float_e = np.frexp(float_x) + + # Both should be NaN + assert np.isnan(float_m), f"NumPy mantissa should also be NaN for frexp({x_val})" + + # Exponent should match + assert quad_e == float_e, \ + f"Exponent mismatch with NumPy for frexp({x_val}): {quad_e} != {float_e}" + + def test_frexp_very_large(self): + """Test frexp with very large values""" + # Large value that's still finite + quad_x = QuadPrecision("1e100") + + quad_m, quad_e = np.frexp(quad_x) + + # Verify mantissa range + abs_m = abs(quad_m) + half = QuadPrecision("0.5") + one = QuadPrecision("1.0") + assert abs_m >= half and abs_m < one, \ + f"Mantissa {quad_m} for large value not in [0.5, 1.0) range" + + # Verify reconstruction using ldexp (preserves full quad precision) + reconstructed = np.ldexp(quad_m, int(quad_e)) + assert reconstructed == quad_x, \ + f"Reconstruction failed for large value: {reconstructed} != {quad_x}" + + def test_frexp_very_small(self): + """Test frexp with very small positive values""" + # Small positive value + quad_x = QuadPrecision("1e-100") + + quad_m, quad_e = np.frexp(quad_x) + + # Verify mantissa range + abs_m = abs(quad_m) + half = QuadPrecision("0.5") + one = QuadPrecision("1.0") + assert abs_m >= half and abs_m < one, \ + f"Mantissa {quad_m} for small value not in [0.5, 1.0) range" + + # Verify reconstruction using ldexp (preserves full quad precision) + reconstructed = np.ldexp(quad_m, int(quad_e)) + assert reconstructed == quad_x, \ + f"Reconstruction failed for small value: {reconstructed} != {quad_x}" \ No newline at end of file From 6727239952812a9fe07e50f372b647ad23682ada Mon Sep 17 00:00:00 2001 From: SwayamInSync Date: Fri, 24 Oct 2025 10:33:49 +0000 Subject: [PATCH 2/2] special case tests --- quaddtype/tests/test_quaddtype.py | 73 ++++++++++++------------------- 1 file changed, 29 insertions(+), 44 deletions(-) diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index 68f645ad..0ceb56c6 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -3000,60 +3000,45 @@ def test_frexp_zero(self, x_val): @pytest.mark.parametrize("x_val", [ "inf", "-inf", - ]) - def test_frexp_inf(self, x_val): - """Test frexp with infinity (should return ±inf mantissa, exponent 0)""" - quad_x = QuadPrecision(x_val) - - quad_m, quad_e = np.frexp(quad_x) - - # Mantissa should be infinity with same sign - assert np.isinf(quad_m), f"Mantissa should be infinity for frexp({x_val})" - assert np.signbit(quad_m) == np.signbit(quad_x), \ - f"Sign mismatch for frexp({x_val}) mantissa" - - # Exponent should be 0 - assert quad_e == 0, f"Exponent should be 0 for frexp({x_val})" - - # Compare with NumPy float64 - float_x = np.float64(x_val) - float_m, float_e = np.frexp(float_x) - - # Both should be infinity with same sign - assert np.isinf(float_m), f"NumPy mantissa should also be infinity for frexp({x_val})" - assert np.signbit(quad_m) == np.signbit(float_m), \ - f"Sign mismatch with NumPy for frexp({x_val})" - - # Exponent should match - assert quad_e == float_e, \ - f"Exponent mismatch with NumPy for frexp({x_val}): {quad_e} != {float_e}" - - @pytest.mark.parametrize("x_val", [ "nan", "-nan", ]) - def test_frexp_nan(self, x_val): - """Test frexp with NaN (should return NaN mantissa, exponent 0)""" + def test_frexp_special_values(self, x_val): + """Test frexp with special values (inf, nan) + + For these edge cases, the C standard specifies that the exponent value + is unspecified/implementation-defined. We only verify that: + 1. The mantissa matches the expected value (±inf or NaN) + 2. The mantissa behavior matches NumPy's float64 + 3. The exponent is an integer type + + We do NOT compare exponent values as they can differ across platforms + (e.g., Linux returns 0, Windows returns -1). + """ quad_x = QuadPrecision(x_val) - quad_m, quad_e = np.frexp(quad_x) - # Mantissa should be NaN - assert np.isnan(quad_m), f"Mantissa should be NaN for frexp({x_val})" - - # Exponent should be 0 - assert quad_e == 0, f"Exponent should be 0 for frexp({x_val})" - # Compare with NumPy float64 float_x = np.float64(x_val) float_m, float_e = np.frexp(float_x) - # Both should be NaN - assert np.isnan(float_m), f"NumPy mantissa should also be NaN for frexp({x_val})" - - # Exponent should match - assert quad_e == float_e, \ - f"Exponent mismatch with NumPy for frexp({x_val}): {quad_e} != {float_e}" + # Exponent should be an integer type (but we don't check the value) + assert isinstance(quad_e, (int, np.integer)), \ + f"Exponent should be integer type for frexp({x_val})" + + # Check mantissa behavior + if "inf" in x_val: + # Mantissa should be infinity with same sign + assert np.isinf(quad_m), f"Mantissa should be infinity for frexp({x_val})" + assert np.isinf(float_m), f"NumPy mantissa should also be infinity for frexp({x_val})" + assert np.signbit(quad_m) == np.signbit(quad_x), \ + f"Sign mismatch for frexp({x_val}) mantissa" + assert np.signbit(quad_m) == np.signbit(float_m), \ + f"Sign mismatch with NumPy for frexp({x_val})" + else: # nan + # Mantissa should be NaN + assert np.isnan(quad_m), f"Mantissa should be NaN for frexp({x_val})" + assert np.isnan(float_m), f"NumPy mantissa should also be NaN for frexp({x_val})" def test_frexp_very_large(self): """Test frexp with very large values"""