Merge branch 'main' into frexp

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);
 }
@@ -1059,6 +1059,7 @@ 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)) {
@@ -1084,6 +1085,7 @@ 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)) {

quaddtype/numpy_quaddtype/src/scalar.c

@@ -14,6 +14,13 @@
 #include "scalar.h"
 #include "scalar_ops.h"
 #include "dragon4.h"
+#include "dtype.h"
+
+// For IEEE 754 binary128 (quad precision), we need 36 decimal digits 
+// to guarantee round-trip conversion (string -> parse -> equals original value)
+// Formula: ceil(1 + MANT_DIG * log10(2)) = ceil(1 + 113 * 0.30103) = 36
+// src: https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format
+#define SLEEF_QUAD_DECIMAL_DIG 36
 
 
 QuadPrecisionObject *
@@ -36,7 +43,77 @@ QuadPrecision_raw_new(QuadBackendType backend)
 
 QuadPrecisionObject *
 QuadPrecision_from_object(PyObject *value, QuadBackendType backend)
-{
+{   
+    // Handle numpy scalars (np.int32, np.float32, etc.) before arrays
+    // We need to check this before PySequence_Check because some numpy scalars are sequences
+    if (PyArray_CheckScalar(value)) {
+        QuadPrecisionObject *self = QuadPrecision_raw_new(backend);
+        if (!self)
+            return NULL;
+        
+        // Try as floating point first
+        if (PyArray_IsScalar(value, Floating)) {
+            PyObject *py_float = PyNumber_Float(value);
+            if (py_float == NULL) {
+                Py_DECREF(self);
+                return NULL;
+            }
+            double dval = PyFloat_AsDouble(py_float);
+            Py_DECREF(py_float);
+            
+            if (backend == BACKEND_SLEEF) {
+                self->value.sleef_value = Sleef_cast_from_doubleq1(dval);
+            }
+            else {
+                self->value.longdouble_value = (long double)dval;
+            }
+            return self;
+        }
+        // Try as integer
+        else if (PyArray_IsScalar(value, Integer)) {
+            PyObject *py_int = PyNumber_Long(value);
+            if (py_int == NULL) {
+                Py_DECREF(self);
+                return NULL;
+            }
+            long long lval = PyLong_AsLongLong(py_int);
+            Py_DECREF(py_int);
+            
+            if (backend == BACKEND_SLEEF) {
+                self->value.sleef_value = Sleef_cast_from_int64q1(lval);
+            }
+            else {
+                self->value.longdouble_value = (long double)lval;
+            }
+            return self;
+        }
+        // For other scalar types, fall through to error handling
+        Py_DECREF(self);
+    }
+    
+    // this checks arrays and sequences (array, tuple)
+    // rejects strings; they're parsed below
+    if (PyArray_Check(value) || (PySequence_Check(value) && !PyUnicode_Check(value) && !PyBytes_Check(value))) 
+    {
+        QuadPrecDTypeObject *dtype_descr = new_quaddtype_instance(backend);
+        if (dtype_descr == NULL) {
+            return NULL;
+        }
+        
+        // steals reference to the descriptor
+        PyObject *result = PyArray_FromAny(
+            value,
+            (PyArray_Descr *)dtype_descr,
+            0,
+            0,
+            NPY_ARRAY_ENSUREARRAY, // this should handle the casting if possible
+            NULL
+        );
+        
+        // PyArray_FromAny steals the reference to dtype_descr, so no need to DECREF
+        return (QuadPrecisionObject *)result;
+    }
+
     QuadPrecisionObject *self = QuadPrecision_raw_new(backend);
     if (!self)
         return NULL;
@@ -99,21 +176,21 @@ QuadPrecision_from_object(PyObject *value, QuadBackendType backend)
             const char *type_cstr = PyUnicode_AsUTF8(type_str);
             if (type_cstr != NULL) {
                 PyErr_Format(PyExc_TypeError,
-                             "QuadPrecision value must be a quad, float, int or string, but got %s "
+                             "QuadPrecision value must be a quad, float, int, string, array or sequence, but got %s "
                              "instead",
                              type_cstr);
             }
             else {
                 PyErr_SetString(
                         PyExc_TypeError,
-                        "QuadPrecision value must be a quad, float, int or string, but got an "
+                        "QuadPrecision value must be a quad, float, int, string, array or sequence, but got an "
                         "unknown type instead");
             }
             Py_DECREF(type_str);
         }
         else {
             PyErr_SetString(PyExc_TypeError,
-                            "QuadPrecision value must be a quad, float, int or string, but got an "
+                            "QuadPrecision value must be a quad, float, int, string, array or sequence, but got an "
                             "unknown type instead");
         }
         Py_DECREF(self);
@@ -152,7 +229,7 @@ QuadPrecision_str_dragon4(QuadPrecisionObject *self)
     Dragon4_Options opt = {.scientific = 0,
                            .digit_mode = DigitMode_Unique,
                            .cutoff_mode = CutoffMode_TotalLength,
-                           .precision = SLEEF_QUAD_DIG,
+                           .precision = SLEEF_QUAD_DECIMAL_DIG,
                            .sign = 1,
                            .trim_mode = TrimMode_LeaveOneZero,
                            .digits_left = 1,
@@ -203,7 +280,7 @@ QuadPrecision_repr_dragon4(QuadPrecisionObject *self)
     Dragon4_Options opt = {.scientific = 1,
                            .digit_mode = DigitMode_Unique,
                            .cutoff_mode = CutoffMode_TotalLength,
-                           .precision = SLEEF_QUAD_DIG,
+                           .precision = SLEEF_QUAD_DECIMAL_DIG,
                            .sign = 1,
                            .trim_mode = TrimMode_LeaveOneZero,
                            .digits_left = 1,

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_INT32);
-    } 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];
@@ -345,14 +342,17 @@ 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;
+    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, in2_ptr, sizeof(int));
+        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, &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);
 
@@ -371,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];
@@ -382,10 +382,13 @@ quad_ldexp_strided_loop_aligned(PyArrayMethod_Context *context, char *const data
     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, (int *)in2_ptr);
+            *(Sleef_quad *)out_ptr = sleef_op((Sleef_quad *)in1_ptr, &exp_value);
         } 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_value);
         }
 
         in1_ptr += in1_stride;