Skip to content

WIP: QuadDtype finfo support #165

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 3 commits into
base: main
Choose a base branch
from
Draft

WIP: QuadDtype finfo support #165

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
fea68bd
small init fixes
SwayamInSync Oct 6, 2025
317a961
small changes to reinstall.sh
SwayamInSync Oct 6, 2025
8ba5167
constant refactor and initial finfo support
SwayamInSync Oct 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions quaddtype/meson.build
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
133 changes: 133 additions & 0 deletions quaddtype/numpy_quaddtype/src/constants.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
87 changes: 86 additions & 1 deletion quaddtype/numpy_quaddtype/src/dtype.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you support unaligned, should use memcpy, unfortunately. (Or we change that we promise to not call it for unaligned stuff, which I can get on-board with.)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually thinking about it, I am not sure if this is true or not :), we needed copyswap anyway due to byte-swapping. I am actually happy to say that the data must be aligned.

return 1;
}

static PyType_Slot QuadPrecDType_Slots[] = {
{NPY_DT_ensure_canonical, &ensure_canonical},
{NPY_DT_common_instance, &common_instance},
Expand All @@ -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 *
Expand Down
102 changes: 21 additions & 81 deletions quaddtype/numpy_quaddtype/src/quaddtype_main.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
#include "umath/umath.h"
#include "quad_common.h"
#include "quadblas_interface.h"
#include "constants.hpp"
#include "float.h"

static PyObject *
Expand All @@ -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)
{
Expand All @@ -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[] = {
Expand Down
1 change: 1 addition & 0 deletions quaddtype/numpy_quaddtype/src/scalar.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -253,5 +253,6 @@ PyTypeObject QuadPrecision_Type = {
int
init_quadprecision_scalar(void)
{
QuadPrecision_Type.tp_base = &PyFloatingArrType_Type;
return PyType_Ready(&QuadPrecision_Type);
}
6 changes: 3 additions & 3 deletions quaddtype/reinstall.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading
Loading