Implement cast support for half

Author: juntyr

quaddtype/meson.build

@@ -23,9 +23,13 @@ incdir_numpy = run_command(py,
   check : true
 ).stdout().strip()
 
+npymath_path = incdir_numpy / '..' / 'lib'
+npymath_lib = cc.find_library('npymath', dirs: npymath_path)
+
+dependencies = [sleef_dep, py_dep, npymath_lib]
+
 # Add OpenMP dependency (optional, for threading)
 openmp_dep = dependency('openmp', required: false)
-dependencies = [sleef_dep, py_dep]
 if openmp_dep.found()
     dependencies += openmp_dep
 endif

quaddtype/numpy_quaddtype/src/casts.cpp

@@ -9,6 +9,7 @@ extern "C" {
 #include <Python.h>
 
 #include "numpy/arrayobject.h"
+#include "numpy/halffloat.h"
 #include "numpy/ndarraytypes.h"
 #include "numpy/dtype_api.h"
 }
@@ -20,7 +21,7 @@ extern "C" {
 #include "casts.h"
 #include "dtype.h"
 
-#define NUM_CASTS 31  // 15 to_casts + 15 from_casts + 1 quad_to_quad
+#define NUM_CASTS 33  // 16 to_casts + 16 from_casts + 1 quad_to_quad
 
 static NPY_CASTING
 quad_to_quad_resolve_descriptors(PyObject *NPY_UNUSED(self),
@@ -150,19 +151,21 @@ quad_to_quad_strided_loop_aligned(PyArrayMethod_Context *context, char *const da
     return 0;
 }
 
-// Template magic to ensure npy_bool and npy_ubyte do not alias in templates
+// Template magic to ensure npy_bool/npy_ubyte and npy_half/npy_ushort do not alias in templates
 struct my_npy_bool {}
+struct my_npy_half {}
 
 template<typename T>
 struct NpyType { typedef T TYPE; };
 template<>
 struct NpyType<my_npy_bool>{ typedef npy_bool TYPE; };
+struct NpyType<my_npy_half>{ typedef npy_half TYPE; };
 
 // Casting from other types to QuadDType
 
 template <typename T>
 static inline quad_value
-to_quad(NpyType<T>::TYPE x, QuadBackendType backend);
+to_quad(typename NpyType<T>::TYPE x, QuadBackendType backend);
 
 template <>
 inline quad_value
@@ -317,6 +320,21 @@ to_quad<npy_ulonglong>(npy_ulonglong x, QuadBackendType backend)
     }
     return result;
 }
+
+template <>
+inline quad_value
+to_quad<my_npy_half>(npy_half x, QuadBackendType backend)
+{
+    quad_value result;
+    if (backend == BACKEND_SLEEF) {
+        result.sleef_value = Sleef_cast_from_doubleq1(npy_half_to_double(x));
+    }
+    else {
+        result.longdouble_value = (long double)npy_half_to_double(x);
+    }
+    return result;
+}
+
 template <>
 inline quad_value
 to_quad<float>(float x, QuadBackendType backend)
@@ -396,10 +414,10 @@ numpy_to_quad_strided_loop_unaligned(PyArrayMethod_Context *context, char *const
     size_t elem_size = (backend == BACKEND_SLEEF) ? sizeof(Sleef_quad) : sizeof(long double);
 
     while (N--) {
-        NpyType<T>::TYPE in_val;
+        typename NpyType<T>::TYPE in_val;
         quad_value out_val;
 
-        memcpy(&in_val, in_ptr, sizeof(NpyType<T>::TYPE));
+        memcpy(&in_val, in_ptr, sizeof(typename NpyType<T>::TYPE));
         out_val = to_quad<T>(in_val, backend);
         memcpy(out_ptr, &out_val, elem_size);
 
@@ -423,7 +441,7 @@ numpy_to_quad_strided_loop_aligned(PyArrayMethod_Context *context, char *const d
     QuadBackendType backend = descr_out->backend;
 
     while (N--) {
-        NpyType<T>::TYPE in_val = *(NpyType<T>::TYPE *)in_ptr;
+        typename NpyType<T>::TYPE in_val = *(typename NpyType<T>::TYPE *)in_ptr;
         quad_value out_val = to_quad<T>(in_val, backend);
 
         if (backend == BACKEND_SLEEF) {
@@ -442,7 +460,7 @@ numpy_to_quad_strided_loop_aligned(PyArrayMethod_Context *context, char *const d
 // Casting from QuadDType to other types
 
 template <typename T>
-static inline NpyType<T>::TYPE
+static inline typename NpyType<T>::TYPE
 from_quad(quad_value x, QuadBackendType backend);
 
 template <>
@@ -577,6 +595,18 @@ from_quad<npy_ulonglong>(quad_value x, QuadBackendType backend)
     }
 }
 
+template <>
+inline npy_half
+from_quad<my_npy_half>(quad_value x, QuadBackendType backend)
+{
+    if (backend == BACKEND_SLEEF) {
+        return npy_double_to_half(Sleef_cast_to_doubleq1(x.sleef_value));
+    }
+    else {
+        return npy_double_to_half((double)x.longdouble_value);
+    }
+}
+
 template <>
 inline float
 from_quad<float>(quad_value x, QuadBackendType backend)
@@ -645,8 +675,8 @@ quad_to_numpy_strided_loop_unaligned(PyArrayMethod_Context *context, char *const
         quad_value in_val;
         memcpy(&in_val, in_ptr, elem_size);
 
-        NpyType<T>::TYPE out_val = from_quad<T>(in_val, backend);
-        memcpy(out_ptr, &out_val, sizeof(NpyType<T>::TYPE));
+        typename NpyType<T>::TYPE out_val = from_quad<T>(in_val, backend);
+        memcpy(out_ptr, &out_val, sizeof(typename NpyType<T>::TYPE));
 
         in_ptr += strides[0];
         out_ptr += strides[1];
@@ -676,8 +706,8 @@ quad_to_numpy_strided_loop_aligned(PyArrayMethod_Context *context, char *const d
             in_val.longdouble_value = *(long double *)in_ptr;
         }
 
-        NpyType<T>::TYPE out_val = from_quad<T>(in_val, backend);
-        *(NpyType<T>::TYPE *)(out_ptr) = out_val;
+        typename NpyType<T>::TYPE out_val = from_quad<T>(in_val, backend);
+        *(typename NpyType<T>::TYPE *)(out_ptr) = out_val;
 
         in_ptr += strides[0];
         out_ptr += strides[1];
@@ -784,6 +814,7 @@ init_casts_internal(void)
     add_cast_to<npy_ulong>(&PyArray_ULongDType);
     add_cast_to<npy_longlong>(&PyArray_LongLongDType);
     add_cast_to<npy_ulonglong>(&PyArray_ULongLongDType);
+    add_cast_to<my_npy_half>(&PyArray_HalfDType);
     add_cast_to<float>(&PyArray_FloatDType);
     add_cast_to<double>(&PyArray_DoubleDType);
     add_cast_to<long double>(&PyArray_LongDoubleDType);
@@ -799,6 +830,7 @@ init_casts_internal(void)
     add_cast_from<npy_ulong>(&PyArray_ULongDType);
     add_cast_from<npy_longlong>(&PyArray_LongLongDType);
     add_cast_from<npy_ulonglong>(&PyArray_ULongLongDType);
+    add_cast_from<my_npy_half>(&PyArray_HalfDType);
     add_cast_from<float>(&PyArray_FloatDType);
     add_cast_from<double>(&PyArray_DoubleDType);
     add_cast_from<long double>(&PyArray_LongDoubleDType);

quaddtype/tests/test_quaddtype.py

@@ -41,9 +41,6 @@ def test_finfo_int_constant(name, value):
 
 @pytest.mark.parametrize("dtype", ["bool", "byte", "int8", "ubyte", "uint8", "short", "int16", "ushort", "uint16", "int", "int32", "uint", "uint32", "long", "ulong", "longlong", "int64", "ulonglong", "uint64", "half", "float16", "float", "float32", "double", "float64", "longdouble"])
 def test_astype(dtype):
-    if dtype in ("half", "float16"):
-        pytest.xfail(f"{dtype} astype not yet supported")
-
     orig = np.array(1, dtype=dtype)
     quad = orig.astype(QuadPrecDType, casting="safe")
     back = quad.astype(dtype, casting="unsafe")