FEAT: modf ufunc implementation

quaddtype/numpy_quaddtype/src/ops.hpp

@@ -10,6 +10,8 @@
 
 // Unary Quad Operations
 typedef Sleef_quad (*unary_op_quad_def)(const Sleef_quad *);
+// Unary Quad operations with 2 outputs (for modf, frexp)
+typedef void (*unary_op_2out_quad_def)(const Sleef_quad *, Sleef_quad *, Sleef_quad *);
 
 static Sleef_quad
 quad_negative(const Sleef_quad *op)
@@ -1361,6 +1363,23 @@ ld_spacing(const long double *x)
     return result;
 }
 
+// Unary operations with 2 outputs
+static inline void
+quad_modf(const Sleef_quad *a, Sleef_quad *out_fractional, Sleef_quad *out_integral)
+{
+    // int part stored in out_integral
+    *out_fractional = Sleef_modfq1(*a, out_integral);
+}
+
+// Unary long double operations with 2 outputs  
+typedef void (*unary_op_2out_longdouble_def)(const long double *, long double *, long double *);
+
+static inline void
+ld_modf(const long double *a, long double *out_fractional, long double *out_integral)
+{
+    *out_fractional = modfl(*a, out_integral);
+}
+
 // comparison quad functions
 typedef npy_bool (*cmp_quad_def)(const Sleef_quad *, const Sleef_quad *);
 

quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp

@@ -267,6 +267,147 @@ create_quad_logical_not_ufunc(PyObject *numpy, const char *ufunc_name)
     return 0;
 }
 
+// Resolver for unary ufuncs with 2 outputs (like modf)
+static NPY_CASTING
+quad_unary_op_2out_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtypes[],
+                                       PyArray_Descr *const given_descrs[], PyArray_Descr *loop_descrs[],
+                                       npy_intp *NPY_UNUSED(view_offset))
+{
+    // Input descriptor
+    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];
+        }
+    }
+
+    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;
+}
+
+// Strided loop for unary ops with 2 outputs (unaligned)
+template <unary_op_2out_quad_def sleef_op, unary_op_2out_longdouble_def longdouble_op>
+int
+quad_generic_unary_op_2out_strided_loop_unaligned(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 *out1_ptr = data[1];
+    char *out2_ptr = data[2];
+    npy_intp in_stride = strides[0];
+    npy_intp out1_stride = strides[1];
+    npy_intp out2_stride = strides[2];
+
+    QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0];
+    QuadBackendType backend = descr->backend;
+    size_t elem_size = (backend == BACKEND_SLEEF) ? sizeof(Sleef_quad) : sizeof(long double);
+
+    quad_value in, out1, out2;
+    while (N--) {
+        memcpy(&in, in_ptr, elem_size);
+        if (backend == BACKEND_SLEEF) {
+            sleef_op(&in.sleef_value, &out1.sleef_value, &out2.sleef_value);
+        }
+        else {
+            longdouble_op(&in.longdouble_value, &out1.longdouble_value, &out2.longdouble_value);
+        }
+        memcpy(out1_ptr, &out1, elem_size);
+        memcpy(out2_ptr, &out2, elem_size);
+
+        in_ptr += in_stride;
+        out1_ptr += out1_stride;
+        out2_ptr += out2_stride;
+    }
+    return 0;
+}
+
+// Strided loop for unary ops with 2 outputs (aligned)
+template <unary_op_2out_quad_def sleef_op, unary_op_2out_longdouble_def longdouble_op>
+int
+quad_generic_unary_op_2out_strided_loop_aligned(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 *out1_ptr = data[1];
+    char *out2_ptr = data[2];
+    npy_intp in_stride = strides[0];
+    npy_intp out1_stride = strides[1];
+    npy_intp out2_stride = strides[2];
+
+    QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0];
+    QuadBackendType backend = descr->backend;
+
+    while (N--) {
+        if (backend == BACKEND_SLEEF) {
+            sleef_op((Sleef_quad *)in_ptr, (Sleef_quad *)out1_ptr, (Sleef_quad *)out2_ptr);
+        }
+        else {
+            longdouble_op((long double *)in_ptr, (long double *)out1_ptr, (long double *)out2_ptr);
+        }
+        in_ptr += in_stride;
+        out1_ptr += out1_stride;
+        out2_ptr += out2_stride;
+    }
+    return 0;
+}
+
+// Create unary ufunc with 2 outputs
+template <unary_op_2out_quad_def sleef_op, unary_op_2out_longdouble_def longdouble_op>
+int
+create_quad_unary_2out_ufunc(PyObject *numpy, const char *ufunc_name)
+{
+    PyObject *ufunc = PyObject_GetAttrString(numpy, ufunc_name);
+    if (ufunc == NULL) {
+        return -1;
+    }
+
+    // 1 input, 2 outputs
+    PyArray_DTypeMeta *dtypes[3] = {&QuadPrecDType, &QuadPrecDType, &QuadPrecDType};
+
+    PyType_Slot slots[] = {
+            {NPY_METH_resolve_descriptors, (void *)&quad_unary_op_2out_resolve_descriptors},
+            {NPY_METH_strided_loop,
+             (void *)&quad_generic_unary_op_2out_strided_loop_aligned<sleef_op, longdouble_op>},
+            {NPY_METH_unaligned_strided_loop,
+             (void *)&quad_generic_unary_op_2out_strided_loop_unaligned<sleef_op, longdouble_op>},
+            {0, NULL}};
+
+    PyArrayMethod_Spec Spec = {
+            .name = "quad_unary_2out",
+            .nin = 1,
+            .nout = 2,
+            .casting = NPY_NO_CASTING,
+            .flags = NPY_METH_SUPPORTS_UNALIGNED,
+            .dtypes = dtypes,
+            .slots = slots,
+    };
+
+    if (PyUFunc_AddLoopFromSpec(ufunc, &Spec) < 0) {
+        return -1;
+    }
+
+    return 0;
+}
+
 int
 init_quad_unary_ops(PyObject *numpy)
 {
@@ -394,5 +535,10 @@ init_quad_unary_ops(PyObject *numpy)
         return -1;
     }
 
+    // Unary operations with 2 outputs
+    if (create_quad_unary_2out_ufunc<quad_modf, ld_modf>(numpy, "modf") < 0) {
+        return -1;
+    }
+
     return 0;
 }

quaddtype/release_tracker.md

@@ -79,7 +79,7 @@
 | copysign      | ✅    | ✅                                                                   |
 | nextafter     | ✅    | ✅                                                                   |
 | spacing       | ✅    | ✅                                                                   |
-| modf          |       |                                                                      |
+| modf          | ✅    | ✅                                                                   |
 | ldexp         |       |                                                                      |
 | frexp         |       |                                                                      |
 | floor         | ✅    | ✅                                                                   |

quaddtype/tests/test_quaddtype.py

@@ -2429,7 +2429,7 @@ def test_subnormal_range(self, x):
         if result != q_zero:
             assert np.signbit(result) == np.signbit(q_x), \
                 f"spacing({x}) should have same sign as {x}"
-
+            
     def test_smallest_normal_spacing(self):
         """Test spacing for smallest normal value and 2*smallest normal"""
         finfo = np.finfo(QuadPrecDType())
@@ -2457,3 +2457,76 @@ def test_smallest_normal_spacing(self):
             f"spacing(2*smallest_normal) should be {expected2}, got {result2}"
         assert result2 > q_zero, "Result should be positive"
 
+
+class TestModf:
+    """Test cases for np.modf function with QuadPrecision dtype"""
+
+    @pytest.mark.parametrize("x", [
+        # Basic positive/negative numbers
+        "3.14159", "-3.14159", "2.71828", "-2.71828", "1.5", "-1.5", "0.75", "-0.75",
+        # Integers (fractional part should be zero)
+        "0.0", "-0.0", "1.0", "-1.0", "5.0", "-5.0", "42.0", "-42.0",
+        # Small numbers
+        "0.001", "-0.001", "0.000123", "-0.000123",
+        # Large numbers  
+        "1e10", "-1e10", "1e15", "-1e15",
+        # Numbers close to integers
+        "0.999999999999", "-0.999999999999", "1.000000000001", "-1.000000000001",
+        # Special values
+        "inf", "-inf", "nan",
+        # Edge cases for sign consistency
+        "5.7", "-5.7", "0.3", "-0.3"
+    ])
+    def test_modf(self, x):
+        """Test modf against NumPy's behavior"""
+        quad_x = QuadPrecision(x)
+
+        # Compute modf for both QuadPrecision and float64
+        quad_frac, quad_int = np.modf(quad_x)
+
+        # Create numpy float64 for reference
+        try:
+            float_x = np.float64(x)
+            np_frac, np_int = np.modf(float_x)
+        except (ValueError, OverflowError):
+            # Handle cases where string can't convert to float64 (like "nan")
+            float_x = np.float64(float(x))
+            np_frac, np_int = np.modf(float_x)
+
+        # Check return types
+        assert isinstance(quad_frac, QuadPrecision), f"Fractional part should be QuadPrecision for {x}"
+        assert isinstance(quad_int, QuadPrecision), f"Integral part should be QuadPrecision for {x}"
+
+        # Direct comparison with NumPy results
+        if np.isnan(np_frac):
+            assert np.isnan(quad_frac), f"Expected NaN fractional part for modf({x})"
+        else:
+            np.testing.assert_allclose(
+                quad_frac, np_frac, rtol=1e-12, atol=1e-15,
+                err_msg=f"Fractional part mismatch for modf({x})"
+            )
+
+        if np.isnan(np_int):
+            assert np.isnan(quad_int), f"Expected NaN integral part for modf({x})"  
+        elif np.isinf(np_int):
+            assert np.isinf(quad_int), f"Expected inf integral part for modf({x})"
+            assert np.signbit(quad_int) == np.signbit(np_int), f"Sign mismatch for inf integral part modf({x})"
+        else:
+            np.testing.assert_allclose(
+                quad_int, np_int, rtol=1e-12, atol=0,
+                err_msg=f"Integral part mismatch for modf({x})"
+            )
+
+        # Check sign preservation for zero values
+        if np_frac == 0.0:
+            assert np.signbit(quad_frac) == np.signbit(np_frac), f"Zero fractional sign mismatch for modf({x})"
+        if np_int == 0.0:
+            assert np.signbit(quad_int) == np.signbit(np_int), f"Zero integral sign mismatch for modf({x})"
+
+        # Verify reconstruction property for finite values
+        if np.isfinite(float_x) and not np.isnan(float_x):
+            reconstructed = quad_int + quad_frac
+            np.testing.assert_allclose(
+                reconstructed, quad_x, rtol=1e-12, atol=1e-15,
+                err_msg=f"Reconstruction failed for modf({x}): {quad_int} + {quad_frac} != {quad_x}"
+            )