diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index 11c56e3..5005320 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/quaddtype/numpy_quaddtype/src/ops.hpp @@ -467,6 +467,39 @@ quad_div(const Sleef_quad *a, const Sleef_quad *b) return Sleef_divq1_u05(*a, *b); } +static inline Sleef_quad +quad_floor_divide(const Sleef_quad *a, const Sleef_quad *b) +{ + // Handle NaN inputs + if (Sleef_iunordq1(*a, *b)) { + return Sleef_iunordq1(*a, *a) ? *a : *b; + } + + // inf / finite_nonzero or -inf / finite_nonzero -> NaN + // But inf / 0 -> inf + if (quad_isinf(a) && quad_isfinite(b) && !Sleef_icmpeqq1(*b, QUAD_ZERO)) { + return QUAD_NAN; + } + + // 0 / 0 (including -0.0 / 0.0, 0.0 / -0.0, -0.0 / -0.0) -> NaN + if (Sleef_icmpeqq1(*a, QUAD_ZERO) && Sleef_icmpeqq1(*b, QUAD_ZERO)) { + return QUAD_NAN; + } + + Sleef_quad quotient = Sleef_divq1_u05(*a, *b); + Sleef_quad result = Sleef_floorq1(quotient); + + // floor_divide semantics: when result is -0.0 from non-zero numerator, convert to -1.0 + // This happens when: (negative & non-zero)/+inf, (positive & non-zero)/-inf + // But NOT when numerator is ±0.0 (then result stays as ±0.0) + if (Sleef_icmpeqq1(result, QUAD_ZERO) && quad_signbit(&result) && + !Sleef_icmpeqq1(*a, QUAD_ZERO)) { + return Sleef_negq1(QUAD_ONE); // -1.0 + } + + return result; +} + static inline Sleef_quad quad_pow(const Sleef_quad *a, const Sleef_quad *b) { @@ -696,6 +729,38 @@ ld_div(const long double *a, const long double *b) return (*a) / (*b); } +static inline long double +ld_floor_divide(const long double *a, const long double *b) +{ + // Handle NaN inputs + if (isnan(*a) || isnan(*b)) { + return isnan(*a) ? *a : *b; + } + + // inf / finite_nonzero or -inf / finite_nonzero -> NaN + // But inf / 0 -> inf + if (isinf(*a) && isfinite(*b) && *b != 0.0L) { + return NAN; + } + + // 0 / 0 (including -0.0 / 0.0, 0.0 / -0.0, -0.0 / -0.0) -> NaN + if (*a == 0.0L && *b == 0.0L) { + return NAN; + } + + // Compute a / b and apply floor + long double result = floorl((*a) / (*b)); + + // floor_divide semantics: when result is -0.0 from non-zero numerator, convert to -1.0 + // This happens when: (negative & non-zero)/+inf, (positive & non-zero)/-inf + // But NOT when numerator is ±0.0 (then result stays as ±0.0) + if (result == 0.0L && signbit(result) && *a != 0.0L) { + return -1.0L; + } + + return result; +} + static inline long double ld_pow(const long double *a, const long double *b) { diff --git a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp index b6b91d9..45735ea 100644 --- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp @@ -218,6 +218,9 @@ init_quad_binary_ops(PyObject *numpy) } // Note: true_divide is an alias to divide in NumPy for floating-point types // No need to register separately + if (create_quad_binary_ufunc(numpy, "floor_divide") < 0) { + return -1; + } if (create_quad_binary_ufunc(numpy, "power") < 0) { return -1; } diff --git a/quaddtype/release_tracker.md b/quaddtype/release_tracker.md index edcac8e..e28f1c9 100644 --- a/quaddtype/release_tracker.md +++ b/quaddtype/release_tracker.md @@ -13,7 +13,7 @@ | logaddexp | ✅ | ✅ | | logaddexp2 | ✅ | ✅ | | true_divide | ✅ | ✅ | -| floor_divide | | | +| floor_divide | ✅ | ✅ | | negative | ✅ | ✅ | | positive | ✅ | ✅ | | power | ✅ | ✅ | diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index 758a40d..292c895 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -732,6 +732,111 @@ def test_true_divide_special_properties(): np.testing.assert_allclose(float(result), 3.0, rtol=1e-30) +@pytest.mark.parametrize( + "x_val", + [ + 0.0, 1.0, 2.0, -1.0, -2.0, + 0.5, -0.5, + 100.0, 1000.0, -100.0, -1000.0, + 1e-10, -1e-10, 1e-20, -1e-20, + float("inf"), float("-inf"), float("nan"), float("-nan"), -0.0 + ] +) +@pytest.mark.parametrize( + "y_val", + [ + 0.0, 1.0, 2.0, -1.0, -2.0, + 0.5, -0.5, + 100.0, 1000.0, -100.0, -1000.0, + 1e-10, -1e-10, 1e-20, -1e-20, + float("inf"), float("-inf"), float("nan"), float("-nan"), -0.0 + ] +) +def test_floor_divide(x_val, y_val): + """Test floor_divide ufunc with comprehensive edge cases""" + x_quad = QuadPrecision(str(x_val)) + y_quad = QuadPrecision(str(y_val)) + + # Compute using QuadPrecision + result_quad = np.floor_divide(x_quad, y_quad) + + # Compute using float64 for comparison + result_float64 = np.floor_divide(np.float64(x_val), np.float64(y_val)) + + # Compare results + if np.isnan(result_float64): + assert np.isnan(float(result_quad)), f"Expected NaN for floor_divide({x_val}, {y_val})" + elif np.isinf(result_float64): + assert np.isinf(float(result_quad)), f"Expected inf for floor_divide({x_val}, {y_val})" + assert np.sign(float(result_quad)) == np.sign(result_float64), f"Sign mismatch for floor_divide({x_val}, {y_val})" + else: + # For finite results, check relative tolerance + # Use absolute tolerance for large numbers due to float64 precision limits + atol = max(1e-10, abs(result_float64) * 1e-9) if abs(result_float64) > 1e6 else 1e-10 + np.testing.assert_allclose( + float(result_quad), result_float64, rtol=1e-12, atol=atol, + err_msg=f"Mismatch for floor_divide({x_val}, {y_val})" + ) +def test_floor_divide_special_properties(): + """Test special mathematical properties of floor_divide""" + # floor_divide(x, 1) = floor(x) + x = QuadPrecision("42.7") + result = np.floor_divide(x, QuadPrecision("1.0")) + np.testing.assert_allclose(float(result), 42.0, rtol=1e-30) + + # floor_divide(0, non-zero) = 0 + result = np.floor_divide(QuadPrecision("0.0"), QuadPrecision("5.0")) + assert float(result) == 0.0 + + # floor_divide by 0 gives inf (with appropriate sign) + result = np.floor_divide(QuadPrecision("1.0"), QuadPrecision("0.0")) + assert np.isinf(float(result)) and float(result) > 0 + + result = np.floor_divide(QuadPrecision("-1.0"), QuadPrecision("0.0")) + assert np.isinf(float(result)) and float(result) < 0 + + # 0 / 0 = NaN + result = np.floor_divide(QuadPrecision("0.0"), QuadPrecision("0.0")) + assert np.isnan(float(result)) + + # inf / inf = NaN + result = np.floor_divide(QuadPrecision("inf"), QuadPrecision("inf")) + assert np.isnan(float(result)) + + # inf / finite_nonzero = NaN (NumPy behavior) + result = np.floor_divide(QuadPrecision("inf"), QuadPrecision("100.0")) + assert np.isnan(float(result)) + + # finite / inf = 0 + result = np.floor_divide(QuadPrecision("100.0"), QuadPrecision("inf")) + assert float(result) == 0.0 + + # floor_divide rounds toward negative infinity + result = np.floor_divide(QuadPrecision("7.0"), QuadPrecision("3.0")) + assert float(result) == 2.0 # floor(7/3) = floor(2.333...) = 2 + + result = np.floor_divide(QuadPrecision("-7.0"), QuadPrecision("3.0")) + assert float(result) == -3.0 # floor(-7/3) = floor(-2.333...) = -3 + + result = np.floor_divide(QuadPrecision("7.0"), QuadPrecision("-3.0")) + assert float(result) == -3.0 # floor(7/-3) = floor(-2.333...) = -3 + + result = np.floor_divide(QuadPrecision("-7.0"), QuadPrecision("-3.0")) + assert float(result) == 2.0 # floor(-7/-3) = floor(2.333...) = 2 + + # floor_divide(x, x) = 1 for positive finite non-zero x + x = QuadPrecision("7.123456789") + result = np.floor_divide(x, x) + np.testing.assert_allclose(float(result), 1.0, rtol=1e-30) + + # Relationship with floor and true_divide + x = QuadPrecision("10.5") + y = QuadPrecision("3.2") + result_floor_divide = np.floor_divide(x, y) + result_floor_true_divide = np.floor(np.true_divide(x, y)) + np.testing.assert_allclose(float(result_floor_divide), float(result_floor_true_divide), rtol=1e-30) + + def test_inf(): assert QuadPrecision("inf") > QuadPrecision("1e1000") assert np.signbit(QuadPrecision("inf")) == 0