Merge pull request #195 from SwayamInSync/ldexp

Author: SwayamInSync

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)

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 <ldexp_op_quad_def sleef_op, ldexp_op_longdouble_def longdouble_op>
+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 <ldexp_op_quad_def sleef_op, ldexp_op_longdouble_def longdouble_op>
+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 <ldexp_op_quad_def sleef_op, ldexp_op_longdouble_def longdouble_op>
+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<sleef_op, longdouble_op>},
+            {NPY_METH_unaligned_strided_loop,
+             (void *)&quad_ldexp_strided_loop_unaligned<sleef_op, longdouble_op>},
+            {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 <binary_op_2out_quad_def sleef_op, binary_op_2out_longdouble_def longdouble_op>
 int
@@ -466,5 +638,8 @@ init_quad_binary_ops(PyObject *numpy)
     if (create_quad_binary_2out_ufunc<quad_divmod, ld_divmod>(numpy, "divmod") < 0) {
         return -1;
     }
+    if (create_quad_ldexp_ufunc<quad_ldexp, ld_ldexp>(numpy, "ldexp") < 0) {
+        return -1;
+    }
     return 0;
 }

quaddtype/release_tracker.md

@@ -80,7 +80,7 @@
 | nextafter     | ✅    | ✅                                                                   |
 | spacing       | ✅    | ✅                                                                   |
 | modf          | ✅    | ✅                                                                   |
-| ldexp         |       |                                                                      |
+| ldexp         | ✅    | ✅                                                                   |
 | frexp         |       |                                                                      |
 | floor         | ✅    | ✅                                                                   |
 | ceil          | ✅    | ✅                                                                   |