From 76bc5e4b2673f59fe2b41bbe8d106998a98bd943 Mon Sep 17 00:00:00 2001 From: swayaminsync Date: Fri, 17 Oct 2025 12:59:53 +0530 Subject: [PATCH] divmod impl --- quaddtype/numpy_quaddtype/src/ops.hpp | 20 ++ .../numpy_quaddtype/src/umath/binary_ops.cpp | 198 +++++++++++++++++ quaddtype/release_tracker.md | 4 +- quaddtype/tests/test_quaddtype.py | 203 ++++++++++++++++++ 4 files changed, 423 insertions(+), 2 deletions(-) diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index 08973f8a..dd385d92 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/quaddtype/numpy_quaddtype/src/ops.hpp @@ -442,6 +442,8 @@ ld_isnan(const long double *op) // Binary Quad operations typedef Sleef_quad (*binary_op_quad_def)(const Sleef_quad *, const Sleef_quad *); +// Binary Quad operations with 2 outputs (for divmod, modf, frexp) +typedef void (*binary_op_2out_quad_def)(const Sleef_quad *, const Sleef_quad *, Sleef_quad *, Sleef_quad *); static inline Sleef_quad quad_add(const Sleef_quad *in1, const Sleef_quad *in2) @@ -593,6 +595,14 @@ quad_fmod(const Sleef_quad *a, const Sleef_quad *b) return result; } +static inline void +quad_divmod(const Sleef_quad *a, const Sleef_quad *b, + Sleef_quad *out_quotient, Sleef_quad *out_remainder) +{ + *out_quotient = quad_floor_divide(a, b); + *out_remainder = quad_mod(a, b); +} + static inline Sleef_quad quad_minimum(const Sleef_quad *in1, const Sleef_quad *in2) { @@ -744,6 +754,8 @@ quad_logaddexp2(const Sleef_quad *x, const Sleef_quad *y) // Binary long double operations typedef long double (*binary_op_longdouble_def)(const long double *, const long double *); +// Binary long double operations with 2 outputs (for divmod, modf, frexp) +typedef void (*binary_op_2out_longdouble_def)(const long double *, const long double *, long double *, long double *); static inline long double ld_add(const long double *in1, const long double *in2) @@ -871,6 +883,14 @@ ld_fmod(const long double *a, const long double *b) return result; } +static inline void +ld_divmod(const long double *a, const long double *b, + long double *out_quotient, long double *out_remainder) +{ + *out_quotient = ld_floor_divide(a, b); + *out_remainder = ld_mod(a, b); +} + static inline long double ld_minimum(const long double *in1, const long double *in2) { diff --git a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp index 4c9634ae..7775b706 100644 --- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp @@ -145,6 +145,201 @@ quad_generic_binop_strided_loop_aligned(PyArrayMethod_Context *context, char *co } + +// Resolve descriptors for binary ops with 2 outputs (2 inputs, 2 outputs) +static NPY_CASTING +quad_binary_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)) +{ + 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; + } + + // Set up input descriptors, casting if necessary + for (int i = 0; i < 2; i++) { + if (((QuadPrecDTypeObject *)given_descrs[i])->backend != target_backend) { + loop_descrs[i] = (PyArray_Descr *)new_quaddtype_instance(target_backend); + if (!loop_descrs[i]) { + return (NPY_CASTING)-1; + } + } + else { + Py_INCREF(given_descrs[i]); + loop_descrs[i] = given_descrs[i]; + } + } + + // Set up output descriptors (2 outputs for divmod) + for (int i = 2; i < 4; i++) { + if (given_descrs[i] == NULL) { + loop_descrs[i] = (PyArray_Descr *)new_quaddtype_instance(target_backend); + if (!loop_descrs[i]) { + return (NPY_CASTING)-1; + } + } + else { + QuadPrecDTypeObject *descr_out = (QuadPrecDTypeObject *)given_descrs[i]; + if (descr_out->backend != target_backend) { + loop_descrs[i] = (PyArray_Descr *)new_quaddtype_instance(target_backend); + if (!loop_descrs[i]) { + return (NPY_CASTING)-1; + } + } + else { + Py_INCREF(given_descrs[i]); + loop_descrs[i] = given_descrs[i]; + } + } + } + return casting; +} + +// Strided loop for binary ops with 2 outputs (unaligned) +template +int +quad_generic_binop_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 *in1_ptr = data[0], *in2_ptr = data[1]; + char *out1_ptr = data[2], *out2_ptr = data[3]; + npy_intp in1_stride = strides[0]; + npy_intp in2_stride = strides[1]; + npy_intp out1_stride = strides[2]; + npy_intp out2_stride = strides[3]; + + 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 in1, in2, out1, out2; + while (N--) { + memcpy(&in1, in1_ptr, elem_size); + memcpy(&in2, in2_ptr, elem_size); + if (backend == BACKEND_SLEEF) { + sleef_op(&in1.sleef_value, &in2.sleef_value, &out1.sleef_value, &out2.sleef_value); + } + else { + longdouble_op(&in1.longdouble_value, &in2.longdouble_value, + &out1.longdouble_value, &out2.longdouble_value); + } + memcpy(out1_ptr, &out1, elem_size); + memcpy(out2_ptr, &out2, elem_size); + + in1_ptr += in1_stride; + in2_ptr += in2_stride; + out1_ptr += out1_stride; + out2_ptr += out2_stride; + } + return 0; +} + +// Strided loop for binary ops with 2 outputs (aligned) +template +int +quad_generic_binop_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 *in1_ptr = data[0], *in2_ptr = data[1]; + char *out1_ptr = data[2], *out2_ptr = data[3]; + npy_intp in1_stride = strides[0]; + npy_intp in2_stride = strides[1]; + npy_intp out1_stride = strides[2]; + npy_intp out2_stride = strides[3]; + + QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0]; + QuadBackendType backend = descr->backend; + + while (N--) { + if (backend == BACKEND_SLEEF) { + sleef_op((Sleef_quad *)in1_ptr, (Sleef_quad *)in2_ptr, + (Sleef_quad *)out1_ptr, (Sleef_quad *)out2_ptr); + } + else { + longdouble_op((long double *)in1_ptr, (long double *)in2_ptr, + (long double *)out1_ptr, (long double *)out2_ptr); + } + + in1_ptr += in1_stride; + in2_ptr += in2_stride; + out1_ptr += out1_stride; + out2_ptr += out2_stride; + } + return 0; +} + +// Create binary ufunc with 2 outputs (generic for divmod, modf, frexp, etc.) +template +int +create_quad_binary_2out_ufunc(PyObject *numpy, const char *ufunc_name) +{ + PyObject *ufunc = PyObject_GetAttrString(numpy, ufunc_name); + if (ufunc == NULL) { + return -1; + } + + // 2 inputs, 2 outputs + PyArray_DTypeMeta *dtypes[4] = {&QuadPrecDType, &QuadPrecDType, &QuadPrecDType, &QuadPrecDType}; + + PyType_Slot slots[] = { + {NPY_METH_resolve_descriptors, (void *)&quad_binary_op_2out_resolve_descriptors}, + {NPY_METH_strided_loop, + (void *)&quad_generic_binop_2out_strided_loop_aligned}, + {NPY_METH_unaligned_strided_loop, + (void *)&quad_generic_binop_2out_strided_loop_unaligned}, + {0, NULL}}; + + PyArrayMethod_Spec Spec = { + .name = "quad_binop_2out", + .nin = 2, + .nout = 2, + .casting = NPY_NO_CASTING, + .flags = (NPY_ARRAYMETHOD_FLAGS)(NPY_METH_SUPPORTS_UNALIGNED | NPY_METH_IS_REORDERABLE), + .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(4, &PyArrayDescr_Type, &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; +} + template int create_quad_binary_ufunc(PyObject *numpy, const char *ufunc_name) @@ -259,5 +454,8 @@ init_quad_binary_ops(PyObject *numpy) if (create_quad_binary_ufunc(numpy, "logaddexp2") < 0) { return -1; } + if (create_quad_binary_2out_ufunc(numpy, "divmod") < 0) { + return -1; + } return 0; } \ No newline at end of file diff --git a/quaddtype/release_tracker.md b/quaddtype/release_tracker.md index 74c75315..8a812bfd 100644 --- a/quaddtype/release_tracker.md +++ b/quaddtype/release_tracker.md @@ -1,6 +1,6 @@ # Plan for `numpy-quaddtype` v1.0.0 -- [ ] High-Endian System support +- [x] High-Endian System support - [ ] Complete Documentation | ufunc name | Added | Edge Cases Tested\* | @@ -21,7 +21,7 @@ | remainder | ✅ | ✅ | | mod | ✅ | ✅ | | fmod | ✅ | ✅ | -| divmod | | | +| divmod | ✅ | ✅ | | absolute | ✅ | ✅ | | fabs | | | | rint | ✅ | ✅ | diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index 3c765b5b..9af9ea78 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -1077,6 +1077,209 @@ def test_mod(a, b, backend, op): assert result_negative == numpy_negative, f"Sign mismatch for {a} % {b}: quad={result_negative}, numpy={numpy_negative}" +@pytest.mark.parametrize("backend", ["sleef", "longdouble"]) +@pytest.mark.parametrize("a,b", [ + # Basic cases - positive/positive + (7.0, 3.0), (10.5, 3.2), (21.0, 4.0), + + # Positive/negative combinations + (-7.0, 3.0), (7.0, -3.0), (-7.0, -3.0), + (-10.5, 3.2), (10.5, -3.2), (-10.5, -3.2), + + # Zero dividend cases + (0.0, 3.0), (-0.0, 3.0), (0.0, -3.0), (-0.0, -3.0), + + # Cases that result in zero remainder (exact division) + (6.0, 3.0), (-6.0, 3.0), (6.0, -3.0), (-6.0, -3.0), + (1.0, 1.0), (-1.0, 1.0), (1.0, -1.0), (-1.0, -1.0), + (10.0, 2.0), (-10.0, 2.0), (10.0, -2.0), (-10.0, -2.0), + + # Fractional cases + (7.5, 2.5), (-7.5, 2.5), (7.5, -2.5), (-7.5, -2.5), + (0.75, 0.25), (-0.1, 0.3), (0.9, -1.0), (-1.1, -1.0), + (3.14159, 1.0), (-3.14159, 1.0), (3.14159, -1.0), (-3.14159, -1.0), + + # Large/small numbers + (1e10, 1e5), (-1e10, 1e5), (1e-10, 1e-5), (-1e-10, 1e-5), + (1e15, 1e10), (1e-15, 1e-10), + + # Finite % infinity cases + (5.0, float('inf')), (-5.0, float('inf')), + (5.0, float('-inf')), (-5.0, float('-inf')), + (0.0, float('inf')), (-0.0, float('-inf')), + + # NaN cases (should return NaN for both quotient and remainder) + (float('nan'), 3.0), (3.0, float('nan')), (float('nan'), float('nan')), + + # Division by zero cases (should return inf/NaN) + (5.0, 0.0), (-5.0, 0.0), (0.0, 0.0), (-0.0, 0.0), + + # Infinity dividend cases (should return NaN for both) + (float('inf'), 3.0), (float('-inf'), 3.0), + (float('inf'), float('inf')), (float('-inf'), float('-inf')), + + # Cases with dividend < divisor + (1.0, 10.0), (-1.0, 10.0), (1.0, -10.0), (-1.0, -10.0), + (0.5, 1.0), (0.1, 1.0), (0.001, 0.01), +]) +def test_divmod(a, b, backend): + """Comprehensive test for divmod operation against NumPy behavior""" + if backend == "sleef": + quad_a = QuadPrecision(str(a)) + quad_b = QuadPrecision(str(b)) + elif backend == "longdouble": + quad_a = QuadPrecision(a, backend='longdouble') + quad_b = QuadPrecision(b, backend='longdouble') + + float_a = np.float64(a) + float_b = np.float64(b) + + # Compute divmod + quad_quotient, quad_remainder = np.divmod(quad_a, quad_b) + numpy_quotient, numpy_remainder = np.divmod(float_a, float_b) + + # Verify quotient + if np.isnan(numpy_quotient): + assert np.isnan(float(quad_quotient)), \ + f"Expected NaN quotient for divmod({a}, {b})" + elif np.isinf(numpy_quotient): + assert np.isinf(float(quad_quotient)) and \ + np.sign(numpy_quotient) == np.sign(float(quad_quotient)), \ + f"Expected inf quotient with matching sign for divmod({a}, {b})" + else: + # Adaptive tolerance for large quotients due to float64 conversion precision loss + atol_q = abs(numpy_quotient) * 1e-8 if abs(numpy_quotient) > 1e6 else 1e-15 + np.testing.assert_allclose( + float(quad_quotient), numpy_quotient, rtol=1e-9, atol=atol_q, + err_msg=f"Quotient mismatch for divmod({a}, {b})" + ) + if numpy_quotient == 0.0: + assert np.signbit(numpy_quotient) == np.signbit(quad_quotient), \ + f"Zero quotient sign mismatch for divmod({a}, {b})" + + # Verify remainder + if np.isnan(numpy_remainder): + assert np.isnan(float(quad_remainder)), \ + f"Expected NaN remainder for divmod({a}, {b})" + elif np.isinf(numpy_remainder): + assert np.isinf(float(quad_remainder)) and \ + np.sign(numpy_remainder) == np.sign(float(quad_remainder)), \ + f"Expected inf remainder with matching sign for divmod({a}, {b})" + else: + # Standard tolerance for remainder comparison + np.testing.assert_allclose( + float(quad_remainder), numpy_remainder, rtol=1e-9, atol=1e-15, + err_msg=f"Remainder mismatch for divmod({a}, {b})" + ) + if numpy_remainder == 0.0: + assert np.signbit(numpy_remainder) == np.signbit(quad_remainder), \ + f"Zero remainder sign mismatch for divmod({a}, {b})" + elif not np.isnan(b) and not np.isinf(b) and b != 0.0: + assert (float(quad_remainder) < 0) == (numpy_remainder < 0), \ + f"Remainder sign mismatch for divmod({a}, {b})" + + # Verify the fundamental property: a = quotient * b + remainder (for finite values) + if not np.isnan(numpy_quotient) and not np.isinf(numpy_quotient) and \ + not np.isnan(numpy_remainder) and not np.isinf(numpy_remainder) and \ + not np.isnan(b) and not np.isinf(b) and b != 0.0: + reconstructed = float(quad_quotient) * float(quad_b) + float(quad_remainder) + np.testing.assert_allclose( + reconstructed, float(quad_a), rtol=1e-10, atol=1e-15, + err_msg=f"Property a = q*b + r failed for divmod({a}, {b})" + ) + + +def test_divmod_special_properties(): + """Test special mathematical properties of divmod""" + # divmod(x, 1) should give (floor(x), 0) + x = QuadPrecision("42.7") + quotient, remainder = np.divmod(x, QuadPrecision("1.0")) + np.testing.assert_allclose(float(quotient), 42.0, rtol=1e-30) + np.testing.assert_allclose(float(remainder), 0.7, rtol=1e-14) + + # divmod(0, non-zero) should give (0, 0) + quotient, remainder = np.divmod(QuadPrecision("0.0"), QuadPrecision("5.0")) + assert float(quotient) == 0.0 + assert float(remainder) == 0.0 + + # divmod by 0 gives (inf, NaN) for positive dividend + quotient, remainder = np.divmod(QuadPrecision("1.0"), QuadPrecision("0.0")) + assert np.isinf(float(quotient)) and float(quotient) > 0 + assert np.isnan(float(remainder)) + + quotient, remainder = np.divmod(QuadPrecision("-1.0"), QuadPrecision("0.0")) + assert np.isinf(float(quotient)) and float(quotient) < 0 + assert np.isnan(float(remainder)) + + # divmod(inf, finite) gives (NaN, NaN) + quotient, remainder = np.divmod(QuadPrecision("inf"), QuadPrecision("5.0")) + assert np.isnan(float(quotient)) + assert np.isnan(float(remainder)) + + # divmod(finite, inf) gives (0, dividend) + quotient, remainder = np.divmod(QuadPrecision("5.0"), QuadPrecision("inf")) + np.testing.assert_allclose(float(quotient), 0.0, rtol=1e-30) + np.testing.assert_allclose(float(remainder), 5.0, rtol=1e-30) + + # Verify equivalence with floor_divide and mod + a = QuadPrecision("10.5") + b = QuadPrecision("3.2") + quotient, remainder = np.divmod(a, b) + expected_quotient = np.floor_divide(a, b) + expected_remainder = np.mod(a, b) + np.testing.assert_allclose(float(quotient), float(expected_quotient), rtol=1e-30) + np.testing.assert_allclose(float(remainder), float(expected_remainder), rtol=1e-30) + + +def test_divmod_array(): + """Test divmod with arrays""" + a = np.array([10.5, 21.0, -7.5, 0.0], dtype=QuadPrecDType()) + b = np.array([3.2, 4.0, 2.5, 5.0], dtype=QuadPrecDType()) + + quotients, remainders = np.divmod(a, b) + + # Check dtype + assert quotients.dtype.name == "QuadPrecDType128" + assert remainders.dtype.name == "QuadPrecDType128" + + # Check against NumPy float64 + a_float = np.array([10.5, 21.0, -7.5, 0.0], dtype=np.float64) + b_float = np.array([3.2, 4.0, 2.5, 5.0], dtype=np.float64) + expected_quotients, expected_remainders = np.divmod(a_float, b_float) + + for i in range(len(a)): + np.testing.assert_allclose( + float(quotients[i]), expected_quotients[i], rtol=1e-10, atol=1e-15, + err_msg=f"Quotient mismatch at index {i}" + ) + np.testing.assert_allclose( + float(remainders[i]), expected_remainders[i], rtol=1e-10, atol=1e-15, + err_msg=f"Remainder mismatch at index {i}" + ) + + +def test_divmod_broadcasting(): + """Test divmod with broadcasting""" + # Scalar with array + a = np.array([10.5, 21.0, 31.5], dtype=QuadPrecDType()) + b = QuadPrecision("3.0") + + quotients, remainders = np.divmod(a, b) + + assert quotients.dtype.name == "QuadPrecDType128" + assert remainders.dtype.name == "QuadPrecDType128" + assert len(quotients) == 3 + assert len(remainders) == 3 + + # Check values + expected_quotients = [3.0, 7.0, 10.0] + expected_remainders = [1.5, 0.0, 1.5] + + for i in range(3): + np.testing.assert_allclose(float(quotients[i]), expected_quotients[i], rtol=1e-14) + np.testing.assert_allclose(float(remainders[i]), expected_remainders[i], rtol=1e-14) + + @pytest.mark.parametrize("op", ["sinh", "cosh", "tanh", "arcsinh", "arccosh", "arctanh"]) @pytest.mark.parametrize("val", [ # Basic cases