API: Raise when mixing SLEEF and longdouble backends

.github/workflows/big_endian.yml

@@ -137,7 +137,8 @@ jobs:
               python -m pip install --break-system-packages --no-deps . -v --no-build-isolation --force-reinstall &&
 
               # Install test dependencies separately
-              python -m pip install --break-system-packages .[test] &&
+              # pandas 3.1.0 fails to build on s390x: pandas-dev/pandas#71548
+              python -m pip install --break-system-packages .[test] pandas!=3.1.0 &&
 
               python -m pytest -vvv --color=yes --timeout=600 --tb=short tests/
             '"

README.md

@@ -26,6 +26,7 @@ A cross-platform Quad (128-bit) float Data-Type for NumPy.
     - [Serving the documentation](#serving-the-documentation)
   - [Development Tips](#development-tips)
     - [Cleaning the Build Directory](#cleaning-the-build-directory)
+  - [Changelog](#changelog)
 
 ## Installation
 
@@ -257,3 +258,46 @@ The subproject folders (`subprojects/sleef`, `subprojects/qblas`) are cloned as
 ```bash
 git clean -ffxd
 ```
+
+## Changelog
+
+### Unreleased
+
+#### Backwards-incompatible changes
+
+- Values no longer move between the SLEEF and longdouble backends. Every
+  operation that would combine or convert them now raises `TypeError`, with
+  the message `Cannot mix QuadPrecision backends 'sleef' and 'longdouble'`,
+  or NumPy's own `Cannot cast ...` message where NumPy checks the casting
+  rule first (for example `astype(..., casting="same_kind")`). Previously
+  these silently converted through `float64` (losing precision) or picked a
+  backend inconsistently. This covers:
+  - ufuncs with operands or `out=` arrays of different backends (before:
+    arithmetic and comparisons silently used longdouble);
+  - `astype`, `np.copyto`, slice and item assignment between backends
+    (before: silent conversion through `float64`);
+  - `np.result_type`, `np.promote_types` and `np.concatenate` across
+    backends (before: silently picked SLEEF); `np.array` of mixed
+    `QuadPrecision` scalars now gives an object array;
+  - `QuadPrecision(value, backend=...)` with a `QuadPrecision` of the other
+    backend, including `QuadPrecision(x)` for a longdouble `x`, since the
+    default backend is SLEEF.
+- Python `int` and `float` operands use the SLEEF backend, so with longdouble
+  arrays they now raise as well. This includes NumPy functions that pass
+  Python scalars internally, such as `np.ones`, `np.full`, `np.ones_like`,
+  `mean`, `var`, `std`, `np.median`, `np.isclose`, `np.allclose` and
+  `np.round`. Use NumPy scalars or arrays
+  instead, for example `ld / np.float64(2)`,
+  `np.full(3, np.longdouble(1), dtype=ld.dtype)` or
+  `np.ones(3).astype(ld.dtype)`. SLEEF arrays and NumPy's own dtypes,
+  including `np.longdouble`, are unaffected.
+
+#### Bug fixes
+
+- Assigning a `QuadPrecision` scalar into an array of the other backend, or
+  casting an object array of them, stored the raw bytes under the wrong
+  backend (for example `0.0` instead of `1.5`). It now raises.
+- Unary ufuncs such as `np.sqrt` with an `out=` array of the other backend
+  and `casting="unsafe"` wrote values in the wrong format. They now raise.
+- `QuadPrecision(q)` for a `QuadPrecision` `q` of the same backend rounded
+  through `float64`; it now copies the value exactly.

docs/api/core.md

@@ -136,6 +136,39 @@ The fundamental types provided by NumPy QuadDType.
        True
 ```
 
+### Mixing backends
+
+Values never move between the SLEEF and longdouble backends: any operation
+that would combine or convert them raises `TypeError`: either
+`Cannot mix QuadPrecision backends ...` or, where NumPy checks the casting rule
+first, NumPy's own `Cannot cast ...` message. This covers ufuncs,
+`out=` arrays, `astype`, assignment, `np.concatenate`, `np.result_type` and
+`QuadPrecision(value, backend=...)`. NumPy's own dtypes, including
+`np.longdouble`, work with both backends.
+
+Python `int` and `float` operands use the SLEEF backend, so with longdouble
+arrays they raise too. This includes NumPy functions that pass Python scalars
+internally, such as `np.ones`, `np.full`, `np.ones_like`, `mean`, `var`,
+`std`, `np.median`, `np.isclose`, `np.allclose` and `np.round`. Pass NumPy scalars or arrays, or
+`QuadPrecision(..., backend="longdouble")`, instead, for example
+`np.full(3, np.longdouble(1), dtype=dt)` or `np.ones(3).astype(dt)`:
+
+```python
+>>> import numpy as np
+>>> from numpy_quaddtype import QuadPrecDType
+>>> ld = np.array([1.5, 2.5], dtype=QuadPrecDType(backend="longdouble"))
+>>> ld + np.float64(1)
+array([2.5, 3.5], dtype=QuadPrecDType(backend='longdouble'))
+>>> ld + 1
+Traceback (most recent call last):
+  ...
+TypeError: Cannot mix QuadPrecision backends 'sleef' and 'longdouble'
+>>> ld.astype(QuadPrecDType())
+Traceback (most recent call last):
+  ...
+TypeError: Cannot mix QuadPrecision backends 'sleef' and 'longdouble'
+```
+
 ## Convenience Functions
 
 ```{eval-rst}

src/csrc/casts.cpp

@@ -73,14 +73,9 @@ quad_to_quad_resolve_descriptors(PyObject *NPY_UNUSED(self),
     loop_descrs[1] = given_descrs[1];
 
     if (given_descrs[0]->backend != given_descrs[1]->backend) {
-        // Different backends require actual conversion, no view possible
-        *view_offset = NPY_MIN_INTP;
-        if (given_descrs[0]->backend == BACKEND_SLEEF) {
-            // SLEEF -> long double may lose precision
-            return static_cast<NPY_CASTING>(NPY_SAME_KIND_CASTING | NPY_SAME_VALUE_CASTING_FLAG);
-        }
-        // long double -> SLEEF preserves value exactly
-        return static_cast<NPY_CASTING>(NPY_SAFE_CASTING | NPY_SAME_VALUE_CASTING_FLAG);
+        // gh-112: the loop already handles cross-backend casts; allow them here once they are exact.
+        quad_set_backend_mismatch_error();
+        return quad_resolve_descrs_fail((PyArray_Descr **)loop_descrs, 2);
     }
 
     *view_offset = 0;

src/csrc/dtype.c

@@ -55,15 +55,9 @@ ensure_canonical(QuadPrecDTypeObject *self)
 static QuadPrecDTypeObject *
 common_instance(QuadPrecDTypeObject *dtype1, QuadPrecDTypeObject *dtype2)
 {
-    // if backend mismatch then return SLEEF one (safe to cast ld to quad)
     if (dtype1->backend != dtype2->backend) {
-        if (dtype1->backend == BACKEND_SLEEF) {
-            Py_INCREF(dtype1);
-            return dtype1;
-        }
-
-        Py_INCREF(dtype2);
-        return dtype2;
+        quad_set_backend_mismatch_error();
+        return NULL;
     }
     Py_INCREF(dtype1);
     return dtype1;
@@ -120,6 +114,10 @@ quadprec_setitem(QuadPrecDTypeObject *descr, PyObject *obj, char *dataptr)
 {
     QuadPrecisionObject *value;
     if (PyObject_TypeCheck(obj, &QuadPrecision_Type)) {
+        if (((QuadPrecisionObject *)obj)->backend != descr->backend) {
+            quad_set_backend_mismatch_error();
+            return -1;
+        }
         Py_INCREF(obj);
         value = (QuadPrecisionObject *)obj;
     }

src/csrc/scalar.c

@@ -96,7 +96,23 @@ quad_from_py_int(PyObject *py_int, QuadBackendType backend, QuadPrecisionObject
 
 QuadPrecisionObject *
 QuadPrecision_from_object(PyObject *value, QuadBackendType backend)
-{   
+{
+    // QuadPrecision is a NumPy floating scalar, so it must be handled before the
+    // NumPy scalar path, which converts through double.
+    if (PyObject_TypeCheck(value, &QuadPrecision_Type)) {
+        QuadPrecisionObject *quad_obj = (QuadPrecisionObject *)value;
+        if (quad_obj->backend != backend) {
+            quad_set_backend_mismatch_error();
+            return NULL;
+        }
+        QuadPrecisionObject *self = QuadPrecision_raw_new(backend);
+        if (self == NULL) {
+            return NULL;
+        }
+        self->value = quad_obj->value;
+        return self;
+    }
+
     // Handle numpy scalars (np.int32, np.float32, etc.) before arrays
     // We need to check this before PySequence_Check because some numpy scalars are sequences
     if (PyArray_CheckScalar(value)) {
@@ -247,23 +263,6 @@ QuadPrecision_from_object(PyObject *value, QuadBackendType backend)
     else if (PyLong_Check(value)) {
         return quad_from_py_int(value, backend, self);
     }
-    else if (Py_TYPE(value) == &QuadPrecision_Type) {
-        Py_DECREF(self);  // discard the default one
-        QuadPrecisionObject *quad_obj = (QuadPrecisionObject *)value;
-        // create a new one with the same backend
-        QuadPrecisionObject *self = QuadPrecision_raw_new(quad_obj->backend);
-        if (self == NULL) {
-            return NULL;
-        }
-        if (quad_obj->backend == BACKEND_SLEEF) {
-            self->value.sleef_value = quad_obj->value.sleef_value;
-        }
-        else {
-            self->value.longdouble_value = quad_obj->value.longdouble_value;
-        }
-
-        return self;
-    }
     else {
         PyObject *type_str = PyObject_Str((PyObject *)Py_TYPE(value));
         if (type_str != NULL) {

src/csrc/scalar_ops.cpp

@@ -19,6 +19,7 @@ extern "C" {
 #include "ops.hpp"
 #include "scalar_ops.h"
 #include "quad_common.h"
+#include "dtype.h"
 #include "constants.hpp"
 
 template <unary_op_quad_def sleef_op, unary_op_longdouble_def longdouble_op>
@@ -77,7 +78,7 @@ quad_binary_func(PyObject *op1, PyObject *op2)
         Py_INCREF(other);
         other_quad = (QuadPrecisionObject *)other;
         if (other_quad->backend != backend) {
-            PyErr_SetString(PyExc_TypeError, "Cannot mix QuadPrecision backends");
+            quad_set_backend_mismatch_error();
             Py_DECREF(other_quad);
             Py_DECREF(other);
             return NULL;
@@ -138,11 +139,7 @@ quad_richcompare(QuadPrecisionObject *self, PyObject *other, int cmp_op)
         Py_INCREF(other);
         other_quad = (QuadPrecisionObject *)other;
         if (other_quad->backend != backend) {
-            // we could allow, but this will be bad
-            // Two values that are different in quad precision, 
-            // might appear equal when converted to double.
-            PyErr_SetString(PyExc_TypeError,
-                            "Cannot compare QuadPrecision objects with different backends");
+            quad_set_backend_mismatch_error();
             Py_DECREF(other_quad);
             return NULL;
         }

src/csrc/umath/binary_ops.cpp

@@ -32,17 +32,9 @@ quad_binary_op_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtyp
 
     QuadPrecDTypeObject *descr_in1 = (QuadPrecDTypeObject *)given_descrs[0];
     QuadPrecDTypeObject *descr_in2 = (QuadPrecDTypeObject *)given_descrs[1];
-    QuadBackendType target_backend;
-
-    // Determine target backend and if casting is needed
-    NPY_CASTING casting = NPY_NO_CASTING;
-    if (descr_in1->backend != descr_in2->backend) {
-        target_backend = BACKEND_LONGDOUBLE;
-        casting = NPY_SAFE_CASTING;
-    }
-    else {
-        target_backend = descr_in1->backend;
-    }
+    QuadBackendType target_backend = quad_loop_backend(given_descrs, 2);
+    NPY_CASTING casting = (descr_in1->backend == descr_in2->backend) ? NPY_NO_CASTING
+                                                                     : NPY_SAFE_CASTING;
 
     // Set up input descriptors, casting if necessary
     for (int i = 0; i < 2; i++) {
@@ -161,17 +153,9 @@ quad_binary_op_2out_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const
 
     QuadPrecDTypeObject *descr_in1 = (QuadPrecDTypeObject *)given_descrs[0];
     QuadPrecDTypeObject *descr_in2 = (QuadPrecDTypeObject *)given_descrs[1];
-    QuadBackendType target_backend;
-
-    // Determine target backend and if casting is needed
-    NPY_CASTING casting = NPY_NO_CASTING;
-    if (descr_in1->backend != descr_in2->backend) {
-        target_backend = BACKEND_LONGDOUBLE;
-        casting = NPY_SAFE_CASTING;
-    }
-    else {
-        target_backend = descr_in1->backend;
-    }
+    QuadBackendType target_backend = quad_loop_backend(given_descrs, 2);
+    NPY_CASTING casting = (descr_in1->backend == descr_in2->backend) ? NPY_NO_CASTING
+                                                                     : NPY_SAFE_CASTING;
 
     // Set up input descriptors, casting if necessary
     for (int i = 0; i < 2; i++) {

src/csrc/umath/comparison_ops.cpp

@@ -36,17 +36,9 @@ quad_comparison_op_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const
 
     QuadPrecDTypeObject *descr_in1 = (QuadPrecDTypeObject *)given_descrs[0];
     QuadPrecDTypeObject *descr_in2 = (QuadPrecDTypeObject *)given_descrs[1];
-    QuadBackendType target_backend;
-
-    // As dealing with different backends then cast to boolean
-    NPY_CASTING casting = NPY_NO_CASTING;
-    if (descr_in1->backend != descr_in2->backend) {
-        target_backend = BACKEND_LONGDOUBLE;
-        casting = NPY_SAFE_CASTING;
-    }
-    else {
-        target_backend = descr_in1->backend;
-    }
+    QuadBackendType target_backend = quad_loop_backend(given_descrs, 2);
+    NPY_CASTING casting = (descr_in1->backend == descr_in2->backend) ? NPY_NO_CASTING
+                                                                     : NPY_SAFE_CASTING;
 
     // Set up input descriptors, casting if necessary
     for (int i = 0; i < 2; i++) {

src/csrc/umath/unary_ops.cpp

@@ -29,22 +29,7 @@ quad_unary_op_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtype
     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];
-    }
-
-    QuadPrecDTypeObject *descr_in = (QuadPrecDTypeObject *)given_descrs[0];
-    QuadPrecDTypeObject *descr_out = (QuadPrecDTypeObject *)loop_descrs[1];
-
-    if (descr_in->backend != descr_out->backend) {
-        return NPY_UNSAFE_CASTING;
-    }
-
+    loop_descrs[1] = quad_out_loop_descr(given_descrs[0], given_descrs[1]);
     return NPY_NO_CASTING;
 }
 
@@ -288,26 +273,9 @@ quad_unary_op_2out_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const
     Py_INCREF(given_descrs[0]);
     loop_descrs[0] = given_descrs[0];
 
-    // Output descriptors (2 outputs)
     for (int i = 1; i < 3; i++) {
-        if (given_descrs[i] == NULL) {
-            Py_INCREF(given_descrs[0]);
-            loop_descrs[i] = given_descrs[0];
-        }
-        else {
-            Py_INCREF(given_descrs[i]);
-            loop_descrs[i] = given_descrs[i];
-        }
+        loop_descrs[i] = quad_out_loop_descr(given_descrs[0], given_descrs[i]);
     }
-
-    QuadPrecDTypeObject *descr_in = (QuadPrecDTypeObject *)given_descrs[0];
-    QuadPrecDTypeObject *descr_out1 = (QuadPrecDTypeObject *)loop_descrs[1];
-    QuadPrecDTypeObject *descr_out2 = (QuadPrecDTypeObject *)loop_descrs[2];
-
-    if (descr_in->backend != descr_out1->backend || descr_in->backend != descr_out2->backend) {
-        return NPY_UNSAFE_CASTING;
-    }
-
     return NPY_NO_CASTING;
 }
 
@@ -442,14 +410,7 @@ quad_frexp_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtypes[]
     loop_descrs[0] = given_descrs[0];
 
     // Output 1: mantissa (same dtype as input)
-    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];
-    }
+    loop_descrs[1] = quad_out_loop_descr(given_descrs[0], given_descrs[1]);
 
     // Output 2: exponent (int32)
     if (given_descrs[2] == NULL) {
@@ -462,14 +423,6 @@ quad_frexp_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtypes[]
         Py_INCREF(given_descrs[2]);
         loop_descrs[2] = given_descrs[2];
     }
-
-    QuadPrecDTypeObject *descr_in = (QuadPrecDTypeObject *)given_descrs[0];
-    QuadPrecDTypeObject *descr_out1 = (QuadPrecDTypeObject *)loop_descrs[1];
-
-    if (descr_in->backend != descr_out1->backend) {
-        return NPY_UNSAFE_CASTING;
-    }
-
     return NPY_NO_CASTING;
 }