diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index 638a5706..26c67acf 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/quaddtype/numpy_quaddtype/src/ops.hpp @@ -28,7 +28,7 @@ quad_positive(const Sleef_quad *op) static inline Sleef_quad quad_sign(const Sleef_quad *op) { - int32_t sign = Sleef_icmpq1(*op, QUAD_ZERO); + int sign = Sleef_icmpq1(*op, QUAD_ZERO); // sign(x=NaN) = x; otherwise sign(x) in { -1.0; 0.0; +1.0 } return Sleef_iunordq1(*op, *op) ? *op : Sleef_cast_from_int64q1(sign); } @@ -1047,6 +1047,62 @@ quad_spacing(const Sleef_quad *x) return result; } +// Mixed-type operations (quad, int) -> quad +typedef Sleef_quad (*ldexp_op_quad_def)(const Sleef_quad *, const int *); +typedef long double (*ldexp_op_longdouble_def)(const long double *, const int *); + +static inline Sleef_quad +quad_ldexp(const Sleef_quad *x, const int *exp) +{ + // ldexp(x, exp) returns x * 2^exp + // SLEEF expects: Sleef_quad, int + + // NaN input -> NaN output (with sign preserved) + if (Sleef_iunordq1(*x, *x)) { + return *x; + } + + // ±0 * 2^exp = ±0 (preserves sign of zero) + if (Sleef_icmpeqq1(*x, QUAD_ZERO)) { + return *x; + } + + // ±inf * 2^exp = ±inf (preserves sign of infinity) + if (quad_isinf(x)) { + return *x; + } + + Sleef_quad result = Sleef_ldexpq1(*x, *exp); + + return result; +} + +static inline long double +ld_ldexp(const long double *x, const int *exp) +{ + // ldexp(x, exp) returns x * 2^exp + // stdlib ldexpl expects: long double, int + + // NaN input -> NaN output + if (isnan(*x)) { + return *x; + } + + // ±0 * 2^exp = ±0 (preserves sign of zero) + if (*x == 0.0L) { + return *x; + } + + // ±inf * 2^exp = ±inf (preserves sign of infinity) + if (isinf(*x)) { + return *x; + } + + long double result = ldexpl(*x, *exp); + + 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) diff --git a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp index 562f277b..91d8d3c1 100644 --- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp @@ -282,6 +282,178 @@ quad_generic_binop_2out_strided_loop_aligned(PyArrayMethod_Context *context, cha return 0; } +// todo: I'll preferrable get all this code duplication in templates later +// resolve descriptors for ldexp (QuadPrecDType, int) -> QuadPrecDType +static NPY_CASTING +quad_ldexp_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtypes[], + PyArray_Descr *const given_descrs[], + PyArray_Descr *loop_descrs[], npy_intp *NPY_UNUSED(view_offset)) +{ + QuadPrecDTypeObject *descr_in1 = (QuadPrecDTypeObject *)given_descrs[0]; + QuadBackendType target_backend = descr_in1->backend; + + // Input 0: QuadPrecDType + Py_INCREF(given_descrs[0]); + loop_descrs[0] = given_descrs[0]; + + // Input 1: Use NPY_INTP (int64 on 64-bit, int32 on 32-bit) to match platform integer size + // This ensures we can handle the full range of PyArray_PyLongDType without data loss + loop_descrs[1] = PyArray_DescrFromType(NPY_INTP); + + // Output: QuadPrecDType with same backend as input + if (given_descrs[2] == NULL) { + loop_descrs[2] = (PyArray_Descr *)new_quaddtype_instance(target_backend); + if (!loop_descrs[2]) { + return (NPY_CASTING)-1; + } + } else { + QuadPrecDTypeObject *descr_out = (QuadPrecDTypeObject *)given_descrs[2]; + if (descr_out->backend != target_backend) { + loop_descrs[2] = (PyArray_Descr *)new_quaddtype_instance(target_backend); + if (!loop_descrs[2]) { + return (NPY_CASTING)-1; + } + } else { + Py_INCREF(given_descrs[2]); + loop_descrs[2] = given_descrs[2]; + } + } + // Return SAFE_CASTING to allow conversion from other integer types to intp + return NPY_SAFE_CASTING; +} + +// Strided loop for ldexp (unaligned) +template +int +quad_ldexp_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 *in1_ptr = data[0]; + char *in2_ptr = data[1]; + char *out_ptr = data[2]; + npy_intp in1_stride = strides[0]; + npy_intp in2_stride = strides[1]; + npy_intp out_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 in1, out; + npy_intp in2_intp; // Platform-native integer (int64 on 64-bit, int32 on 32-bit) + while (N--) { + memcpy(&in1, in1_ptr, elem_size); + memcpy(&in2_intp, in2_ptr, sizeof(npy_intp)); + + int exp_value = (int)in2_intp; + + if (backend == BACKEND_SLEEF) { + out.sleef_value = sleef_op(&in1.sleef_value, &exp_value); + } else { + out.longdouble_value = longdouble_op(&in1.longdouble_value, &exp_value); + } + memcpy(out_ptr, &out, elem_size); + + in1_ptr += in1_stride; + in2_ptr += in2_stride; + out_ptr += out_stride; + } + return 0; +} + +// Strided loop for ldexp (aligned) +template +int +quad_ldexp_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 *in1_ptr = data[0]; + char *in2_ptr = data[1]; + char *out_ptr = data[2]; + npy_intp in1_stride = strides[0]; + npy_intp in2_stride = strides[1]; + npy_intp out_stride = strides[2]; + + QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0]; + QuadBackendType backend = descr->backend; + + while (N--) { + npy_intp exp_intp = *(npy_intp *)in2_ptr; + int exp_value = (int)exp_intp; + + if (backend == BACKEND_SLEEF) { + *(Sleef_quad *)out_ptr = sleef_op((Sleef_quad *)in1_ptr, &exp_value); + } else { + *(long double *)out_ptr = longdouble_op((long double *)in1_ptr, &exp_value); + } + + in1_ptr += in1_stride; + in2_ptr += in2_stride; + out_ptr += out_stride; + } + return 0; +} + +// Create ldexp ufunc +template +int +create_quad_ldexp_ufunc(PyObject *numpy, const char *ufunc_name) +{ + PyObject *ufunc = PyObject_GetAttrString(numpy, ufunc_name); + if (ufunc == NULL) { + return -1; + } + + PyArray_DTypeMeta *dtypes[3] = {&QuadPrecDType, &PyArray_PyLongDType, &QuadPrecDType}; + + PyType_Slot slots[] = { + {NPY_METH_resolve_descriptors, (void *)&quad_ldexp_resolve_descriptors}, + {NPY_METH_strided_loop, + (void *)&quad_ldexp_strided_loop_aligned}, + {NPY_METH_unaligned_strided_loop, + (void *)&quad_ldexp_strided_loop_unaligned}, + {0, NULL}}; + + PyArrayMethod_Spec Spec = { + .name = "quad_ldexp", + .nin = 2, + .nout = 1, + .casting = NPY_NO_CASTING, + .flags = (NPY_ARRAYMETHOD_FLAGS)(NPY_METH_SUPPORTS_UNALIGNED | NPY_METH_IS_REORDERABLE), + .dtypes = dtypes, + .slots = slots, + }; + + if (PyUFunc_AddLoopFromSpec(ufunc, &Spec) < 0) { + return -1; + } + + PyObject *promoter_capsule = + PyCapsule_New((void *)&quad_ufunc_promoter, "numpy._ufunc_promoter", NULL); + if (promoter_capsule == NULL) { + return -1; + } + + PyObject *DTypes = PyTuple_Pack(3, &PyArrayDescr_Type, &PyArray_PyLongDType, &PyArrayDescr_Type); + if (DTypes == 0) { + Py_DECREF(promoter_capsule); + return -1; + } + + if (PyUFunc_AddPromoter(ufunc, DTypes, promoter_capsule) < 0) { + Py_DECREF(promoter_capsule); + Py_DECREF(DTypes); + return -1; + } + Py_DECREF(promoter_capsule); + Py_DECREF(DTypes); + return 0; +} + // Create binary ufunc with 2 outputs (generic for divmod, modf, frexp, etc.) template int @@ -466,5 +638,8 @@ init_quad_binary_ops(PyObject *numpy) if (create_quad_binary_2out_ufunc(numpy, "divmod") < 0) { return -1; } + if (create_quad_ldexp_ufunc(numpy, "ldexp") < 0) { + return -1; + } return 0; } \ No newline at end of file diff --git a/quaddtype/release_tracker.md b/quaddtype/release_tracker.md index e5be6e88..5fd413db 100644 --- a/quaddtype/release_tracker.md +++ b/quaddtype/release_tracker.md @@ -80,7 +80,7 @@ | nextafter | ✅ | ✅ | | spacing | ✅ | ✅ | | modf | ✅ | ✅ | -| ldexp | | | +| ldexp | ✅ | ✅ | | frexp | | | | floor | ✅ | ✅ | | ceil | ✅ | ✅ | diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index 5662a643..0fc04eab 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -2529,4 +2529,145 @@ def test_modf(self, x): np.testing.assert_allclose( reconstructed, quad_x, rtol=1e-12, atol=1e-15, err_msg=f"Reconstruction failed for modf({x}): {quad_int} + {quad_frac} != {quad_x}" - ) \ No newline at end of file + ) + + +class TestLdexp: + """Tests for ldexp function (x * 2**exp)""" + + @pytest.mark.parametrize("x_val,exp_val", [ + ("1.0", 0), + ("1.0", 1), + ("1.0", 2), + ("1.0", 10), + ("1.0", -1), + ("1.0", -2), + ("1.0", -10), + ("2.5", 3), + ("0.5", 5), + ("-3.0", 4), + ("-0.75", -2), + ("1.5", 100), + ("1.5", -100), + ]) + def test_ldexp_basic(self, x_val, exp_val): + """Test ldexp with basic values""" + quad_x = QuadPrecision(x_val) + float_x = np.float64(x_val) + + quad_result = np.ldexp(quad_x, exp_val) + float_result = np.ldexp(float_x, exp_val) + + assert isinstance(quad_result, QuadPrecision), f"Result should be QuadPrecision for ldexp({x_val}, {exp_val})" + + np.testing.assert_allclose( + quad_result, float_result, rtol=1e-12, atol=1e-15, + err_msg=f"ldexp({x_val}, {exp_val}) mismatch" + ) + + # Verify against direct calculation for finite results + if np.isfinite(float_result): + expected = float(x_val) * (2 ** exp_val) + np.testing.assert_allclose( + quad_result, expected, rtol=1e-10, + err_msg=f"ldexp({x_val}, {exp_val}) doesn't match x * 2^exp" + ) + + @pytest.mark.parametrize("x_val,exp_val", [ + ("0.0", 0), + ("0.0", 1), + ("0.0", -1), + ("0.0", 100), + ("0.0", -100), + ("-0.0", 0), + ("-0.0", 1), + ("-0.0", -1), + ("-0.0", 100), + ("-0.0", -100), + ]) + def test_ldexp_zero(self, x_val, exp_val): + """Test ldexp with zero values (should preserve sign)""" + quad_x = QuadPrecision(x_val) + float_x = np.float64(x_val) + + quad_result = np.ldexp(quad_x, exp_val) + float_result = np.ldexp(float_x, exp_val) + + # Zero * 2^exp = zero (with sign preserved) + assert quad_result == 0.0, f"ldexp({x_val}, {exp_val}) should be zero" + assert np.signbit(quad_result) == np.signbit(float_result), \ + f"Sign mismatch for ldexp({x_val}, {exp_val})" + + @pytest.mark.parametrize("x_val,exp_val", [ + ("inf", 0), + ("inf", 1), + ("inf", -1), + ("inf", 100), + ("-inf", 0), + ("-inf", 1), + ("-inf", -1), + ("-inf", 100), + ]) + def test_ldexp_inf(self, x_val, exp_val): + """Test ldexp with infinity (should preserve infinity and sign)""" + quad_x = QuadPrecision(x_val) + float_x = np.float64(x_val) + + quad_result = np.ldexp(quad_x, exp_val) + float_result = np.ldexp(float_x, exp_val) + + assert np.isinf(quad_result), f"ldexp({x_val}, {exp_val}) should be infinity" + assert np.signbit(quad_result) == np.signbit(float_result), \ + f"Sign mismatch for ldexp({x_val}, {exp_val})" + + @pytest.mark.parametrize("x_val,exp_val", [ + ("nan", 0), + ("nan", 1), + ("nan", -1), + ("nan", 100), + ("-nan", 0), + ]) + def test_ldexp_nan(self, x_val, exp_val): + """Test ldexp with NaN (should return NaN)""" + quad_x = QuadPrecision(x_val) + + quad_result = np.ldexp(quad_x, exp_val) + + assert np.isnan(quad_result), f"ldexp({x_val}, {exp_val}) should be NaN" + + @pytest.mark.parametrize("x_val,exp_val", [ + ("1.5", 16384), # Large positive exponent (likely overflow) + ("2.0", 20000), + ]) + def test_ldexp_overflow(self, x_val, exp_val): + """Test ldexp with overflow to infinity""" + quad_x = QuadPrecision(x_val) + float_x = np.float64(x_val) + + quad_result = np.ldexp(quad_x, exp_val) + float_result = np.ldexp(float_x, exp_val) + + # Both should overflow to infinity + assert np.isinf(quad_result), f"ldexp({x_val}, {exp_val}) should overflow to infinity" + assert np.isinf(float_result), f"numpy ldexp({x_val}, {exp_val}) should overflow to infinity" + assert np.signbit(quad_result) == np.signbit(float_result), \ + f"Sign mismatch for overflow ldexp({x_val}, {exp_val})" + + @pytest.mark.parametrize("x_val,exp_val", [ + ("1.5", -16500), # Large negative exponent (likely underflow) + ("2.0", -20000), + ]) + def test_ldexp_underflow(self, x_val, exp_val): + """Test ldexp with underflow to zero""" + quad_x = QuadPrecision(x_val) + float_x = np.float64(x_val) + + quad_result = np.ldexp(quad_x, exp_val) + float_result = np.ldexp(float_x, exp_val) + + # Both should underflow to zero + assert quad_result == 0.0, f"ldexp({x_val}, {exp_val}) should underflow to zero" + 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