QuadDtype finfo support & Nightly CI

quaddtype/meson.build

@@ -117,6 +117,7 @@ srcs = [
     'numpy_quaddtype/src/umath/promoters.hpp',
     'numpy_quaddtype/src/umath/matmul.h',
     'numpy_quaddtype/src/umath/matmul.cpp',
+    'numpy_quaddtype/src/constants.hpp',
 ]
 
 py.install_sources(

quaddtype/numpy_quaddtype/src/constants.hpp

@@ -0,0 +1,133 @@
+#ifndef QUAD_CONSTANTS_HPP
+#define QUAD_CONSTANTS_HPP
+
+#include <sleef.h>
+#include <sleefquad.h>
+#include <stdint.h>
+#include <string.h>
+
+// Quad precision constants using sleef_q macro
+#define QUAD_PRECISION_ZERO sleef_q(+0x0000000000000LL, 0x0000000000000000ULL, -16383)
+#define QUAD_PRECISION_ONE sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 0)
+#define QUAD_PRECISION_INF sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 16384)
+#define QUAD_PRECISION_NINF sleef_q(-0x1000000000000LL, 0x0000000000000000ULL, 16384)
+#define QUAD_PRECISION_NAN sleef_q(+0x1ffffffffffffLL, 0xffffffffffffffffULL, 16384)
+
+// Additional constants
+#define QUAD_PRECISION_MAX_FINITE SLEEF_QUAD_MAX
+#define QUAD_PRECISION_MIN_FINITE Sleef_negq1(SLEEF_QUAD_MAX)
+#define QUAD_PRECISION_RADIX sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 1)  // 2.0
+
+#ifdef SLEEF_QUAD_C
+static const Sleef_quad SMALLEST_SUBNORMAL_VALUE = SLEEF_QUAD_DENORM_MIN;
+#else
+static const union {
+    struct {
+#if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_BIG_ENDIAN__)
+        uint64_t h, l;
+#else
+        uint64_t l, h;
+#endif
+    } parts;
+    Sleef_quad value;
+} smallest_subnormal_const = {.parts = {
+#if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_BIG_ENDIAN__)
+                                      .h = 0x0000000000000000ULL, .l = 0x0000000000000001ULL
+#else
+                                      .l = 0x0000000000000001ULL, .h = 0x0000000000000000ULL
+#endif
+                              }};
+#define SMALLEST_SUBNORMAL_VALUE (smallest_subnormal_const.value)
+#endif
+
+// Integer constants for finfo
+#define QUAD_NMANT 112           // mantissa bits (excluding implicit bit)
+#define QUAD_MIN_EXP -16382      // minimum exponent for normalized numbers
+#define QUAD_MAX_EXP 16384       // maximum exponent
+#define QUAD_DECIMAL_DIGITS 33   // decimal digits of precision
+
+typedef enum ConstantResultType {
+    CONSTANT_QUAD,      // Sleef_quad value
+    CONSTANT_INT64,     // int64_t value
+    CONSTANT_ERROR      // Error occurred
+} ConstantResultType;
+
+typedef struct ConstantResult {
+    ConstantResultType type;
+    union {
+        Sleef_quad quad_value;
+        int64_t int_value;
+    } data;
+} ConstantResult;
+
+
+static inline ConstantResult get_sleef_constant_by_name(const char* constant_name) {
+    ConstantResult result;
+
+    if (strcmp(constant_name, "pi") == 0) {
+        result.type = CONSTANT_QUAD;
+        result.data.quad_value = SLEEF_M_PIq;
+    }
+    else if (strcmp(constant_name, "e") == 0) {
+        result.type = CONSTANT_QUAD;
+        result.data.quad_value = SLEEF_M_Eq;
+    }
+    else if (strcmp(constant_name, "log2e") == 0) {
+        result.type = CONSTANT_QUAD;
+        result.data.quad_value = SLEEF_M_LOG2Eq;
+    }
+    else if (strcmp(constant_name, "log10e") == 0) {
+        result.type = CONSTANT_QUAD;
+        result.data.quad_value = SLEEF_M_LOG10Eq;
+    }
+    else if (strcmp(constant_name, "ln2") == 0) {
+        result.type = CONSTANT_QUAD;
+        result.data.quad_value = SLEEF_M_LN2q;
+    }
+    else if (strcmp(constant_name, "ln10") == 0) {
+        result.type = CONSTANT_QUAD;
+        result.data.quad_value = SLEEF_M_LN10q;
+    }
+    else if (strcmp(constant_name, "max_value") == 0) {
+        result.type = CONSTANT_QUAD;
+        result.data.quad_value = SLEEF_QUAD_MAX;
+    }
+    else if (strcmp(constant_name, "epsilon") == 0) {
+        result.type = CONSTANT_QUAD;
+        result.data.quad_value = SLEEF_QUAD_EPSILON;
+    }
+    else if (strcmp(constant_name, "smallest_normal") == 0) {
+        result.type = CONSTANT_QUAD;
+        result.data.quad_value = SLEEF_QUAD_MIN;
+    }
+    else if (strcmp(constant_name, "smallest_subnormal") == 0) {
+        result.type = CONSTANT_QUAD;
+        result.data.quad_value = SMALLEST_SUBNORMAL_VALUE;
+    }
+    else if (strcmp(constant_name, "bits") == 0) {
+        result.type = CONSTANT_INT64;
+        result.data.int_value = sizeof(Sleef_quad) * 8;
+    }
+    else if (strcmp(constant_name, "precision") == 0) {
+        result.type = CONSTANT_INT64;
+        // precision = int(-log10(epsilon))
+        result.data.int_value = 
+            Sleef_cast_to_int64q1(Sleef_negq1(Sleef_log10q1_u10(SLEEF_QUAD_EPSILON)));
+    }
+    else if (strcmp(constant_name, "resolution") == 0) {
+        result.type = CONSTANT_QUAD;
+        // precision = int(-log10(epsilon))
+        int64_t precision = 
+            Sleef_cast_to_int64q1(Sleef_negq1(Sleef_log10q1_u10(SLEEF_QUAD_EPSILON)));
+        // resolution = 10 ** (-precision)
+        result.data.quad_value = 
+            Sleef_powq1_u10(Sleef_cast_from_int64q1(10), Sleef_cast_from_int64q1(-precision));
+    }
+    else {
+        result.type = CONSTANT_ERROR;
+    }
+
+    return result;
+}
+
+#endif // QUAD_CONSTANTS_HPP

quaddtype/numpy_quaddtype/src/dtype.c

@@ -16,6 +16,7 @@
 #include "casts.h"
 #include "dtype.h"
 #include "dragon4.h"
+#include "constants.hpp"
 
 static inline int
 quad_load(void *x, char *data_ptr, QuadBackendType backend)
@@ -176,6 +177,90 @@ quadprec_default_descr(PyArray_DTypeMeta *cls)
     return (PyArray_Descr *)temp;
 }
 
+static int
+quadprec_get_constant(PyArray_Descr *descr, int constant_id, void *ptr)
+{
+    QuadPrecDTypeObject *quad_descr = (QuadPrecDTypeObject *)descr;
+
+    if (quad_descr->backend != BACKEND_SLEEF) {
+        /* Long double backend not yet implemented */
+        return 0;
+    }
+
+    Sleef_quad val;
+
+    switch (constant_id) {
+        case NPY_CONSTANT_zero:
+            val = QUAD_PRECISION_ZERO;
+            break;
+
+        case NPY_CONSTANT_one:
+            val = QUAD_PRECISION_ONE;
+            break;
+
+        case NPY_CONSTANT_minimum_finite:
+            val = QUAD_PRECISION_MIN_FINITE;
+            break;
+
+        case NPY_CONSTANT_maximum_finite:
+            val = QUAD_PRECISION_MAX_FINITE;
+            break;
+
+        case NPY_CONSTANT_inf:
+            val = QUAD_PRECISION_INF;
+            break;
+
+        case NPY_CONSTANT_ninf:
+            val = QUAD_PRECISION_NINF;
+            break;
+
+        case NPY_CONSTANT_nan:
+            val = QUAD_PRECISION_NAN;
+            break;
+
+        case NPY_CONSTANT_finfo_radix:
+            val = QUAD_PRECISION_RADIX;
+            break;
+
+        case NPY_CONSTANT_finfo_eps:
+            val = SLEEF_QUAD_EPSILON;
+            break;
+
+        case NPY_CONSTANT_finfo_smallest_normal:
+            val = SLEEF_QUAD_MIN;
+            break;
+
+        case NPY_CONSTANT_finfo_smallest_subnormal:
+            val = SMALLEST_SUBNORMAL_VALUE;
+            break;
+
+        /* Integer constants - these return npy_intp values */
+        case NPY_CONSTANT_finfo_nmant:
+            *(npy_intp *)ptr = QUAD_NMANT;
+            return 1;
+
+        case NPY_CONSTANT_finfo_min_exp:
+            *(npy_intp *)ptr = QUAD_MIN_EXP;
+            return 1;
+
+        case NPY_CONSTANT_finfo_max_exp:
+            *(npy_intp *)ptr = QUAD_MAX_EXP;
+            return 1;
+
+        case NPY_CONSTANT_finfo_decimal_digits:
+            *(npy_intp *)ptr = QUAD_DECIMAL_DIGITS;
+            return 1;
+
+        default:
+            /* Constant not supported */
+            return 0;
+    }
+
+    /* Store the Sleef_quad value to the provided pointer */
+    *(Sleef_quad *)ptr = val;
+    return 1;
+}
+
 static PyType_Slot QuadPrecDType_Slots[] = {
         {NPY_DT_ensure_canonical, &ensure_canonical},
         {NPY_DT_common_instance, &common_instance},
@@ -184,7 +269,7 @@ static PyType_Slot QuadPrecDType_Slots[] = {
         {NPY_DT_setitem, &quadprec_setitem},
         {NPY_DT_getitem, &quadprec_getitem},
         {NPY_DT_default_descr, &quadprec_default_descr},
-        {NPY_DT_PyArray_ArrFuncs_dotfunc, NULL},
+        {NPY_DT_get_constant, &quadprec_get_constant},
         {0, NULL}};
 
 static PyObject *

quaddtype/numpy_quaddtype/src/quaddtype_main.c

@@ -17,6 +17,7 @@
 #include "umath/umath.h"
 #include "quad_common.h"
 #include "quadblas_interface.h"
+#include "constants.hpp"
 #include "float.h"
 
 static PyObject *
@@ -30,28 +31,6 @@ py_is_longdouble_128(PyObject *self, PyObject *args)
     }
 }
 
-#ifdef SLEEF_QUAD_C
-static const Sleef_quad SMALLEST_SUBNORMAL_VALUE = SLEEF_QUAD_DENORM_MIN;
-#else
-static const union {
-    struct {
-#if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_BIG_ENDIAN__)
-        uint64_t h, l;
-#else
-        uint64_t l, h;
-#endif
-    } parts;
-    Sleef_quad value;
-} smallest_subnormal_const = {.parts = {
-#if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_BIG_ENDIAN__)
-                                      .h = 0x0000000000000000ULL, .l = 0x0000000000000001ULL
-#else
-                                      .l = 0x0000000000000001ULL, .h = 0x0000000000000000ULL
-#endif
-                              }};
-#define SMALLEST_SUBNORMAL_VALUE (smallest_subnormal_const.value)
-#endif
-
 static PyObject *
 get_sleef_constant(PyObject *self, PyObject *args)
 {
@@ -60,68 +39,29 @@ get_sleef_constant(PyObject *self, PyObject *args)
         return NULL;
     }
 
-    QuadPrecisionObject *result = QuadPrecision_raw_new(BACKEND_SLEEF);
-    if (result == NULL) {
-        return NULL;
-    }
-
-    if (strcmp(constant_name, "pi") == 0) {
-        result->value.sleef_value = SLEEF_M_PIq;
-    }
-    else if (strcmp(constant_name, "e") == 0) {
-        result->value.sleef_value = SLEEF_M_Eq;
-    }
-    else if (strcmp(constant_name, "log2e") == 0) {
-        result->value.sleef_value = SLEEF_M_LOG2Eq;
-    }
-    else if (strcmp(constant_name, "log10e") == 0) {
-        result->value.sleef_value = SLEEF_M_LOG10Eq;
-    }
-    else if (strcmp(constant_name, "ln2") == 0) {
-        result->value.sleef_value = SLEEF_M_LN2q;
-    }
-    else if (strcmp(constant_name, "ln10") == 0) {
-        result->value.sleef_value = SLEEF_M_LN10q;
-    }
-    else if (strcmp(constant_name, "max_value") == 0) {
-        result->value.sleef_value = SLEEF_QUAD_MAX;
-    }
-    else if (strcmp(constant_name, "epsilon") == 0) {
-        result->value.sleef_value = SLEEF_QUAD_EPSILON;
-    }
-    else if (strcmp(constant_name, "smallest_normal") == 0) {
-        result->value.sleef_value = SLEEF_QUAD_MIN;
-    }
-    else if (strcmp(constant_name, "smallest_subnormal") == 0) {
-        // or just use sleef_q(+0x0000000000000LL, 0x0000000000000001ULL, -16383);
-        result->value.sleef_value = SMALLEST_SUBNORMAL_VALUE;
-    }
-    else if (strcmp(constant_name, "bits") == 0) {
-        Py_DECREF(result);
-        return PyLong_FromLong(sizeof(Sleef_quad) * CHAR_BIT);
-    }
-    else if (strcmp(constant_name, "precision") == 0) {
-        Py_DECREF(result);
-        // precision = int(-log10(epsilon))
-        int64_t precision =
-                Sleef_cast_to_int64q1(Sleef_negq1(Sleef_log10q1_u10(SLEEF_QUAD_EPSILON)));
-        return PyLong_FromLong(precision);
-    }
-    else if (strcmp(constant_name, "resolution") == 0) {
-        // precision = int(-log10(epsilon))
-        int64_t precision =
-                Sleef_cast_to_int64q1(Sleef_negq1(Sleef_log10q1_u10(SLEEF_QUAD_EPSILON)));
-        // resolution = 10 ** (-precision)
-        result->value.sleef_value =
-                Sleef_powq1_u10(Sleef_cast_from_int64q1(10), Sleef_cast_from_int64q1(-precision));
-    }
-    else {
+    ConstantResult const_result = get_sleef_constant_by_name(constant_name);
+
+    if (const_result.type == CONSTANT_ERROR) {
         PyErr_SetString(PyExc_ValueError, "Unknown constant name");
-        Py_DECREF(result);
         return NULL;
     }
-
-    return (PyObject *)result;
+
+    if (const_result.type == CONSTANT_INT64) {
+        return PyLong_FromLongLong(const_result.data.int_value);
+    }
+
+    if (const_result.type == CONSTANT_QUAD) {
+        QuadPrecisionObject *result = QuadPrecision_raw_new(BACKEND_SLEEF);
+        if (result == NULL) {
+            return NULL;
+        }
+        result->value.sleef_value = const_result.data.quad_value;
+        return (PyObject *)result;
+    }
+
+    // Should never reach here
+    PyErr_SetString(PyExc_RuntimeError, "Unexpected constant result type");
+    return NULL;
 }
 
 static PyMethodDef module_methods[] = {

quaddtype/numpy_quaddtype/src/scalar.c

@@ -253,5 +253,6 @@ PyTypeObject QuadPrecision_Type = {
 int
 init_quadprecision_scalar(void)
 {
+    QuadPrecision_Type.tp_base = &PyFloatingArrType_Type;
     return PyType_Ready(&QuadPrecision_Type);
 }

quaddtype/reinstall.sh

@@ -8,7 +8,7 @@ if [ -d "build/" ]; then
     rm -rf subprojects/sleef
 fi
 
-export CFLAGS="-g -O0" 
-export CXXFLAGS="-g -O0"
+# export CFLAGS="-g -O0" 
+# export CXXFLAGS="-g -O0"
 python -m pip uninstall -y numpy_quaddtype
-python -m pip install . -v
+python -m pip install . -v --no-build-isolation 2>&1 | tee build_log.txt