diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 286290c0..d0d32495 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -55,11 +55,21 @@ jobs: working-directory: unytdtype run: | pytest -vvv --color=yes + - name: Install quaddtype dependencies + run: | + sudo apt-get update + sudo apt-get install -y libmpfr-dev libssl-dev libfftw3-dev + - name: Install SLEEF + run: | + git clone https://github.com/shibatch/sleef.git + cd sleef + cmake -S . -B build -DSLEEF_BUILD_QUAD:BOOL=ON -DSLEEF_BUILD_SHARED_LIBS:BOOL=ON -DCMAKE_POSITION_INDEPENDENT_CODE=ON + cmake --build build/ --clean-first -j + sudo cmake --install build --prefix /usr - name: Install quaddtype working-directory: quaddtype run: | - python -m build --no-isolation --wheel -Cbuilddir=build - find ./dist/*.whl | xargs python -m pip install + LDFLAGS="-Wl,-rpath,/usr/lib" python -m pip install . -v --no-build-isolation -Cbuilddir=build -C'compile-args=-v' -Csetup-args="-Dbuildtype=debug" - name: Run quaddtype tests working-directory: quaddtype run: | diff --git a/.gitignore b/.gitignore index 0bc2e657..65005454 100644 --- a/.gitignore +++ b/.gitignore @@ -133,3 +133,4 @@ compile_commands.json .ruff-cache/ .asv +.vscode/ \ No newline at end of file diff --git a/quaddtype/.flake8 b/quaddtype/.flake8 deleted file mode 100644 index e82426e8..00000000 --- a/quaddtype/.flake8 +++ /dev/null @@ -1,6 +0,0 @@ -[flake8] -max-line-length = 100 -ignore = D107,D104,W503 -per-file-ignores = __init__.py:F401,tests/*:D100 -docstring-convention = numpy -docstring_style = numpy diff --git a/quaddtype/README.md b/quaddtype/README.md index f1b9322e..e69de29b 100644 --- a/quaddtype/README.md +++ b/quaddtype/README.md @@ -1,16 +0,0 @@ -# quaddtype - -Quad (128-bit) float dtype for numpy - -## Installation - -To install, make sure you have `numpy` nightly installed. Then build without -isolation so that the `quaddtype` can link against the experimental dtype API -headers, which aren't in the latest releases of `numpy`: - -```bash -pip install -i https://pypi.anaconda.org/scientific-python-nightly-wheels/simple numpy -pip install . --no-build-isolation -``` - -Developed with Python 3.11, but 3.9 and 3.10 will probably also work. diff --git a/quaddtype/meson.build b/quaddtype/meson.build index 3c1b0712..f1f9eab1 100644 --- a/quaddtype/meson.build +++ b/quaddtype/meson.build @@ -1,47 +1,56 @@ -project( - 'quaddtype', - 'c', -) +project('quaddtype', 'c', 'cpp', default_options : ['cpp_std=c++17', 'b_pie=true']) py_mod = import('python') py = py_mod.find_installation() +c = meson.get_compiler('c') + +sleef_dep = c.find_library('sleef') +sleefquad_dep = c.find_library('sleefquad') + incdir_numpy = run_command(py, [ '-c', - 'import numpy; print(numpy.get_include())' + 'import numpy; import os; print(os.path.relpath(numpy.get_include()))' ], check: true ).stdout().strip() includes = include_directories( - [ - incdir_numpy, - 'quaddtype/src' - ] + [ + incdir_numpy, + 'quaddtype/src', + ] ) srcs = [ - 'quaddtype/src/umath.c', - 'quaddtype/src/casts.c', - 'quaddtype/src/dtype.c', - 'quaddtype/src/quaddtype_main.c', + 'quaddtype/src/casts.h', + 'quaddtype/src/casts.cpp', + 'quaddtype/src/scalar.h', + 'quaddtype/src/scalar.c', + 'quaddtype/src/dtype.h', + 'quaddtype/src/dtype.c', + 'quaddtype/src/quaddtype_main.c', + 'quaddtype/src/scalar_ops.h', + 'quaddtype/src/scalar_ops.cpp', + 'quaddtype/src/ops.hpp', + 'quaddtype/src/umath.h', + 'quaddtype/src/umath.cpp' ] py.install_sources( - [ - 'quaddtype/__init__.py', - 'quaddtype/quadscalar.py' - ], - subdir: 'quaddtype', - pure: false + [ + 'quaddtype/__init__.py', + ], + subdir: 'quaddtype', + pure: false ) -py.extension_module( - '_quaddtype_main', - srcs, - c_args: ['-g', '-O0'], - install: true, - subdir: 'quaddtype', - include_directories: includes -) +py.extension_module('_quaddtype_main', +srcs, +c_args: ['-g', '-O0', '-lsleef', '-lsleefquad'], +dependencies: [sleef_dep, sleefquad_dep], +install: true, +subdir: 'quaddtype', +include_directories: includes +) \ No newline at end of file diff --git a/quaddtype/pyproject.toml b/quaddtype/pyproject.toml index a2a6fca1..68b61e7f 100644 --- a/quaddtype/pyproject.toml +++ b/quaddtype/pyproject.toml @@ -1,6 +1,6 @@ [build-system] requires = [ - "meson>=0.63.0", + "meson>=1.3.2", "meson-python", "patchelf", "wheel", @@ -13,7 +13,7 @@ name = "quaddtype" description = "Quad (128-bit) float dtype for numpy" version = "0.0.1" readme = 'README.md' -author = "Peyton Murray" +author = "Swayam Singh" requires-python = ">=3.9.0" dependencies = [ "numpy" diff --git a/quaddtype/quaddtype/__init__.py b/quaddtype/quaddtype/__init__.py index c6fe3db1..9f246f73 100644 --- a/quaddtype/quaddtype/__init__.py +++ b/quaddtype/quaddtype/__init__.py @@ -1,8 +1 @@ -# Scalar quantity must be defined _before_ the dtype, so don't isort it. -# During initialization of _quaddtype_main, QuadScalar is imported from this -# (partially initialized) -# module, and therefore has to be defined first. -from .quadscalar import QuadScalar # isort: skip -from ._quaddtype_main import QuadDType - -__all__ = ["QuadScalar", "QuadDType"] +from ._quaddtype_main import QuadPrecDType, QuadPrecision \ No newline at end of file diff --git a/quaddtype/quaddtype/quadscalar.py b/quaddtype/quaddtype/quadscalar.py deleted file mode 100644 index 117c2b34..00000000 --- a/quaddtype/quaddtype/quadscalar.py +++ /dev/null @@ -1,11 +0,0 @@ -"""Quad scalar floating point type for numpy.""" - - -class QuadScalar: - """Quad scalar floating point type.""" - - def __init__(self, value): - self.value = value - - def __repr__(self): - return f"{self.value}" diff --git a/quaddtype/quaddtype/src/casts.c b/quaddtype/quaddtype/src/casts.c deleted file mode 100644 index 999344f5..00000000 --- a/quaddtype/quaddtype/src/casts.c +++ /dev/null @@ -1,157 +0,0 @@ -#include - -#define PY_ARRAY_UNIQUE_SYMBOL quaddtype_ARRAY_API -#define PY_UFUNC_UNIQUE_SYMBOL quaddtype_UFUNC_API -#define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION -#define NPY_TARGET_VERSION NPY_2_0_API_VERSION -#define NO_IMPORT_ARRAY -#define NO_IMPORT_UFUNC -#include "numpy/ndarraytypes.h" -#include "numpy/arrayobject.h" -#include "numpy/dtype_api.h" - -#include "dtype.h" -#include "casts.h" - -// And now the actual cast code! Starting with the "resolver" which tells -// us about cast safety. -// Note also the `view_offset`! It is a way for you to tell NumPy, that this -// cast does not require anything at all, but the cast can simply be done as -// a view. -// For `arr.astype()` it might mean returning a view (eventually, not yet). -// For ufuncs, it already means that they don't have to do a cast at all! -// -// From https://numpy.org/neps/nep-0043-extensible-ufuncs.html#arraymethod: -// resolve_descriptors returns the safety of the operation (casting safety) -static NPY_CASTING -quad_to_quad_resolve_descriptors(PyObject *NPY_UNUSED(self), - PyArray_DTypeMeta *NPY_UNUSED(dtypes[2]), - PyArray_Descr *given_descrs[2], PyArray_Descr *loop_descrs[2], - npy_intp *view_offset) -{ - Py_INCREF(given_descrs[0]); - loop_descrs[0] = given_descrs[0]; - - if (given_descrs[1] == NULL) { - Py_INCREF(given_descrs[0]); - loop_descrs[1] = given_descrs[0]; - } - else { - Py_INCREF(given_descrs[1]); - loop_descrs[1] = given_descrs[1]; - } - - return NPY_SAME_KIND_CASTING; -} - -// Each element is a __float128 element; no casting needed -static int -quad_to_quad_contiguous(PyArrayMethod_Context *NPY_UNUSED(context), char *const data[], - npy_intp const dimensions[], npy_intp const strides[], - NpyAuxData *NPY_UNUSED(auxdata)) -{ - npy_intp N = dimensions[0]; - __float128 *in = (__float128 *)data[0]; - __float128 *out = (__float128 *)data[1]; - - while (N--) { - *out = *in; - out++; - in++; - } - - return 0; -} - -// Elements are strided, e.g. -// -// x = np.linspace(40) -// x[::3] -// -// Therefore the stride needs to be used to increment the pointers inside the loop. -static int -quad_to_quad_strided(PyArrayMethod_Context *NPY_UNUSED(context), char *const data[], - npy_intp const dimensions[], npy_intp const strides[], - NpyAuxData *NPY_UNUSED(auxdata)) -{ - npy_intp N = dimensions[0]; - char *in = data[0]; - char *out = data[1]; - npy_intp in_stride = strides[0]; - npy_intp out_stride = strides[1]; - - while (N--) { - *(__float128 *)out = *(__float128 *)in; - in += in_stride; - out += out_stride; - } - - return 0; -} - -// Arrays are unaligned. -static int -quad_to_quad_unaligned(PyArrayMethod_Context *NPY_UNUSED(context), char *const data[], - npy_intp const dimensions[], npy_intp const strides[], - NpyAuxData *NPY_UNUSED(auxdata)) -{ - npy_intp N = dimensions[0]; - char *in = data[0]; - char *out = data[1]; - npy_intp in_stride = strides[0]; - npy_intp out_stride = strides[1]; - - while (N--) { - memcpy(out, in, sizeof(__float128)); // NOLINT - in += in_stride; - out += out_stride; - } - - return 0; -} - -// Returns the low-level C (strided inner-loop) function which -// performs the actual operation. This method may initially be private, users will be -// able to provide a set of optimized inner-loop functions instead: -// * `strided_inner_loop` -// * `contiguous_inner_loop` -// * `unaligned_strided_loop` -// * ... -static int -quad_to_quad_get_loop(PyArrayMethod_Context *context, int aligned, int NPY_UNUSED(move_references), - const npy_intp *strides, PyArrayMethod_StridedLoop **out_loop, - NpyAuxData **out_transferdata, NPY_ARRAYMETHOD_FLAGS *flags) -{ - int contig = (strides[0] == sizeof(__float128) && strides[1] == sizeof(__float128)); - - if (aligned && contig) - *out_loop = (PyArrayMethod_StridedLoop *)&quad_to_quad_contiguous; - else if (aligned) - *out_loop = (PyArrayMethod_StridedLoop *)&quad_to_quad_strided; - else - *out_loop = (PyArrayMethod_StridedLoop *)&quad_to_quad_unaligned; - - *flags = 0; - return 0; -} - -/* - * NumPy currently allows NULL for the own DType/"cls". For other DTypes - * we would have to fill it in here: - */ -static PyArray_DTypeMeta *QuadToQuadDtypes[2] = {NULL, NULL}; - -static PyType_Slot QuadToQuadSlots[] = { - {NPY_METH_resolve_descriptors, &quad_to_quad_resolve_descriptors}, - {NPY_METH_get_loop, &quad_to_quad_get_loop}, - {0, NULL}}; - -PyArrayMethod_Spec QuadToQuadCastSpec = { - .name = "cast_QuadDType_to_QuadDType", - .nin = 1, - .nout = 1, - .flags = NPY_METH_SUPPORTS_UNALIGNED, - .casting = NPY_SAME_KIND_CASTING, - .dtypes = QuadToQuadDtypes, - .slots = QuadToQuadSlots, -}; diff --git a/quaddtype/quaddtype/src/casts.cpp b/quaddtype/quaddtype/src/casts.cpp new file mode 100644 index 00000000..c096fe5c --- /dev/null +++ b/quaddtype/quaddtype/src/casts.cpp @@ -0,0 +1,468 @@ +#define PY_ARRAY_UNIQUE_SYMBOL QuadPrecType_ARRAY_API +#define PY_UFUNC_UNIQUE_SYMBOL QuadPrecType_UFUNC_API +#define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION +#define NPY_TARGET_VERSION NPY_2_0_API_VERSION +#define NO_IMPORT_ARRAY +#define NO_IMPORT_UFUNC + +extern "C" { +#include + +#include "numpy/arrayobject.h" +#include "numpy/ndarraytypes.h" +#include "numpy/dtype_api.h" +} +#include "sleef.h" +#include "sleefquad.h" + +#include "scalar.h" +#include "casts.h" +#include "dtype.h" + +#define NUM_CASTS 29 // 14 to_casts + 14 from_casts + 1 quad_to_quad + +static NPY_CASTING +quad_to_quad_resolve_descriptors(PyObject *NPY_UNUSED(self), + PyArray_DTypeMeta *NPY_UNUSED(dtypes[2]), + QuadPrecDTypeObject *given_descrs[2], + QuadPrecDTypeObject *loop_descrs[2], npy_intp *view_offset) +{ + Py_INCREF(given_descrs[0]); + loop_descrs[0] = given_descrs[0]; + + if (given_descrs[1] == NULL) { + Py_INCREF(given_descrs[0]); + loop_descrs[1] = given_descrs[0]; + } + else { + Py_INCREF(given_descrs[1]); + loop_descrs[1] = given_descrs[1]; + } + + *view_offset = 0; + return NPY_NO_CASTING; +} + +static int +quad_to_quad_strided_loop(PyArrayMethod_Context *context, char *const data[], + npy_intp const dimensions[], npy_intp const strides[], + void *NPY_UNUSED(auxdata)) +{ + npy_intp N = dimensions[0]; + char *in_ptr = data[0]; + char *out_ptr = data[1]; + + npy_intp in_stride = strides[0]; + npy_intp out_stride = strides[1]; + + while (N--) { + memcpy(out_ptr, in_ptr, sizeof(Sleef_quad)); + in_ptr += in_stride; + out_ptr += out_stride; + } + return 0; +} + +// Casting from other types to QuadDType + +template +static inline Sleef_quad +to_quad(T x); + +template <> +inline Sleef_quad +to_quad(npy_bool x) +{ + return x ? Sleef_cast_from_doubleq1(1.0) : Sleef_cast_from_doubleq1(0.0); +} +template <> +inline Sleef_quad +to_quad(npy_byte x) +{ + return Sleef_cast_from_int64q1(x); +} +// template <> +// inline Sleef_quad +// to_quad(npy_ubyte x) +// { +// return Sleef_cast_from_uint64q1(x); +// } +template <> +inline Sleef_quad +to_quad(npy_short x) +{ + return Sleef_cast_from_int64q1(x); +} +template <> +inline Sleef_quad +to_quad(npy_ushort x) +{ + return Sleef_cast_from_uint64q1(x); +} +template <> +inline Sleef_quad +to_quad(npy_int x) +{ + return Sleef_cast_from_int64q1(x); +} +template <> +inline Sleef_quad +to_quad(npy_uint x) +{ + return Sleef_cast_from_uint64q1(x); +} +template <> +inline Sleef_quad +to_quad(npy_long x) +{ + return Sleef_cast_from_int64q1(x); +} +template <> +inline Sleef_quad +to_quad(npy_ulong x) +{ + return Sleef_cast_from_uint64q1(x); +} +template <> +inline Sleef_quad +to_quad(npy_longlong x) +{ + return Sleef_cast_from_int64q1(x); +} +template <> +inline Sleef_quad +to_quad(npy_ulonglong x) +{ + return Sleef_cast_from_uint64q1(x); +} +template <> +inline Sleef_quad +to_quad(float x) +{ + return Sleef_cast_from_doubleq1(x); +} +template <> +inline Sleef_quad +to_quad(double x) +{ + return Sleef_cast_from_doubleq1(x); +} +template <> +inline Sleef_quad +to_quad(long double x) +{ + return Sleef_cast_from_doubleq1(x); +} + +template +static NPY_CASTING +numpy_to_quad_resolve_descriptors(PyObject *NPY_UNUSED(self), PyArray_DTypeMeta *dtypes[2], + PyArray_Descr *given_descrs[2], PyArray_Descr *loop_descrs[2], + npy_intp *view_offset) +{ + if (given_descrs[1] == NULL) { + loop_descrs[1] = (PyArray_Descr *)new_quaddtype_instance(); + if (loop_descrs[1] == nullptr) { + return (NPY_CASTING)-1; + } + } + else { + Py_INCREF(given_descrs[1]); + loop_descrs[1] = given_descrs[1]; + } + + loop_descrs[0] = PyArray_GetDefaultDescr(dtypes[0]); + // *view_offset = 0; + return NPY_SAFE_CASTING; +} + +template +static int +numpy_to_quad_strided_loop(PyArrayMethod_Context *context, char *const data[], + npy_intp const dimensions[], npy_intp const strides[], + void *NPY_UNUSED(auxdata)) +{ + npy_intp N = dimensions[0]; + char *in_ptr = data[0]; + char *out_ptr = data[1]; + + while (N--) { + T in_val; + Sleef_quad out_val; + + memcpy(&in_val, in_ptr, sizeof(T)); + out_val = to_quad(in_val); + memcpy(out_ptr, &out_val, sizeof(Sleef_quad)); + + in_ptr += strides[0]; + out_ptr += strides[1]; + } + return 0; +} + +// Casting from QuadDType to other types + +template +static inline T +from_quad(Sleef_quad x); + +template <> +inline npy_bool +from_quad(Sleef_quad x) +{ + return Sleef_cast_to_int64q1(x) != 0; +} +template <> +inline npy_byte +from_quad(Sleef_quad x) +{ + return (npy_byte)Sleef_cast_to_int64q1(x); +} +// template <> +// inline npy_ubyte +// from_quad(Sleef_quad x) +// { +// return (npy_ubyte)Sleef_cast_to_uint64q1(x); +// } +template <> +inline npy_short +from_quad(Sleef_quad x) +{ + return (npy_short)Sleef_cast_to_int64q1(x); +} +template <> +inline npy_ushort +from_quad(Sleef_quad x) +{ + return (npy_ushort)Sleef_cast_to_uint64q1(x); +} +template <> +inline npy_int +from_quad(Sleef_quad x) +{ + return (npy_int)Sleef_cast_to_int64q1(x); +} +template <> +inline npy_uint +from_quad(Sleef_quad x) +{ + return (npy_uint)Sleef_cast_to_uint64q1(x); +} +template <> +inline npy_long +from_quad(Sleef_quad x) +{ + return (npy_long)Sleef_cast_to_int64q1(x); +} +template <> +inline npy_ulong +from_quad(Sleef_quad x) +{ + return (npy_ulong)Sleef_cast_to_uint64q1(x); +} +template <> +inline npy_longlong +from_quad(Sleef_quad x) +{ + return Sleef_cast_to_int64q1(x); +} +template <> +inline npy_ulonglong +from_quad(Sleef_quad x) +{ + return Sleef_cast_to_uint64q1(x); +} +template <> +inline float +from_quad(Sleef_quad x) +{ + return Sleef_cast_to_doubleq1(x); +} +template <> +inline double +from_quad(Sleef_quad x) +{ + return Sleef_cast_to_doubleq1(x); +} +template <> +inline long double +from_quad(Sleef_quad x) +{ + return Sleef_cast_to_doubleq1(x); +} + +template +static NPY_CASTING +quad_to_numpy_resolve_descriptors(PyObject *NPY_UNUSED(self), PyArray_DTypeMeta *dtypes[2], + PyArray_Descr *given_descrs[2], PyArray_Descr *loop_descrs[2], + npy_intp *view_offset) +{ + Py_INCREF(given_descrs[0]); + loop_descrs[0] = given_descrs[0]; + + loop_descrs[1] = PyArray_GetDefaultDescr(dtypes[1]); + // *view_offset = 0; + return NPY_UNSAFE_CASTING; +} + +template +static int +quad_to_numpy_strided_loop(PyArrayMethod_Context *context, char *const data[], + npy_intp const dimensions[], npy_intp const strides[], + void *NPY_UNUSED(auxdata)) +{ + npy_intp N = dimensions[0]; + char *in_ptr = data[0]; + char *out_ptr = data[1]; + + while (N--) { + Sleef_quad in_val = *(Sleef_quad *)in_ptr; + T *out_val = (T *)out_ptr; + *out_val = from_quad(in_val); + + in_ptr += strides[0]; + out_ptr += strides[1]; + } + return 0; +} + +static PyArrayMethod_Spec *specs[NUM_CASTS + 1]; // +1 for NULL terminator +static size_t spec_count = 0; + +void +add_spec(PyArrayMethod_Spec *spec) +{ + if (spec_count < NUM_CASTS) { + specs[spec_count++] = spec; + } +} + +// functions to add casts +template +void +add_cast_from(PyArray_DTypeMeta *to) +{ + PyArray_DTypeMeta **dtypes = new PyArray_DTypeMeta *[2]{nullptr, to}; + + PyType_Slot *slots = new PyType_Slot[3]{ + {NPY_METH_resolve_descriptors, (void *)&quad_to_numpy_resolve_descriptors}, + {NPY_METH_strided_loop, (void *)&quad_to_numpy_strided_loop}, + {0, nullptr}}; + + PyArrayMethod_Spec *spec = new PyArrayMethod_Spec{ + .name = "cast_QuadPrec_to_NumPy", + .nin = 1, + .nout = 1, + .casting = NPY_UNSAFE_CASTING, + .flags = (NPY_ARRAYMETHOD_FLAGS)0, + .dtypes = dtypes, + .slots = slots, + }; + add_spec(spec); +} + +template +void +add_cast_to(PyArray_DTypeMeta *from) +{ + PyArray_DTypeMeta **dtypes = new PyArray_DTypeMeta *[2]{from, nullptr}; + + PyType_Slot *slots = new PyType_Slot[3]{ + {NPY_METH_resolve_descriptors, (void *)&numpy_to_quad_resolve_descriptors}, + {NPY_METH_strided_loop, (void *)&numpy_to_quad_strided_loop}, + {0, nullptr}}; + + PyArrayMethod_Spec *spec = new PyArrayMethod_Spec{ + .name = "cast_NumPy_to_QuadPrec", + .nin = 1, + .nout = 1, + .casting = NPY_SAFE_CASTING, + .flags = (NPY_ARRAYMETHOD_FLAGS)0, + .dtypes = dtypes, + .slots = slots, + }; + + add_spec(spec); +} + +PyArrayMethod_Spec ** +init_casts_internal(void) +{ + PyArray_DTypeMeta **quad2quad_dtypes = new PyArray_DTypeMeta *[2]{nullptr, nullptr}; + PyType_Slot *quad2quad_slots = new PyType_Slot[4]{ + {NPY_METH_resolve_descriptors, (void *)&quad_to_quad_resolve_descriptors}, + {NPY_METH_strided_loop, (void *)&quad_to_quad_strided_loop}, + {NPY_METH_unaligned_strided_loop, (void *)&quad_to_quad_strided_loop}, + {0, nullptr}}; + + PyArrayMethod_Spec *quad2quad_spec = new PyArrayMethod_Spec{ + .name = "cast_QuadPrec_to_QuadPrec", + .nin = 1, + .nout = 1, + .casting = NPY_NO_CASTING, + .flags = NPY_METH_SUPPORTS_UNALIGNED, + .dtypes = quad2quad_dtypes, + .slots = quad2quad_slots, + }; + + add_spec(quad2quad_spec); + + add_cast_to(&PyArray_BoolDType); + add_cast_to(&PyArray_ByteDType); + add_cast_to(&PyArray_UByteDType); + add_cast_to(&PyArray_ShortDType); + add_cast_to(&PyArray_UShortDType); + add_cast_to(&PyArray_IntDType); + add_cast_to(&PyArray_UIntDType); + add_cast_to(&PyArray_LongDType); + add_cast_to(&PyArray_ULongDType); + add_cast_to(&PyArray_LongLongDType); + add_cast_to(&PyArray_ULongLongDType); + add_cast_to(&PyArray_FloatDType); + add_cast_to(&PyArray_DoubleDType); + add_cast_to(&PyArray_LongDoubleDType); + + add_cast_from(&PyArray_BoolDType); + add_cast_from(&PyArray_ByteDType); + add_cast_from(&PyArray_UByteDType); + add_cast_from(&PyArray_ShortDType); + add_cast_from(&PyArray_UShortDType); + add_cast_from(&PyArray_IntDType); + add_cast_from(&PyArray_UIntDType); + add_cast_from(&PyArray_LongDType); + add_cast_from(&PyArray_ULongDType); + add_cast_from(&PyArray_LongLongDType); + add_cast_from(&PyArray_ULongLongDType); + add_cast_from(&PyArray_FloatDType); + add_cast_from(&PyArray_DoubleDType); + add_cast_from(&PyArray_LongDoubleDType); + + specs[spec_count] = nullptr; + return specs; +} + +PyArrayMethod_Spec ** +init_casts(void) +{ + try { + return init_casts_internal(); + } + catch (int e) { + PyErr_NoMemory(); + return nullptr; + } +} + +void +free_casts(void) +{ + for (auto cast : specs) { + if (cast == nullptr) { + continue; + } + delete cast->dtypes; + delete cast->slots; + delete cast; + } + spec_count = 0; +} diff --git a/quaddtype/quaddtype/src/casts.h b/quaddtype/quaddtype/src/casts.h index de775617..7a7cd9f5 100644 --- a/quaddtype/quaddtype/src/casts.h +++ b/quaddtype/quaddtype/src/casts.h @@ -1,7 +1,21 @@ -#ifndef _NPY_CASTS_H -#define _NPY_CASTS_H +#ifndef _QUADDTYPE_CASTS_H +#define _QUADDTYPE_CASTS_H -extern PyArrayMethod_Spec QuadToQuadCastSpec; -extern PyArrayMethod_Spec QuadToFloat128CastSpec; +#include +#include "numpy/dtype_api.h" -#endif /* _NPY_CASTS_H */ +#ifdef __cplusplus +extern "C" { +#endif + +PyArrayMethod_Spec ** +init_casts(void); + +void +free_casts(void); + +#ifdef __cplusplus +} +#endif + +#endif \ No newline at end of file diff --git a/quaddtype/quaddtype/src/dtype.c b/quaddtype/quaddtype/src/dtype.c index 63a6840e..a6f485bf 100644 --- a/quaddtype/quaddtype/src/dtype.c +++ b/quaddtype/quaddtype/src/dtype.c @@ -1,197 +1,216 @@ #include +#include +#include -#define PY_ARRAY_UNIQUE_SYMBOL quaddtype_ARRAY_API -#define PY_UFUNC_UNIQUE_SYMBOL quaddtype_UFUNC_API +#define PY_ARRAY_UNIQUE_SYMBOL QuadPrecType_ARRAY_API +#define PY_UFUNC_UNIQUE_SYMBOL QuadPrecType_UFUNC_API #define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION #define NPY_TARGET_VERSION NPY_2_0_API_VERSION #define NO_IMPORT_ARRAY #define NO_IMPORT_UFUNC -#include "numpy/ndarraytypes.h" #include "numpy/arrayobject.h" +#include "numpy/ndarraytypes.h" #include "numpy/dtype_api.h" -#include "dtype.h" -#include "abstract.h" +#include "scalar.h" #include "casts.h" +#include "dtype.h" -PyTypeObject *QuadScalar_Type = NULL; - -QuadDTypeObject * -new_quaddtype_instance(void) +static inline int quad_load(Sleef_quad *x, char *data_ptr) { - QuadDTypeObject *new = - (QuadDTypeObject *)PyArrayDescr_Type.tp_new((PyTypeObject *)&QuadDType, NULL, NULL); - if (new == NULL) { - return NULL; + if (data_ptr == NULL || x == NULL) + { + return -1; } - - new->base.elsize = sizeof(__float128); - new->base.alignment = _Alignof(__float128); - return new; + *x = *(Sleef_quad *)data_ptr; + return 0; } -// Take an python double and put a copy into the array -static int -quad_setitem(QuadDTypeObject *descr, PyObject *obj, char *dataptr) +static inline int quad_store(char *data_ptr, Sleef_quad x) { - __float128 val = (__float128)PyFloat_AsDouble(obj); - memcpy(dataptr, &val, sizeof(__float128)); + if (data_ptr == NULL) + { + return -1; + } + *(Sleef_quad *)data_ptr = x; return 0; } -static PyObject * -quad_getitem(QuadDTypeObject *descr, char *dataptr) +QuadPrecDTypeObject * new_quaddtype_instance(void) { - __float128 val; - memcpy(&val, dataptr, sizeof(__float128)); - - PyObject *val_obj = PyFloat_FromDouble((double)val); - if (val_obj == NULL) { + QuadPrecDTypeObject *new = (QuadPrecDTypeObject *)PyArrayDescr_Type.tp_new((PyTypeObject *)&QuadPrecDType, NULL, NULL); + if (new == NULL) { return NULL; } + new->base.elsize = sizeof(Sleef_quad); + new->base.alignment = _Alignof(Sleef_quad); - // Need to create a new QuadScalar instance here and return that... - PyObject *res = PyObject_CallFunctionObjArgs((PyObject *)QuadScalar_Type, val_obj, NULL); - if (res == NULL) { - return NULL; - } - Py_DECREF(val_obj); - return res; + return new; } -// For two instances of the same dtype, both have the same precision. Return -// self. -static QuadDTypeObject * -common_instance(QuadDTypeObject *self, QuadDTypeObject *other) +static QuadPrecDTypeObject * ensure_canonical(QuadPrecDTypeObject *self) { + Py_INCREF(self); return self; } -// When dtypes are mixed, find a "common" dtype for the two which can hold -// content of both without loss of information. I guess this should return a -// 256-bit float dtype? Since this isn't natively supported by any platform, -// just return another 128-bit float dtype. -static PyArray_DTypeMeta * -common_dtype(PyArray_DTypeMeta *self, PyArray_DTypeMeta *other) +static QuadPrecDTypeObject * common_instance(QuadPrecDTypeObject *dtype1, QuadPrecDTypeObject *dtype2) +{ + Py_INCREF(dtype1); + return dtype1; +} + + +static PyArray_DTypeMeta * common_dtype(PyArray_DTypeMeta *cls, PyArray_DTypeMeta *other) { - /* - * Typenum is useful for NumPy, but there it can still be convenient. - * (New-style user dtypes will probably get -1 as type number...) - */ - if (other->type_num >= 0 && PyTypeNum_ISNUMBER(other->type_num) && - !PyTypeNum_ISCOMPLEX(other->type_num)) { - // float128 is the biggest natively supported float. Return it in all - // cases where other is a number (and not complex). - Py_INCREF(self); - return self; + // Promote integer and floating-point types to QuadPrecDType + if (other->type_num >= 0 && + (PyTypeNum_ISINTEGER(other->type_num) || + PyTypeNum_ISFLOAT(other->type_num))) { + Py_INCREF(cls); + return cls; + } + // Don't promote complex types + if (PyTypeNum_ISCOMPLEX(other->type_num)) { + Py_INCREF(Py_NotImplemented); + return (PyArray_DTypeMeta *)Py_NotImplemented; } - // Revert to object dtype in all other cases. Py_INCREF(Py_NotImplemented); return (PyArray_DTypeMeta *)Py_NotImplemented; } -// Expected to have this, and that it does an incref; see NEP42 -// Without this you'll get weird memory corruption bugs in the casting code -static QuadDTypeObject * -quaddtype_ensure_canonical(QuadDTypeObject *self) +static PyArray_Descr * +quadprec_discover_descriptor_from_pyobject(PyArray_DTypeMeta *NPY_UNUSED(cls), PyObject *obj) { - Py_INCREF(self); - return self; + if (Py_TYPE(obj) != &QuadPrecision_Type) + { + PyErr_SetString(PyExc_TypeError, "Can only store QuadPrecision in a QuadPrecDType array."); + return NULL; + } + return (PyArray_Descr *)new_quaddtype_instance(); } -static PyArray_Descr * -quad_discover_descriptor_from_pyobject(PyArray_DTypeMeta *NPY_UNUSED(cls), PyObject *obj) +static int quadprec_setitem(QuadPrecDTypeObject *descr, PyObject *obj, char *dataptr) { - if (Py_TYPE(obj) != QuadScalar_Type) { - PyErr_SetString(PyExc_TypeError, "Can only store QuadScalars in a QuadDType array."); - return NULL; + QuadPrecisionObject *value; + if (PyObject_TypeCheck(obj, &QuadPrecision_Type)) + { + Py_INCREF(obj); + value = (QuadPrecisionObject *)obj; + } + else + { + value = QuadPrecision_from_object(obj); + if (value == NULL) { + return -1; + } } - // Get the dtype attribute from the object. - return (PyArray_Descr *)PyObject_GetAttrString(obj, "dtype"); + if (quad_store(dataptr, value->quad.value) < 0) + { + Py_DECREF(value); + char error_msg[100]; + snprintf(error_msg, sizeof(error_msg), "Invalid memory location %p", (void*)dataptr); + PyErr_SetString(PyExc_ValueError, error_msg); + return -1; + } + + Py_DECREF(value); + return 0; } -static PyType_Slot QuadDType_Slots[] = { - {NPY_DT_common_instance, &common_instance}, - {NPY_DT_common_dtype, &common_dtype}, - {NPY_DT_discover_descr_from_pyobject, &quad_discover_descriptor_from_pyobject}, - /* The header is wrong on main :(, so we add 1 */ - {NPY_DT_setitem, &quad_setitem}, - {NPY_DT_getitem, &quad_getitem}, - {NPY_DT_ensure_canonical, &quaddtype_ensure_canonical}, - {0, NULL}}; - -/* - * The following defines everything type object related (i.e. not NumPy - * specific). - * - * Note that this function is by default called without any arguments to fetch - * a default version of the descriptor (in principle at least). During init - * we fill in `cls->singleton` though for the dimensionless unit. - */ -static PyObject * -quaddtype_new(PyTypeObject *NPY_UNUSED(cls), PyObject *args, PyObject *kwargs) +static PyObject * quadprec_getitem(QuadPrecDTypeObject *descr, char *dataptr) { - return (PyObject *)new_quaddtype_instance(); + QuadPrecisionObject *new = QuadPrecision_raw_new(); + if (!new) + { + return NULL; + } + if (quad_load(&new->quad.value, dataptr) < 0) + { + Py_DECREF(new); + char error_msg[100]; + snprintf(error_msg, sizeof(error_msg), "Invalid memory location %p", (void*)dataptr); + PyErr_SetString(PyExc_ValueError, error_msg); + return NULL; + } + return (PyObject *)new; +} + +static PyArray_Descr *quadprec_default_descr(PyArray_DTypeMeta *NPY_UNUSED(cls)) +{ + return (PyArray_Descr *)new_quaddtype_instance(); } -static void -quaddtype_dealloc(QuadDTypeObject *self) +static PyType_Slot QuadPrecDType_Slots[] = { - PyArrayDescr_Type.tp_dealloc((PyObject *)self); + {NPY_DT_ensure_canonical, &ensure_canonical}, + {NPY_DT_common_instance, &common_instance}, + {NPY_DT_common_dtype, &common_dtype}, + {NPY_DT_discover_descr_from_pyobject, &quadprec_discover_descriptor_from_pyobject}, + {NPY_DT_setitem, &quadprec_setitem}, + {NPY_DT_getitem, &quadprec_getitem}, + {NPY_DT_default_descr, &quadprec_default_descr}, + {0, NULL} +}; + + +static PyObject * QuadPrecDType_new(PyTypeObject *NPY_UNUSED(cls), PyObject *args, PyObject *kwds) +{ + if (PyTuple_GET_SIZE(args) != 0 || (kwds != NULL && PyDict_Size(kwds) != 0)) { + PyErr_SetString(PyExc_TypeError, + "QuadPrecDType takes no arguments"); + return NULL; + } + + return (PyObject *)new_quaddtype_instance(); } -static PyObject * -quaddtype_repr(QuadDTypeObject *self) +static PyObject * QuadPrecDType_repr(QuadPrecDTypeObject *self) { - PyObject *res = PyUnicode_FromString("This is a quad (128-bit float) dtype."); - return res; + return PyUnicode_FromString("QuadPrecDType()"); } -// These are the basic things that you need to create a Python Type/Class in C. -// However, there is a slight difference here because we create a -// PyArray_DTypeMeta, which is a larger struct than a typical type. -// (This should get a bit nicer eventually with Python >3.11.) -PyArray_DTypeMeta QuadDType = { - {{ - PyVarObject_HEAD_INIT(NULL, 0).tp_name = "quaddtype.QuadDType", - .tp_basicsize = sizeof(QuadDTypeObject), - .tp_new = quaddtype_new, - .tp_dealloc = (destructor)quaddtype_dealloc, - .tp_repr = (reprfunc)quaddtype_repr, - .tp_str = (reprfunc)quaddtype_repr, - }}, - /* rest, filled in during DTypeMeta initialization */ +PyArray_DTypeMeta QuadPrecDType = { + {{ + PyVarObject_HEAD_INIT(NULL, 0) + .tp_name = "QuadPrecDType.QuadPrecDType", + .tp_basicsize = sizeof(QuadPrecDTypeObject), + .tp_new = QuadPrecDType_new, + .tp_repr = (reprfunc)QuadPrecDType_repr, + .tp_str = (reprfunc)QuadPrecDType_repr, + }}, }; -int -init_quad_dtype(void) +int init_quadprec_dtype(void) { - // To create our DType, we have to use a "Spec" that tells NumPy how to - // do it. You first have to create a static type, but see the note there! - PyArrayMethod_Spec *casts[] = { - &QuadToQuadCastSpec, - NULL, - }; + PyArrayMethod_Spec **casts = init_casts(); + if (!casts) + return -1; - PyArrayDTypeMeta_Spec QuadDType_DTypeSpec = { - .flags = NPY_DT_NUMERIC, - .casts = casts, - .typeobj = QuadScalar_Type, - .slots = QuadDType_Slots, + PyArrayDTypeMeta_Spec QuadPrecDType_DTypeSpec = { + .flags = NPY_DT_NUMERIC, + .casts = casts, + .typeobj = &QuadPrecision_Type, + .slots = QuadPrecDType_Slots, }; - ((PyObject *)&QuadDType)->ob_type = &PyArrayDTypeMeta_Type; - ((PyTypeObject *)&QuadDType)->tp_base = &PyArrayDescr_Type; - if (PyType_Ready((PyTypeObject *)&QuadDType) < 0) { + ((PyObject *)&QuadPrecDType)->ob_type = &PyArrayDTypeMeta_Type; + + ((PyTypeObject *)&QuadPrecDType)->tp_base = &PyArrayDescr_Type; + + if (PyType_Ready((PyTypeObject *)&QuadPrecDType) < 0) + { return -1; } - if (PyArrayInitDTypeMeta_FromSpec(&QuadDType, &QuadDType_DTypeSpec) < 0) { + if (PyArrayInitDTypeMeta_FromSpec(&QuadPrecDType, &QuadPrecDType_DTypeSpec) < 0) + { return -1; } - QuadDType.singleton = PyArray_GetDefaultDescr(&QuadDType); + free_casts(); + return 0; -} +} \ No newline at end of file diff --git a/quaddtype/quaddtype/src/dtype.h b/quaddtype/quaddtype/src/dtype.h index 104ed3e3..162d5714 100644 --- a/quaddtype/quaddtype/src/dtype.h +++ b/quaddtype/quaddtype/src/dtype.h @@ -1,17 +1,30 @@ -#ifndef _NPY_DTYPE_H -#define _NPY_DTYPE_H +#ifndef _QUADDTYPE_DTYPE_H +#define _QUADDTYPE_DTYPE_H +#ifdef __cplusplus +extern "C" { +#endif -typedef struct { +#include +#include +#include +#include + +#include "scalar.h" + +typedef struct +{ PyArray_Descr base; -} QuadDTypeObject; + +} QuadPrecDTypeObject; + +extern PyArray_DTypeMeta QuadPrecDType; -extern PyArray_DTypeMeta QuadDType; -extern PyTypeObject *QuadScalar_Type; +QuadPrecDTypeObject * new_quaddtype_instance(void); -QuadDTypeObject * -new_quaddtype_instance(void); +int init_quadprec_dtype(void); -int -init_quad_dtype(void); +#ifdef __cplusplus +} +#endif -#endif /*_NPY_DTYPE_H*/ +#endif \ No newline at end of file diff --git a/quaddtype/quaddtype/src/ops.hpp b/quaddtype/quaddtype/src/ops.hpp new file mode 100644 index 00000000..6a3511c4 --- /dev/null +++ b/quaddtype/quaddtype/src/ops.hpp @@ -0,0 +1,207 @@ +#include +#include + +typedef int (*unary_op_def)(Sleef_quad *, Sleef_quad *); + +static inline int +quad_negative(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_negq1(*op); + return 0; +} + +static int +quad_positive(Sleef_quad *op, Sleef_quad *out) +{ + *out = *op; + return 0; +} + +static inline int +quad_absolute(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_fabsq1(*op); + return 0; +} + +static inline int +quad_rint(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_rintq1(*op); + return 0; +} + +static inline int +quad_trunc(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_truncq1(*op); + return 0; +} + +static inline int +quad_floor(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_floorq1(*op); + return 0; +} + +static inline int +quad_ceil(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_ceilq1(*op); + return 0; +} + +static inline int +quad_sqrt(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_sqrtq1_u05(*op); + return 0; +} + +static inline int +quad_square(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_mulq1_u05(*op, *op); + return 0; +} + +static inline int +quad_log(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_logq1_u10(*op); + return 0; +} + +static inline int +quad_log2(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_log2q1_u10(*op); + return 0; +} + +static inline int +quad_log10(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_log10q1_u10(*op); + return 0; +} + +static inline int +quad_log1p(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_log1pq1_u10(*op); + return 0; +} + +static inline int +quad_exp(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_expq1_u10(*op); + return 0; +} + +static inline int +quad_exp2(Sleef_quad *op, Sleef_quad *out) +{ + *out = Sleef_exp2q1_u10(*op); + return 0; +} + +// binary ops +typedef int (*binop_def)(Sleef_quad *, Sleef_quad *, Sleef_quad *); + +static inline int +quad_add(Sleef_quad *out, Sleef_quad *in1, Sleef_quad *in2) +{ + *out = Sleef_addq1_u05(*in1, *in2); + return 0; +} + +static inline int +quad_sub(Sleef_quad *out, Sleef_quad *in1, Sleef_quad *in2) +{ + *out = Sleef_subq1_u05(*in1, *in2); + return 0; +} + +static inline int +quad_mul(Sleef_quad *res, Sleef_quad *a, Sleef_quad *b) +{ + *res = Sleef_mulq1_u05(*a, *b); + return 0; +} + +static inline int +quad_div(Sleef_quad *res, Sleef_quad *a, Sleef_quad *b) +{ + *res = Sleef_divq1_u05(*a, *b); + return 0; +} + +static inline int +quad_pow(Sleef_quad *res, Sleef_quad *a, Sleef_quad *b) +{ + *res = Sleef_powq1_u10(*a, *b); + return 0; +} + +static inline int +quad_mod(Sleef_quad *res, Sleef_quad *a, Sleef_quad *b) +{ + *res = Sleef_fmodq1(*a, *b); + return 0; +} + +static inline int +quad_minimum(Sleef_quad *out, Sleef_quad *in1, Sleef_quad *in2) +{ + *out = Sleef_icmpleq1(*in1, *in2) ? *in1 : *in2; + return 0; +} + +static inline int +quad_maximum(Sleef_quad *out, Sleef_quad *in1, Sleef_quad *in2) +{ + *out = Sleef_icmpgeq1(*in1, *in2) ? *in1 : *in2; + return 0; +} + +// comparison functions +typedef npy_bool (*cmp_def)(const Sleef_quad *, const Sleef_quad *); + +static inline npy_bool +quad_equal(const Sleef_quad *a, const Sleef_quad *b) +{ + return Sleef_icmpeqq1(*a, *b); +} + +static inline npy_bool +quad_notequal(const Sleef_quad *a, const Sleef_quad *b) +{ + return Sleef_icmpneq1(*a, *b); +} + +static inline npy_bool +quad_less(const Sleef_quad *a, const Sleef_quad *b) +{ + return Sleef_icmpltq1(*a, *b); +} + +static inline npy_bool +quad_lessequal(const Sleef_quad *a, const Sleef_quad *b) +{ + return Sleef_icmpleq1(*a, *b); +} + +static inline npy_bool +quad_greater(const Sleef_quad *a, const Sleef_quad *b) +{ + return Sleef_icmpgtq1(*a, *b); +} + +static inline npy_bool +quad_greaterequal(const Sleef_quad *a, const Sleef_quad *b) +{ + return Sleef_icmpgeq1(*a, *b); +} \ No newline at end of file diff --git a/quaddtype/quaddtype/src/quaddtype_main.c b/quaddtype/quaddtype/src/quaddtype_main.c index d9adeff5..098c0749 100644 --- a/quaddtype/quaddtype/src/quaddtype_main.c +++ b/quaddtype/quaddtype/src/quaddtype_main.c @@ -1,65 +1,55 @@ #include -#define PY_ARRAY_UNIQUE_SYMBOL quaddtype_ARRAY_API -#define PY_UFUNC_UNIQUE_SYMBOL quaddtype_UFUNC_API +#define PY_ARRAY_UNIQUE_SYMBOL QuadPrecType_ARRAY_API +#define PY_UFUNC_UNIQUE_SYMBOL QuadPrecType_UFUNC_API #define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION +#define NPY_TARGET_VERSION NPY_2_0_API_VERSION + #include "numpy/arrayobject.h" -#include "numpy/ufuncobject.h" #include "numpy/dtype_api.h" +#include "numpy/ufuncobject.h" #include "dtype.h" #include "umath.h" static struct PyModuleDef moduledef = { PyModuleDef_HEAD_INIT, - .m_name = "quaddtype_main", - .m_doc = "Quad (128-bit) floating point experimental numpy dtype", + .m_name = "_quaddtype_main", + .m_doc = "Quad (128-bit) floating point Data Type for Numpy", .m_size = -1, }; -// Initialize the python module PyMODINIT_FUNC PyInit__quaddtype_main(void) { import_array(); import_umath(); - PyObject *m = PyModule_Create(&moduledef); - if (m == NULL) { - PyErr_SetString(PyExc_ImportError, "Unable to create the quaddtype_main module."); + if (!m) + { return NULL; } - PyObject *mod = PyImport_ImportModule("quaddtype"); - if (mod == NULL) { - PyErr_SetString(PyExc_ImportError, "Unable to import the quaddtype module."); + if (init_quadprecision_scalar() < 0) goto error; - } - QuadScalar_Type = (PyTypeObject *)PyObject_GetAttrString(mod, "QuadScalar"); - Py_DECREF(mod); - if (QuadScalar_Type == NULL) { - PyErr_SetString(PyExc_AttributeError, - "Unable to find QuadScalar attribute in the " - "quaddtype_main module."); + + if(PyModule_AddObject(m, "QuadPrecision", (PyObject *)&QuadPrecision_Type) < 0) goto error; - } - if (init_quad_dtype() < 0) { - PyErr_SetString(PyExc_AttributeError, "QuadDType failed to initialize."); + + if(init_quadprec_dtype() < 0) goto error; - } - if (PyModule_AddObject(m, "QuadDType", (PyObject *)&QuadDType) < 0) { - PyErr_SetString(PyExc_TypeError, "Failed to add QuadDType to the quaddtype_main module."); + + if(PyModule_AddObject(m, "QuadPrecDType", (PyObject *)&QuadPrecDType) < 0) goto error; - } - if (init_multiply_ufunc() == -1) { - PyErr_SetString(PyExc_TypeError, "Failed to initialize the quadscalar multiply ufunc."); + if (init_quad_umath() < 0) { goto error; } return m; + error: Py_DECREF(m); return NULL; -} +} \ No newline at end of file diff --git a/quaddtype/quaddtype/src/scalar.c b/quaddtype/quaddtype/src/scalar.c new file mode 100644 index 00000000..15f727d8 --- /dev/null +++ b/quaddtype/quaddtype/src/scalar.c @@ -0,0 +1,126 @@ +#include +#include +#include + +#define PY_ARRAY_UNIQUE_SYMBOL QuadPrecType_ARRAY_API +#define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION +#define NO_IMPORT_ARRAY + +#include "numpy/arrayobject.h" +#include "numpy/ndarraytypes.h" +#include "numpy/dtype_api.h" + +#include "scalar.h" +#include "scalar_ops.h" + + +// static PyTypeObject QuadPrecision_Type; + + +QuadPrecisionObject * QuadPrecision_raw_new(void) +{ + QuadPrecisionObject * new = PyObject_New(QuadPrecisionObject, &QuadPrecision_Type); + if(!new) + return NULL; + new->quad.value = Sleef_cast_from_doubleq1(0.0); // initialize to 0 + return new; +} + +QuadPrecisionObject * QuadPrecision_from_object(PyObject * value) +{ + QuadPrecisionObject *self = QuadPrecision_raw_new(); + if(!self) + return NULL; + + if(PyFloat_Check(value)) + self->quad.value = Sleef_cast_from_doubleq1(PyFloat_AsDouble(value)); + + else if(PyUnicode_CheckExact(value)) + { + const char * s = PyUnicode_AsUTF8(value); + char *endptr = NULL; + self->quad.value = Sleef_strtoq(s, &endptr); + if (*endptr != '\0' || endptr == s) + { + PyErr_SetString(PyExc_ValueError, "Unable to parse string to QuadPrecision"); + Py_DECREF(self); + return NULL; + } + } + else if(PyLong_Check(value)) + { + long int val = PyLong_AsLong(value); + if (val == -1) + { + PyErr_SetString(PyExc_OverflowError, "Overflow Error, value out of range"); + Py_DECREF(self); + return NULL; + } + self->quad.value = Sleef_cast_from_int64q1(val); + } + else + { + PyErr_SetString(PyExc_TypeError, "QuadPrecision value must be a float, int or string"); + Py_DECREF(self); + return NULL; + } + + return self; +} + +static PyObject * +QuadPrecision_new(PyTypeObject *cls, PyObject *args, PyObject *kwargs) +{ + PyObject *value; + + if (!PyArg_ParseTuple(args, "O", &value)) { + return NULL; + } + + return (PyObject *)QuadPrecision_from_object(value); +} + +static PyObject * QuadPrecision_str(QuadPrecisionObject * self) +{ + char buffer[128]; + Sleef_snprintf(buffer, sizeof(buffer), "%.*Qe", SLEEF_QUAD_DIG, self->quad.value); + return PyUnicode_FromString(buffer); +} + +static PyObject * QuadPrecision_repr(QuadPrecisionObject* self) +{ + PyObject *str = QuadPrecision_str(self); + if (str == NULL) { + return NULL; + } + PyObject *res = PyUnicode_FromFormat("QuadPrecision('%S')", str); + Py_DECREF(str); + return res; +} + +static void +quad_dealloc(QuadPrecDTypeObject *self) +{ + Py_TYPE(self)->tp_free((PyObject *)self); +} + +PyTypeObject QuadPrecision_Type = +{ + PyVarObject_HEAD_INIT(NULL, 0) + .tp_name = "QuadPrecType.QuadPrecision", + .tp_basicsize = sizeof(QuadPrecisionObject), + .tp_itemsize = 0, + .tp_new = QuadPrecision_new, + .tp_dealloc = (destructor)quad_dealloc, + .tp_repr = (reprfunc)QuadPrecision_repr, + .tp_str = (reprfunc)QuadPrecision_str, + .tp_as_number = &quad_as_scalar, + .tp_richcompare = (richcmpfunc)quad_richcompare + +}; + +int +init_quadprecision_scalar(void) +{ + return PyType_Ready(&QuadPrecision_Type); +} \ No newline at end of file diff --git a/quaddtype/quaddtype/src/scalar.h b/quaddtype/quaddtype/src/scalar.h new file mode 100644 index 00000000..962a1fec --- /dev/null +++ b/quaddtype/quaddtype/src/scalar.h @@ -0,0 +1,35 @@ +#ifndef _QUADDTYPE_SCALAR_H +#define _QUADDTYPE_SCALAR_H + +#ifdef __cplusplus +extern "C" { +#endif + +#include +#include + +typedef struct { + Sleef_quad value; +} quad_field; + +typedef struct { + PyObject_HEAD + quad_field quad; +} QuadPrecisionObject; + +extern PyTypeObject QuadPrecision_Type; + +QuadPrecisionObject * +QuadPrecision_raw_new(void); + +QuadPrecisionObject * +QuadPrecision_from_object(PyObject *value); + +int +init_quadprecision_scalar(void); + +#ifdef __cplusplus +} +#endif + +#endif \ No newline at end of file diff --git a/quaddtype/quaddtype/src/scalar_ops.cpp b/quaddtype/quaddtype/src/scalar_ops.cpp new file mode 100644 index 00000000..a51bc69a --- /dev/null +++ b/quaddtype/quaddtype/src/scalar_ops.cpp @@ -0,0 +1,171 @@ +#define PY_ARRAY_UNIQUE_SYMBOL QuadPrecType_ARRAY_API +#define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION +#define NPY_TARGET_VERSION NPY_2_0_API_VERSION +#define NO_IMPORT_ARRAY + +extern "C" { +#include + +#include "numpy/arrayobject.h" +#include "numpy/ndarraytypes.h" +#include "numpy/ufuncobject.h" + +#include "numpy/dtype_api.h" +} + +#include "scalar.h" +#include "ops.hpp" +#include "scalar_ops.h" + +template +static PyObject * +quad_unary_func(QuadPrecisionObject *self) +{ + QuadPrecisionObject *res = QuadPrecision_raw_new(); + if (!res) { + return NULL; + } + + unary_op(&self->quad.value, &res->quad.value); + return (PyObject *)res; +} + +PyObject * +quad_nonzero(QuadPrecisionObject *self) +{ + return PyBool_FromLong(Sleef_icmpneq1(self->quad.value, Sleef_cast_from_int64q1(0))); +} + +template +static PyObject * +quad_binary_func(PyObject *op1, PyObject *op2) +{ + QuadPrecisionObject *self; + PyObject *other; + QuadPrecisionObject *other_quad = NULL; + int is_forward; + + if (PyObject_TypeCheck(op1, &QuadPrecision_Type)) { + is_forward = 1; + self = (QuadPrecisionObject *)op1; + other = Py_NewRef(op2); + } + else { + is_forward = 0; + self = (QuadPrecisionObject *)op2; + other = Py_NewRef(op1); + } + + if (PyObject_TypeCheck(other, &QuadPrecision_Type)) { + Py_INCREF(other); + other_quad = (QuadPrecisionObject *)other; + } + else if (PyLong_Check(other) || PyFloat_Check(other)) { + other_quad = QuadPrecision_raw_new(); + if (!other_quad) { + Py_DECREF(other); + return NULL; + } + + if (PyLong_Check(other)) { + long long value = PyLong_AsLongLong(other); + if (value == -1 && PyErr_Occurred()) { + Py_DECREF(other); + Py_DECREF(other_quad); + return NULL; + } + other_quad->quad.value = Sleef_cast_from_int64q1(value); + } + else { + double value = PyFloat_AsDouble(other); + if (value == -1.0 && PyErr_Occurred()) { + Py_DECREF(other); + Py_DECREF(other_quad); + return NULL; + } + other_quad->quad.value = Sleef_cast_from_doubleq1(value); + } + } + else { + Py_DECREF(other); + Py_RETURN_NOTIMPLEMENTED; + } + + QuadPrecisionObject *res = QuadPrecision_raw_new(); + if (!res) { + Py_DECREF(other_quad); + Py_DECREF(other); + return NULL; + } + + if (is_forward) { + binary_op(&res->quad.value, &self->quad.value, &other_quad->quad.value); + } + else { + binary_op(&res->quad.value, &other_quad->quad.value, &self->quad.value); + } + + Py_DECREF(other_quad); + Py_DECREF(other); + return (PyObject *)res; +} + +// todo: add support with float and int +PyObject * +quad_richcompare(QuadPrecisionObject *self, PyObject *other, int cmp_op) +{ + QuadPrecisionObject *other_quad = NULL; + + if (PyObject_TypeCheck(other, &QuadPrecision_Type)) { + Py_INCREF(other); + other_quad = (QuadPrecisionObject *)other; + } + else if (PyLong_CheckExact(other) || PyFloat_CheckExact(other)) { + other_quad = QuadPrecision_from_object(other); + if (other_quad == NULL) { + return NULL; + } + } + else { + Py_RETURN_NOTIMPLEMENTED; + } + int cmp; + switch (cmp_op) { + case Py_LT: + cmp = Sleef_icmpltq1(self->quad.value, other_quad->quad.value); + break; + case Py_LE: + cmp = Sleef_icmpleq1(self->quad.value, other_quad->quad.value); + break; + case Py_EQ: + cmp = Sleef_icmpeqq1(self->quad.value, other_quad->quad.value); + break; + case Py_NE: + cmp = Sleef_icmpneq1(self->quad.value, other_quad->quad.value); + break; + case Py_GT: + cmp = Sleef_icmpgtq1(self->quad.value, other_quad->quad.value); + break; + case Py_GE: + cmp = Sleef_icmpgeq1(self->quad.value, other_quad->quad.value); + break; + default: + Py_DECREF(other_quad); + Py_RETURN_NOTIMPLEMENTED; + } + Py_DECREF(other_quad); + + return PyBool_FromLong(cmp); +} + +PyNumberMethods quad_as_scalar = { + .nb_add = (binaryfunc)quad_binary_func, + .nb_subtract = (binaryfunc)quad_binary_func, + .nb_multiply = (binaryfunc)quad_binary_func, + .nb_power = (ternaryfunc)quad_binary_func, + .nb_negative = (unaryfunc)quad_unary_func, + .nb_positive = (unaryfunc)quad_unary_func, + .nb_absolute = (unaryfunc)quad_unary_func, + .nb_bool = (inquiry)quad_nonzero, + .nb_true_divide = (binaryfunc)quad_binary_func, +}; \ No newline at end of file diff --git a/quaddtype/quaddtype/src/scalar_ops.h b/quaddtype/quaddtype/src/scalar_ops.h new file mode 100644 index 00000000..138738fd --- /dev/null +++ b/quaddtype/quaddtype/src/scalar_ops.h @@ -0,0 +1,20 @@ +#ifndef _QUADDTYPE_SCALAR_OPS_H +#define _QUADDTYPE_SCALAR_OPS_H + +#ifdef __cplusplus +extern "C" { +#endif + +#include +#include "dtype.h" + +PyObject * +quad_richcompare(QuadPrecisionObject *self, PyObject *other, int cmp_op); + +extern PyNumberMethods quad_as_scalar; + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/quaddtype/quaddtype/src/umath.c b/quaddtype/quaddtype/src/umath.c deleted file mode 100644 index 6c9fece6..00000000 --- a/quaddtype/quaddtype/src/umath.c +++ /dev/null @@ -1,124 +0,0 @@ -#include - -#define PY_ARRAY_UNIQUE_SYMBOL quaddtype_ARRAY_API -#define PY_UFUNC_UNIQUE_SYMBOL quaddtype_UFUNC_API -#define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION -#define NPY_TARGET_VERSION NPY_2_0_API_VERSION -#define NO_IMPORT_ARRAY -#define NO_IMPORT_UFUNC -#include "numpy/ndarraytypes.h" -#include "numpy/arrayobject.h" -#include "numpy/ufuncobject.h" -#include "numpy/dtype_api.h" - -#include "dtype.h" -#include "umath.h" - -// The multiplication loop, this is very minimal! Look at the cast to see -// some of the more advanced things you can do for optimization! -// Look at the seberg/unitdtype exaple repository for how this can be done -// more generically without implementing a multiply loop! -static int -quad_multiply_strided_loop(PyArrayMethod_Context *context, char *const data[], - npy_intp const dimensions[], npy_intp const strides[], - NpyAuxData *auxdata) -{ - npy_intp N = dimensions[0]; - char *in1 = data[0], *in2 = data[1]; - char *out = data[2]; - npy_intp in1_stride = strides[0]; - npy_intp in2_stride = strides[1]; - npy_intp out_stride = strides[2]; - - while (N--) { - *(__float128 *)out = *(__float128 *)in1 * *(__float128 *)in2; - in1 += in1_stride; - in2 += in2_stride; - out += out_stride; - } - return 0; -} - -// This is the "dtype/descriptor resolver". Its main job is to fill in the -// result dtype in this case, that is: - -// given_descrs[0] = arr1.dtype -// given_descrs[1] = arr2.dtype -// given_descrs[2] = NULL or out.dtype - -// The code now sets `loop_descrs` to what we actually get as a result. -// In practice that is just fill in the last `out.dtype`. But in principle, -// we might need casting (e.g. swap `>f8` to ` +#include + +#include "numpy/arrayobject.h" +#include "numpy/ndarraytypes.h" +#include "numpy/ufuncobject.h" + +#include "numpy/dtype_api.h" +} +#include "dtype.h" +#include "umath.h" +#include "ops.hpp" + +template +int +quad_generic_unary_op_strided_loop(PyArrayMethod_Context *context, char *const data[], + npy_intp const dimensions[], npy_intp const strides[], + NpyAuxData *auxdata) +{ + npy_intp N = dimensions[0]; + char *in_ptr = data[0]; + char *out_ptr = data[1]; + npy_intp in_stride = strides[0]; + npy_intp out_stride = strides[1]; + + Sleef_quad in, out; + while (N--) { + memcpy(&in, in_ptr, sizeof(Sleef_quad)); + unary_op(&in, &out); + memcpy(out_ptr, &out, sizeof(Sleef_quad)); + + in_ptr += in_stride; + out_ptr += out_stride; + } + return 0; +} + +static NPY_CASTING +quad_unary_op_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtypes[], + PyArray_Descr *const given_descrs[], PyArray_Descr *loop_descrs[], + npy_intp *NPY_UNUSED(view_offset)) +{ + Py_INCREF(given_descrs[0]); + loop_descrs[0] = given_descrs[0]; + + if (given_descrs[1] == NULL) { + Py_INCREF(given_descrs[0]); + loop_descrs[1] = given_descrs[0]; + return NPY_NO_CASTING; + } + Py_INCREF(given_descrs[1]); + loop_descrs[1] = given_descrs[1]; + + return NPY_NO_CASTING; +} + +template +int +create_quad_unary_ufunc(PyObject *numpy, const char *ufunc_name) +{ + PyObject *ufunc = PyObject_GetAttrString(numpy, ufunc_name); + if (ufunc == NULL) { + return -1; + } + + PyArray_DTypeMeta *dtypes[2] = {&QuadPrecDType, &QuadPrecDType}; + + PyType_Slot slots[] = { + {NPY_METH_resolve_descriptors, (void *)&quad_unary_op_resolve_descriptors}, + {NPY_METH_strided_loop, (void *)&quad_generic_unary_op_strided_loop}, + {0, NULL}}; + + PyArrayMethod_Spec Spec = { + .name = "quad_unary_op", + .nin = 1, + .nout = 1, + .casting = NPY_NO_CASTING, + .flags = (NPY_ARRAYMETHOD_FLAGS)0, + .dtypes = dtypes, + .slots = slots, + }; + + if (PyUFunc_AddLoopFromSpec(ufunc, &Spec) < 0) { + return -1; + } + + return 0; +} + +int +init_quad_unary_ops(PyObject *numpy) +{ + if (create_quad_unary_ufunc(numpy, "negative") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "absolute") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "rint") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "trunc") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "floor") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "ceil") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "sqrt") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "square") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "log") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "log2") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "log10") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "log1p") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "exp") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "exp2") < 0) { + return -1; + } + return 0; +} + +// Binary ufuncs + +template +int +quad_generic_binop_strided_loop(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], *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]; + + Sleef_quad in1, in2, out; + while (N--) { + memcpy(&in1, in1_ptr, sizeof(Sleef_quad)); + memcpy(&in2, in2_ptr, sizeof(Sleef_quad)); + binop(&out, &in1, &in2); + memcpy(out_ptr, &out, sizeof(Sleef_quad)); + + in1_ptr += in1_stride; + in2_ptr += in2_stride; + out_ptr += out_stride; + } + return 0; +} + +static NPY_CASTING +quad_binary_op_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtypes[], + PyArray_Descr *const given_descrs[], + PyArray_Descr *loop_descrs[], npy_intp *NPY_UNUSED(view_offset)) +{ + Py_INCREF(given_descrs[0]); + loop_descrs[0] = given_descrs[0]; + Py_INCREF(given_descrs[1]); + loop_descrs[1] = given_descrs[1]; + + if (given_descrs[2] == NULL) { + PyArray_Descr *out_descr = (PyArray_Descr *)new_quaddtype_instance(); + if (!out_descr) { + return (NPY_CASTING)-1; + } + Py_INCREF(given_descrs[0]); + loop_descrs[2] = out_descr; + } + else { + Py_INCREF(given_descrs[2]); + loop_descrs[2] = given_descrs[2]; + } + + return NPY_NO_CASTING; +} + +// helper debugging function +static const char * +get_dtype_name(PyArray_DTypeMeta *dtype) +{ + if (dtype == &QuadPrecDType) { + return "QuadPrecDType"; + } + else if (dtype == &PyArray_BoolDType) { + return "BoolDType"; + } + else if (dtype == &PyArray_ByteDType) { + return "ByteDType"; + } + else if (dtype == &PyArray_UByteDType) { + return "UByteDType"; + } + else if (dtype == &PyArray_ShortDType) { + return "ShortDType"; + } + else if (dtype == &PyArray_UShortDType) { + return "UShortDType"; + } + else if (dtype == &PyArray_IntDType) { + return "IntDType"; + } + else if (dtype == &PyArray_UIntDType) { + return "UIntDType"; + } + else if (dtype == &PyArray_LongDType) { + return "LongDType"; + } + else if (dtype == &PyArray_ULongDType) { + return "ULongDType"; + } + else if (dtype == &PyArray_LongLongDType) { + return "LongLongDType"; + } + else if (dtype == &PyArray_ULongLongDType) { + return "ULongLongDType"; + } + else if (dtype == &PyArray_FloatDType) { + return "FloatDType"; + } + else if (dtype == &PyArray_DoubleDType) { + return "DoubleDType"; + } + else if (dtype == &PyArray_LongDoubleDType) { + return "LongDoubleDType"; + } + else { + return "UnknownDType"; + } +} + +static int +quad_ufunc_promoter(PyUFuncObject *ufunc, PyArray_DTypeMeta *op_dtypes[], + PyArray_DTypeMeta *signature[], PyArray_DTypeMeta *new_op_dtypes[]) +{ + // printf("quad_ufunc_promoter called for ufunc: %s\n", ufunc->name); + // printf("Entering quad_ufunc_promoter\n"); + // printf("Ufunc name: %s\n", ufunc->name); + // printf("nin: %d, nargs: %d\n", ufunc->nin, ufunc->nargs); + + int nin = ufunc->nin; + int nargs = ufunc->nargs; + PyArray_DTypeMeta *common = NULL; + bool has_quad = false; + + // Handle the special case for reductions + if (op_dtypes[0] == NULL) { + assert(nin == 2 && ufunc->nout == 1); /* must be reduction */ + for (int i = 0; i < 3; i++) { + Py_INCREF(op_dtypes[1]); + new_op_dtypes[i] = op_dtypes[1]; + // printf("new_op_dtypes[%d] set to %s\n", i, get_dtype_name(new_op_dtypes[i])); + } + return 0; + } + + // Check if any input or signature is QuadPrecision + for (int i = 0; i < nargs; i++) { + if ((i < nin && op_dtypes[i] == &QuadPrecDType) || (signature[i] == &QuadPrecDType)) { + has_quad = true; + // printf("QuadPrecision detected in input %d or signature\n", i); + break; + } + } + + if (has_quad) { + // If QuadPrecision is involved, use it for all arguments + common = &QuadPrecDType; + // printf("Using QuadPrecDType as common type\n"); + } + else { + // Check if output signature is homogeneous + for (int i = nin; i < nargs; i++) { + if (signature[i] != NULL) { + if (common == NULL) { + Py_INCREF(signature[i]); + common = signature[i]; + // printf("Common type set to %s from signature\n", get_dtype_name(common)); + } + else if (common != signature[i]) { + Py_CLEAR(common); // Not homogeneous, unset common + // printf("Output signature not homogeneous, cleared common type\n"); + break; + } + } + } + + // If no common output dtype, use standard promotion for inputs + if (common == NULL) { + // printf("Using standard promotion for inputs\n"); + common = PyArray_PromoteDTypeSequence(nin, op_dtypes); + if (common == NULL) { + if (PyErr_ExceptionMatches(PyExc_TypeError)) { + PyErr_Clear(); // Do not propagate normal promotion errors + } + // printf("Exiting quad_ufunc_promoter (promotion failed)\n"); + return -1; + } + // printf("Common type after promotion: %s\n", get_dtype_name(common)); + } + } + + // Set all new_op_dtypes to the common dtype + for (int i = 0; i < nargs; i++) { + if (signature[i]) { + // If signature is specified for this argument, use it + Py_INCREF(signature[i]); + new_op_dtypes[i] = signature[i]; + // printf("new_op_dtypes[%d] set to %s (from signature)\n", i, + // get_dtype_name(new_op_dtypes[i])); + } + else { + // Otherwise, use the common dtype + Py_INCREF(common); + new_op_dtypes[i] = common; + // printf("new_op_dtypes[%d] set to %s (from common)\n", i, + // get_dtype_name(new_op_dtypes[i])); + } + } + + Py_XDECREF(common); + // printf("Exiting quad_ufunc_promoter\n"); + return 0; +} + +template +int +create_quad_binary_ufunc(PyObject *numpy, const char *ufunc_name) +{ + PyObject *ufunc = PyObject_GetAttrString(numpy, ufunc_name); + if (ufunc == NULL) { + Py_DecRef(ufunc); + return -1; + } + + PyArray_DTypeMeta *dtypes[3] = {&QuadPrecDType, &QuadPrecDType, &QuadPrecDType}; + + PyType_Slot slots[] = { + {NPY_METH_resolve_descriptors, (void *)&quad_binary_op_resolve_descriptors}, + {NPY_METH_strided_loop, (void *)&quad_generic_binop_strided_loop}, + {0, NULL}}; + + PyArrayMethod_Spec Spec = { + .name = "quad_binop", + .nin = 2, + .nout = 1, + .casting = NPY_NO_CASTING, + .flags = (NPY_ARRAYMETHOD_FLAGS)0, + .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, &PyArrayDescr_Type, &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; +} + +int +init_quad_binary_ops(PyObject *numpy) +{ + if (create_quad_binary_ufunc(numpy, "add") < 0) { + return -1; + } + if (create_quad_binary_ufunc(numpy, "subtract") < 0) { + return -1; + } + if (create_quad_binary_ufunc(numpy, "multiply") < 0) { + return -1; + } + if (create_quad_binary_ufunc(numpy, "divide") < 0) { + return -1; + } + if (create_quad_binary_ufunc(numpy, "power") < 0) { + return -1; + } + if (create_quad_binary_ufunc(numpy, "mod") < 0) { + return -1; + } + if (create_quad_binary_ufunc(numpy, "minimum") < 0) { + return -1; + } + if (create_quad_binary_ufunc(numpy, "maximum") < 0) { + return -1; + } + return 0; +} + +// comparison functions + +template +int +quad_generic_comp_strided_loop(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], *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]; + + while (N--) { + *((npy_bool *)out_ptr) = comp((Sleef_quad *)in1_ptr, (Sleef_quad *)in2_ptr); + + in1_ptr += in1_stride; + in2_ptr += in2_stride; + out_ptr += out_stride; + } + return 0; +} + +NPY_NO_EXPORT int +comparison_ufunc_promoter(PyUFuncObject *ufunc, PyArray_DTypeMeta *op_dtypes[], + PyArray_DTypeMeta *signature[], PyArray_DTypeMeta *new_op_dtypes[]) +{ + PyArray_DTypeMeta *new_signature[NPY_MAXARGS]; + + memcpy(new_signature, signature, 3 * sizeof(PyArray_DTypeMeta *)); + new_signature[2] = NULL; + int res = quad_ufunc_promoter(ufunc, op_dtypes, new_signature, new_op_dtypes); + if (res < 0) { + return -1; + } + Py_XSETREF(new_op_dtypes[2], &PyArray_BoolDType); + return 0; +} + +template +int +create_quad_comparison_ufunc(PyObject *numpy, const char *ufunc_name) +{ + PyObject *ufunc = PyObject_GetAttrString(numpy, ufunc_name); + if (ufunc == NULL) { + return -1; + } + + PyArray_DTypeMeta *dtypes[3] = {&QuadPrecDType, &QuadPrecDType, &PyArray_BoolDType}; + + PyType_Slot slots[] = {{NPY_METH_strided_loop, (void *)&quad_generic_comp_strided_loop}, + {0, NULL}}; + + PyArrayMethod_Spec Spec = { + .name = "quad_comp", + .nin = 2, + .nout = 1, + .casting = NPY_NO_CASTING, + .flags = (NPY_ARRAYMETHOD_FLAGS)0, + .dtypes = dtypes, + .slots = slots, + }; + + if (PyUFunc_AddLoopFromSpec(ufunc, &Spec) < 0) { + return -1; + } + + PyObject *promoter_capsule = + PyCapsule_New((void *)&comparison_ufunc_promoter, "numpy._ufunc_promoter", NULL); + if (promoter_capsule == NULL) { + return -1; + } + + PyObject *DTypes = PyTuple_Pack(3, &PyArrayDescr_Type, &PyArrayDescr_Type, &PyArray_BoolDType); + 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; +} + +int +init_quad_comps(PyObject *numpy) +{ + if (create_quad_comparison_ufunc(numpy, "equal") < 0) { + return -1; + } + if (create_quad_comparison_ufunc(numpy, "not_equal") < 0) { + return -1; + } + if (create_quad_comparison_ufunc(numpy, "less") < 0) { + return -1; + } + if (create_quad_comparison_ufunc(numpy, "less_equal") < 0) { + return -1; + } + if (create_quad_comparison_ufunc(numpy, "greater") < 0) { + return -1; + } + if (create_quad_comparison_ufunc(numpy, "greater_equal") < 0) { + return -1; + } + + return 0; +} + +int +init_quad_umath(void) +{ + PyObject *numpy = PyImport_ImportModule("numpy"); + if (!numpy) + return -1; + + if (init_quad_unary_ops(numpy) < 0) { + goto err; + } + + if (init_quad_binary_ops(numpy) < 0) { + goto err; + } + + if (init_quad_comps(numpy) < 0) { + goto err; + } + + Py_DECREF(numpy); + return 0; + +err: + Py_DECREF(numpy); + return -1; +} \ No newline at end of file diff --git a/quaddtype/quaddtype/src/umath.h b/quaddtype/quaddtype/src/umath.h index b883f57a..d64f26be 100644 --- a/quaddtype/quaddtype/src/umath.h +++ b/quaddtype/quaddtype/src/umath.h @@ -1,7 +1,15 @@ -#ifndef _NPY_UFUNC_H -#define _NPY_UFUNC_H +#ifndef _QUADDTYPE_UMATH_H +#define _QUADDTYPE_UMATH_H + +#ifdef __cplusplus +extern "C" { +#endif int -init_multiply_ufunc(void); +init_quad_umath(void); + +#ifdef __cplusplus +} +#endif -#endif /*_NPY_UFUNC_H */ +#endif diff --git a/quaddtype/reinstall.sh b/quaddtype/reinstall.sh index 52125765..31479511 100755 --- a/quaddtype/reinstall.sh +++ b/quaddtype/reinstall.sh @@ -8,7 +8,6 @@ then fi #meson setup build -Db_sanitize=address,undefined -meson setup build python -m pip uninstall -y quaddtype python -m pip install . -v --no-build-isolation -Cbuilddir=build -C'compile-args=-v' -Csetup-args="-Dbuildtype=debug" -#python -m pip install . -v --no-build-isolation -Cbuilddir=build -C'compile-args=-v' +#python -m pip install . -v --no-build-isolation -Cbuilddir=build -C'compile-args=-v' \ No newline at end of file diff --git a/quaddtype/tests/conftest.py b/quaddtype/tests/conftest.py deleted file mode 100644 index f5870e45..00000000 --- a/quaddtype/tests/conftest.py +++ /dev/null @@ -1,3 +0,0 @@ -import os - -os.environ["NUMPY_EXPERIMENTAL_DTYPE_API"] = "1" diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index fc65b232..20ccd91c 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -1,33 +1,108 @@ +import pytest +import sys import numpy as np +import operator -from quaddtype import QuadDType, QuadScalar +from quaddtype import QuadPrecDType, QuadPrecision -def test_dtype_creation(): - assert str(QuadDType()) == "This is a quad (128-bit float) dtype." +def test_create_scalar_simple(): + assert isinstance(QuadPrecision("12.0"), QuadPrecision) + assert isinstance(QuadPrecision(1.63), QuadPrecision) + assert isinstance(QuadPrecision(1), QuadPrecision) + + +def test_basic_equality(): + assert QuadPrecision("12") == QuadPrecision( + "12.0") == QuadPrecision("12.00") + + +@pytest.mark.parametrize("val", ["123532.543", "12893283.5"]) +def test_scalar_repr(val): + expected = f"QuadPrecision('{str(QuadPrecision(val))}')" + assert repr(QuadPrecision(val)) == expected + + +@pytest.mark.parametrize("op", ["add", "sub", "mul", "truediv", "pow"]) +@pytest.mark.parametrize("other", ["3.0", "12.5", "100.0"]) +def test_binary_ops(op, other): + op_func = getattr(operator, op) + quad_a = QuadPrecision("12.5") + quad_b = QuadPrecision(other) + float_a = 12.5 + float_b = float(other) + + quad_result = op_func(quad_a, quad_b) + float_result = op_func(float_a, float_b) + assert np.abs(np.float64(quad_result) - float_result) < 1e-10 -def test_scalar_creation(): - assert str(QuadScalar(3.1)) == "3.1" +@pytest.mark.parametrize("op", ["eq", "ne", "le", "lt", "ge", "gt"]) +@pytest.mark.parametrize("other", ["3.0", "12.5", "100.0"]) +def test_comparisons(op, other): + op_func = getattr(operator, op) + quad_a = QuadPrecision("12.5") + quad_b = QuadPrecision(other) + float_a = 12.5 + float_b = float(other) -def test_create_with_explicit_dtype(): - assert ( - repr(np.array([3.0, 3.1, 3.2], dtype=QuadDType())) - == "array([3.0, 3.1, 3.2], dtype=This is a quad (128-bit float) dtype.)" - ) + assert op_func(quad_a, quad_b) == op_func(float_a, float_b) -def test_multiply(): - x = np.array([3, 8.0], dtype=QuadDType()) - assert str(x * x) == "[9.0 64.0]" +@pytest.mark.parametrize("op, val, expected", [ + ("neg", "3.0", "-3.0"), + ("neg", "-3.0", "3.0"), + ("pos", "3.0", "3.0"), + ("pos", "-3.0", "-3.0"), + ("abs", "3.0", "3.0"), + ("abs", "-3.0", "3.0"), + ("neg", "12.5", "-12.5"), + ("pos", "100.0", "100.0"), + ("abs", "-25.5", "25.5"), +]) +def test_unary_ops(op, val, expected): + quad_val = QuadPrecision(val) + expected_val = QuadPrecision(expected) + + if op == "neg": + result = -quad_val + elif op == "pos": + result = +quad_val + elif op == "abs": + result = abs(quad_val) + else: + raise ValueError(f"Unsupported operation: {op}") + + assert result == expected_val, f"{op}({val}) should be {expected}, but got {result}" + + +def test_nan_and_inf(): + assert (QuadPrecision("nan") != QuadPrecision("nan")) == ( + QuadPrecision("nan") == QuadPrecision("nan")) + assert QuadPrecision("inf") > QuadPrecision("1e1000") + assert QuadPrecision("-inf") < QuadPrecision("-1e1000") + + +def test_dtype_creation(): + dtype = QuadPrecDType() + assert isinstance(dtype, np.dtype) + assert dtype.name == 'QuadPrecDType128' + +def test_array_creation(): + arr = np.array([1, 2, 3], dtype=QuadPrecDType()) + assert arr.dtype.name == 'QuadPrecDType128' + assert all(isinstance(x, QuadPrecision) for x in arr) -def test_bytes(): - """Check that each quad is 16 bytes.""" - x = np.array([3, 8.0, 1.4], dtype=QuadDType()) - assert len(x.tobytes()) == x.size * 16 +def test_array_operations(): + arr1 = np.array( + [QuadPrecision("1.5"), QuadPrecision("2.5"), QuadPrecision("3.5")]) + arr2 = np.array( + [QuadPrecision("0.5"), QuadPrecision("1.0"), QuadPrecision("1.5")]) -def test_is_numeric(): - assert QuadDType._is_numeric + result = arr1 + arr2 + expected = np.array( + [QuadPrecision("2.0"), QuadPrecision("3.5"), QuadPrecision("5.0")]) + assert np.all(result == expected)