From 3d79e9c168a5e472c5a34afc5f63b751ab72737d Mon Sep 17 00:00:00 2001 From: swayaminsync Date: Wed, 22 Oct 2025 20:04:44 +0530 Subject: [PATCH 1/7] ldexp impl --- quaddtype/numpy_quaddtype/src/ops.hpp | 54 ++++++ .../numpy_quaddtype/src/umath/binary_ops.cpp | 172 ++++++++++++++++++ quaddtype/release_tracker.md | 2 +- quaddtype/tests/test_quaddtype.py | 142 ++++++++++++++- 4 files changed, 368 insertions(+), 2 deletions(-) diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index 638a5706..718b68a4 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/quaddtype/numpy_quaddtype/src/ops.hpp @@ -1047,6 +1047,60 @@ 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 + + // 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 + + // 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..4c9dba9c 100644 --- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp @@ -282,6 +282,175 @@ 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 +// Special resolve descriptors for ldexp (QuadPrecDType, int32) -> 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: int (no need to incref, it's a builtin dtype) + if (given_descrs[1] == NULL) { + loop_descrs[1] = PyArray_DescrFromType(NPY_INT32); + } else { + Py_INCREF(given_descrs[1]); + loop_descrs[1] = given_descrs[1]; + } + + // 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 NPY_NO_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]; // quad + char *in2_ptr = data[1]; // int32 + char *out_ptr = data[2]; // quad + 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; + int in2; + while (N--) { + memcpy(&in1, in1_ptr, elem_size); + memcpy(&in2, in2_ptr, sizeof(int)); + if (backend == BACKEND_SLEEF) { + out.sleef_value = sleef_op(&in1.sleef_value, &in2); + } else { + out.longdouble_value = longdouble_op(&in1.longdouble_value, &in2); + } + 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]; // quad + char *in2_ptr = data[1]; // int32 + char *out_ptr = data[2]; // quad + 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--) { + if (backend == BACKEND_SLEEF) { + *(Sleef_quad *)out_ptr = sleef_op((Sleef_quad *)in1_ptr, (int *)in2_ptr); + } else { + *(long double *)out_ptr = longdouble_op((long double *)in1_ptr, (int *)in2_ptr); + } + + 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 +635,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..a00a8dd5 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -2529,4 +2529,144 @@ 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), + ]) + 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 From 31ec9d616dce5cf4719aac5c594f9ed0bcdb01c0 Mon Sep 17 00:00:00 2001 From: swayaminsync Date: Wed, 22 Oct 2025 21:06:13 +0530 Subject: [PATCH 2/7] strict int32_t usage --- quaddtype/numpy_quaddtype/src/ops.hpp | 14 +++++++++----- quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp | 10 ++++++---- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index 718b68a4..229dd532 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/quaddtype/numpy_quaddtype/src/ops.hpp @@ -1047,14 +1047,15 @@ 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 *); +// Mixed-type operations (quad, int32) -> quad +typedef Sleef_quad (*ldexp_op_quad_def)(const Sleef_quad *, const int32_t *); +typedef long double (*ldexp_op_longdouble_def)(const long double *, const int32_t *); static inline Sleef_quad -quad_ldexp(const Sleef_quad *x, const int *exp) +quad_ldexp(const Sleef_quad *x, const int32_t *exp) { // ldexp(x, exp) returns x * 2^exp + // SLEEF expects: Sleef_quad, int32_t // NaN input -> NaN output (with sign preserved) if (Sleef_iunordq1(*x, *x)) { @@ -1071,15 +1072,17 @@ quad_ldexp(const Sleef_quad *x, const int *exp) return *x; } + // Pass dereferenced values directly to SLEEF (no conversion needed) Sleef_quad result = Sleef_ldexpq1(*x, *exp); return result; } static inline long double -ld_ldexp(const long double *x, const int *exp) +ld_ldexp(const long double *x, const int32_t *exp) { // ldexp(x, exp) returns x * 2^exp + // stdlib ldexpl expects: long double, int // NaN input -> NaN output if (isnan(*x)) { @@ -1096,6 +1099,7 @@ ld_ldexp(const long double *x, const int *exp) return *x; } + // Cast int32_t to int for stdlib function long double result = ldexpl(*x, *exp); return result; diff --git a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp index 4c9dba9c..d6bd00bd 100644 --- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp @@ -345,10 +345,10 @@ quad_ldexp_strided_loop_unaligned(PyArrayMethod_Context *context, char *const da size_t elem_size = (backend == BACKEND_SLEEF) ? sizeof(Sleef_quad) : sizeof(long double); quad_value in1, out; - int in2; + int32_t in2; while (N--) { memcpy(&in1, in1_ptr, elem_size); - memcpy(&in2, in2_ptr, sizeof(int)); + memcpy(&in2, in2_ptr, sizeof(int32_t)); if (backend == BACKEND_SLEEF) { out.sleef_value = sleef_op(&in1.sleef_value, &in2); } else { @@ -382,10 +382,12 @@ quad_ldexp_strided_loop_aligned(PyArrayMethod_Context *context, char *const data QuadBackendType backend = descr->backend; while (N--) { + int32_t *exp = (int32_t *)in2_ptr; + if (backend == BACKEND_SLEEF) { - *(Sleef_quad *)out_ptr = sleef_op((Sleef_quad *)in1_ptr, (int *)in2_ptr); + *(Sleef_quad *)out_ptr = sleef_op((Sleef_quad *)in1_ptr, exp); } else { - *(long double *)out_ptr = longdouble_op((long double *)in1_ptr, (int *)in2_ptr); + *(long double *)out_ptr = longdouble_op((long double *)in1_ptr, exp); } in1_ptr += in1_stride; From 633df7d38b8be32c658984b1b7fa0806e72b3ed5 Mon Sep 17 00:00:00 2001 From: swayaminsync Date: Thu, 23 Oct 2025 11:34:19 +0530 Subject: [PATCH 3/7] switch to system int --- quaddtype/numpy_quaddtype/src/ops.hpp | 14 ++++++-------- quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp | 8 ++++---- 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index 229dd532..71802f33 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); } @@ -1048,14 +1048,14 @@ quad_spacing(const Sleef_quad *x) } // Mixed-type operations (quad, int32) -> quad -typedef Sleef_quad (*ldexp_op_quad_def)(const Sleef_quad *, const int32_t *); -typedef long double (*ldexp_op_longdouble_def)(const long double *, const int32_t *); +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 int32_t *exp) +quad_ldexp(const Sleef_quad *x, const int *exp) { // ldexp(x, exp) returns x * 2^exp - // SLEEF expects: Sleef_quad, int32_t + // SLEEF expects: Sleef_quad, int // NaN input -> NaN output (with sign preserved) if (Sleef_iunordq1(*x, *x)) { @@ -1072,14 +1072,13 @@ quad_ldexp(const Sleef_quad *x, const int32_t *exp) return *x; } - // Pass dereferenced values directly to SLEEF (no conversion needed) Sleef_quad result = Sleef_ldexpq1(*x, *exp); return result; } static inline long double -ld_ldexp(const long double *x, const int32_t *exp) +ld_ldexp(const long double *x, const int *exp) { // ldexp(x, exp) returns x * 2^exp // stdlib ldexpl expects: long double, int @@ -1099,7 +1098,6 @@ ld_ldexp(const long double *x, const int32_t *exp) return *x; } - // Cast int32_t to int for stdlib function long double result = ldexpl(*x, *exp); return result; diff --git a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp index d6bd00bd..34b66b6e 100644 --- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp @@ -298,7 +298,7 @@ quad_ldexp_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtypes[] // Input 1: int (no need to incref, it's a builtin dtype) if (given_descrs[1] == NULL) { - loop_descrs[1] = PyArray_DescrFromType(NPY_INT32); + loop_descrs[1] = PyArray_DescrFromType(NPY_INT); } else { Py_INCREF(given_descrs[1]); loop_descrs[1] = given_descrs[1]; @@ -345,10 +345,10 @@ quad_ldexp_strided_loop_unaligned(PyArrayMethod_Context *context, char *const da size_t elem_size = (backend == BACKEND_SLEEF) ? sizeof(Sleef_quad) : sizeof(long double); quad_value in1, out; - int32_t in2; + int in2; while (N--) { memcpy(&in1, in1_ptr, elem_size); - memcpy(&in2, in2_ptr, sizeof(int32_t)); + memcpy(&in2, in2_ptr, sizeof(int)); if (backend == BACKEND_SLEEF) { out.sleef_value = sleef_op(&in1.sleef_value, &in2); } else { @@ -382,7 +382,7 @@ quad_ldexp_strided_loop_aligned(PyArrayMethod_Context *context, char *const data QuadBackendType backend = descr->backend; while (N--) { - int32_t *exp = (int32_t *)in2_ptr; + int *exp = (int *)in2_ptr; if (backend == BACKEND_SLEEF) { *(Sleef_quad *)out_ptr = sleef_op((Sleef_quad *)in1_ptr, exp); From 46878a069be5d1f05390412f7c5fa759acae20d8 Mon Sep 17 00:00:00 2001 From: swayaminsync Date: Thu, 23 Oct 2025 12:10:18 +0530 Subject: [PATCH 4/7] lets see the actual size --- .../numpy_quaddtype/src/umath/binary_ops.cpp | 54 +++++++++++++++++-- 1 file changed, 50 insertions(+), 4 deletions(-) diff --git a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp index 34b66b6e..0bebab86 100644 --- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp @@ -344,11 +344,38 @@ quad_ldexp_strided_loop_unaligned(PyArrayMethod_Context *context, char *const da QuadBackendType backend = descr->backend; size_t elem_size = (backend == BACKEND_SLEEF) ? sizeof(Sleef_quad) : sizeof(long double); + // Get the actual integer descriptor to determine its size + PyArray_Descr *int_descr = context->descriptors[1]; + int int_size = int_descr->elsize; + quad_value in1, out; int in2; while (N--) { memcpy(&in1, in1_ptr, elem_size); - memcpy(&in2, in2_ptr, sizeof(int)); + + // Read the integer value correctly based on its actual size + // to handle endianness properly + if (int_size == sizeof(npy_int32)) { + npy_int32 temp_int; + memcpy(&temp_int, in2_ptr, sizeof(npy_int32)); + in2 = (int)temp_int; + } else if (int_size == sizeof(npy_int64)) { + npy_int64 temp_int; + memcpy(&temp_int, in2_ptr, sizeof(npy_int64)); + in2 = (int)temp_int; + } else if (int_size == sizeof(npy_int16)) { + npy_int16 temp_int; + memcpy(&temp_int, in2_ptr, sizeof(npy_int16)); + in2 = (int)temp_int; + } else if (int_size == sizeof(npy_int8)) { + npy_int8 temp_int; + memcpy(&temp_int, in2_ptr, sizeof(npy_int8)); + in2 = (int)temp_int; + } else { + // Fallback for other sizes + memcpy(&in2, in2_ptr, sizeof(int)); + } + if (backend == BACKEND_SLEEF) { out.sleef_value = sleef_op(&in1.sleef_value, &in2); } else { @@ -381,13 +408,32 @@ quad_ldexp_strided_loop_aligned(PyArrayMethod_Context *context, char *const data QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0]; QuadBackendType backend = descr->backend; + // Get the actual integer descriptor to determine its size + PyArray_Descr *int_descr = context->descriptors[1]; + int int_size = int_descr->elsize; + while (N--) { - int *exp = (int *)in2_ptr; + int exp_value; + + // Read the integer value correctly based on its actual size + // to handle endianness properly + if (int_size == sizeof(npy_int32)) { + exp_value = (int)(*(npy_int32 *)in2_ptr); + } else if (int_size == sizeof(npy_int64)) { + exp_value = (int)(*(npy_int64 *)in2_ptr); + } else if (int_size == sizeof(npy_int16)) { + exp_value = (int)(*(npy_int16 *)in2_ptr); + } else if (int_size == sizeof(npy_int8)) { + exp_value = (int)(*(npy_int8 *)in2_ptr); + } else { + // Fallback: cast directly (may not work on big-endian) + exp_value = *(int *)in2_ptr; + } if (backend == BACKEND_SLEEF) { - *(Sleef_quad *)out_ptr = sleef_op((Sleef_quad *)in1_ptr, exp); + *(Sleef_quad *)out_ptr = sleef_op((Sleef_quad *)in1_ptr, &exp_value); } else { - *(long double *)out_ptr = longdouble_op((long double *)in1_ptr, exp); + *(long double *)out_ptr = longdouble_op((long double *)in1_ptr, &exp_value); } in1_ptr += in1_stride; From 77d3b404fff8c6173fdaa937007477ded1bdc782 Mon Sep 17 00:00:00 2001 From: swayaminsync Date: Thu, 23 Oct 2025 12:49:07 +0530 Subject: [PATCH 5/7] switch to npy_intp --- .../numpy_quaddtype/src/umath/binary_ops.cpp | 83 +++++-------------- 1 file changed, 19 insertions(+), 64 deletions(-) diff --git a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp index 0bebab86..91d8d3c1 100644 --- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp @@ -283,7 +283,7 @@ quad_generic_binop_2out_strided_loop_aligned(PyArrayMethod_Context *context, cha } // todo: I'll preferrable get all this code duplication in templates later -// Special resolve descriptors for ldexp (QuadPrecDType, int32) -> QuadPrecDType +// 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[], @@ -296,13 +296,9 @@ quad_ldexp_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtypes[] Py_INCREF(given_descrs[0]); loop_descrs[0] = given_descrs[0]; - // Input 1: int (no need to incref, it's a builtin dtype) - if (given_descrs[1] == NULL) { - loop_descrs[1] = PyArray_DescrFromType(NPY_INT); - } else { - Py_INCREF(given_descrs[1]); - loop_descrs[1] = given_descrs[1]; - } + // 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) { @@ -322,7 +318,8 @@ quad_ldexp_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtypes[] loop_descrs[2] = given_descrs[2]; } } - return NPY_NO_CASTING; + // Return SAFE_CASTING to allow conversion from other integer types to intp + return NPY_SAFE_CASTING; } // Strided loop for ldexp (unaligned) @@ -333,9 +330,9 @@ quad_ldexp_strided_loop_unaligned(PyArrayMethod_Context *context, char *const da NpyAuxData *auxdata) { npy_intp N = dimensions[0]; - char *in1_ptr = data[0]; // quad - char *in2_ptr = data[1]; // int32 - char *out_ptr = data[2]; // quad + 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]; @@ -344,42 +341,18 @@ quad_ldexp_strided_loop_unaligned(PyArrayMethod_Context *context, char *const da QuadBackendType backend = descr->backend; size_t elem_size = (backend == BACKEND_SLEEF) ? sizeof(Sleef_quad) : sizeof(long double); - // Get the actual integer descriptor to determine its size - PyArray_Descr *int_descr = context->descriptors[1]; - int int_size = int_descr->elsize; - quad_value in1, out; - int in2; + 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)); - // Read the integer value correctly based on its actual size - // to handle endianness properly - if (int_size == sizeof(npy_int32)) { - npy_int32 temp_int; - memcpy(&temp_int, in2_ptr, sizeof(npy_int32)); - in2 = (int)temp_int; - } else if (int_size == sizeof(npy_int64)) { - npy_int64 temp_int; - memcpy(&temp_int, in2_ptr, sizeof(npy_int64)); - in2 = (int)temp_int; - } else if (int_size == sizeof(npy_int16)) { - npy_int16 temp_int; - memcpy(&temp_int, in2_ptr, sizeof(npy_int16)); - in2 = (int)temp_int; - } else if (int_size == sizeof(npy_int8)) { - npy_int8 temp_int; - memcpy(&temp_int, in2_ptr, sizeof(npy_int8)); - in2 = (int)temp_int; - } else { - // Fallback for other sizes - memcpy(&in2, in2_ptr, sizeof(int)); - } + int exp_value = (int)in2_intp; if (backend == BACKEND_SLEEF) { - out.sleef_value = sleef_op(&in1.sleef_value, &in2); + out.sleef_value = sleef_op(&in1.sleef_value, &exp_value); } else { - out.longdouble_value = longdouble_op(&in1.longdouble_value, &in2); + out.longdouble_value = longdouble_op(&in1.longdouble_value, &exp_value); } memcpy(out_ptr, &out, elem_size); @@ -398,9 +371,9 @@ quad_ldexp_strided_loop_aligned(PyArrayMethod_Context *context, char *const data NpyAuxData *auxdata) { npy_intp N = dimensions[0]; - char *in1_ptr = data[0]; // quad - char *in2_ptr = data[1]; // int32 - char *out_ptr = data[2]; // quad + 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]; @@ -408,27 +381,9 @@ quad_ldexp_strided_loop_aligned(PyArrayMethod_Context *context, char *const data QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0]; QuadBackendType backend = descr->backend; - // Get the actual integer descriptor to determine its size - PyArray_Descr *int_descr = context->descriptors[1]; - int int_size = int_descr->elsize; - while (N--) { - int exp_value; - - // Read the integer value correctly based on its actual size - // to handle endianness properly - if (int_size == sizeof(npy_int32)) { - exp_value = (int)(*(npy_int32 *)in2_ptr); - } else if (int_size == sizeof(npy_int64)) { - exp_value = (int)(*(npy_int64 *)in2_ptr); - } else if (int_size == sizeof(npy_int16)) { - exp_value = (int)(*(npy_int16 *)in2_ptr); - } else if (int_size == sizeof(npy_int8)) { - exp_value = (int)(*(npy_int8 *)in2_ptr); - } else { - // Fallback: cast directly (may not work on big-endian) - exp_value = *(int *)in2_ptr; - } + 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); From d7ff2a9a15657212a70247ce20c1aaf76a7ddc7a Mon Sep 17 00:00:00 2001 From: swayaminsync Date: Thu, 23 Oct 2025 13:10:38 +0530 Subject: [PATCH 6/7] suggested pair --- quaddtype/tests/test_quaddtype.py | 1 + 1 file changed, 1 insertion(+) diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index a00a8dd5..0fc04eab 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -2583,6 +2583,7 @@ def test_ldexp_basic(self, x_val, exp_val): ("-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)""" From 3e81b24908cf1d27dcea4e2306fcaefc1e5888df Mon Sep 17 00:00:00 2001 From: Swayam Date: Thu, 23 Oct 2025 16:43:46 +0530 Subject: [PATCH 7/7] Update quaddtype/numpy_quaddtype/src/ops.hpp Co-authored-by: Juniper Tyree <50025784+juntyr@users.noreply.github.com> --- quaddtype/numpy_quaddtype/src/ops.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index 71802f33..26c67acf 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/quaddtype/numpy_quaddtype/src/ops.hpp @@ -1047,7 +1047,7 @@ quad_spacing(const Sleef_quad *x) return result; } -// Mixed-type operations (quad, int32) -> quad +// 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 *);