diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index 40ee029b..638a5706 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/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 *); diff --git a/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp index c2710747..32cae501 100644 --- a/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp +++ b/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 +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 +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 +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}, + {NPY_METH_unaligned_strided_loop, + (void *)&quad_generic_unary_op_2out_strided_loop_unaligned}, + {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(numpy, "modf") < 0) { + return -1; + } + return 0; } \ No newline at end of file diff --git a/quaddtype/release_tracker.md b/quaddtype/release_tracker.md index c5e1ad2b..e5be6e88 100644 --- a/quaddtype/release_tracker.md +++ b/quaddtype/release_tracker.md @@ -79,7 +79,7 @@ | copysign | ✅ | ✅ | | nextafter | ✅ | ✅ | | spacing | ✅ | ✅ | -| modf | | | +| modf | ✅ | ✅ | | ldexp | | | | frexp | | | | floor | ✅ | ✅ | diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index 693b9090..5662a643 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/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}" + ) \ No newline at end of file