From 0a26cfc6d8dea2cfed7a5823f3f0cc73ba1c43a3 Mon Sep 17 00:00:00 2001 From: swayaminsync Date: Thu, 16 Oct 2025 12:30:57 +0530 Subject: [PATCH 1/3] logaddexp2 ufunc --- quaddtype/numpy_quaddtype/src/ops.hpp | 82 ++++++++++++++ .../numpy_quaddtype/src/umath/binary_ops.cpp | 3 + quaddtype/release_tracker.md | 2 +- quaddtype/tests/test_quaddtype.py | 107 ++++++++++++++++++ 4 files changed, 193 insertions(+), 1 deletion(-) diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index 9fbf4b4a..948a5434 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/quaddtype/numpy_quaddtype/src/ops.hpp @@ -626,6 +626,49 @@ quad_logaddexp(const Sleef_quad *x, const Sleef_quad *y) return Sleef_addq1_u05(max_val, log1p_term); } +static inline Sleef_quad +quad_logaddexp2(const Sleef_quad *x, const Sleef_quad *y) +{ + // logaddexp2(x, y) = log2(2^x + 2^y) + // Numerically stable implementation: max(x, y) + log2(1 + 2^(-abs(x - y))) + + // Handle NaN + if (Sleef_iunordq1(*x, *y)) { + return Sleef_iunordq1(*x, *x) ? *x : *y; + } + + // Handle infinities + // If both are -inf, result is -inf + Sleef_quad neg_inf = Sleef_negq1(QUAD_POS_INF); + if (Sleef_icmpeqq1(*x, neg_inf) && Sleef_icmpeqq1(*y, neg_inf)) { + return neg_inf; + } + + // If either is +inf, result is +inf + if (Sleef_icmpeqq1(*x, QUAD_POS_INF) || Sleef_icmpeqq1(*y, QUAD_POS_INF)) { + return QUAD_POS_INF; + } + + // If one is -inf, result is the other value + if (Sleef_icmpeqq1(*x, neg_inf)) { + return *y; + } + if (Sleef_icmpeqq1(*y, neg_inf)) { + return *x; + } + + // log2(2^x + 2^y) = max(x, y) + log2(1 + 2^(-abs(x - y))) + Sleef_quad diff = Sleef_subq1_u05(*x, *y); + Sleef_quad abs_diff = Sleef_fabsq1(diff); + Sleef_quad neg_abs_diff = Sleef_negq1(abs_diff); + Sleef_quad exp2_term = Sleef_exp2q1_u10(neg_abs_diff); + Sleef_quad one_plus_exp2 = Sleef_addq1_u05(QUAD_ONE, exp2_term); + Sleef_quad log2_term = Sleef_log2q1_u10(one_plus_exp2); + + Sleef_quad max_val = Sleef_icmpgtq1(*x, *y) ? *x : *y; + return Sleef_addq1_u05(max_val, log2_term); +} + // Binary long double operations typedef long double (*binary_op_longdouble_def)(const long double *, const long double *); @@ -759,6 +802,45 @@ ld_logaddexp(const long double *x, const long double *y) return max_val + log1pl(expl(-abs_diff)); } +static inline long double +ld_logaddexp2(const long double *x, const long double *y) +{ + // logaddexp2(x, y) = log2(2^x + 2^y) + // Numerically stable implementation: max(x, y) + log2(1 + 2^(-abs(x - y))) + + // Handle NaN + if (isnan(*x) || isnan(*y)) { + return isnan(*x) ? *x : *y; + } + + // Handle infinities + // If both are -inf, result is -inf + if (isinf(*x) && *x < 0 && isinf(*y) && *y < 0) { + return -INFINITY; + } + + // If either is +inf, result is +inf + if ((isinf(*x) && *x > 0) || (isinf(*y) && *y > 0)) { + return INFINITY; + } + + // If one is -inf, result is the other value + if (isinf(*x) && *x < 0) { + return *y; + } + if (isinf(*y) && *y < 0) { + return *x; + } + + // Numerically stable computation + // log2(2^x + 2^y) = max(x, y) + log2(1 + 2^(-abs(x - y))) + long double diff = *x - *y; + long double abs_diff = fabsl(diff); + long double max_val = (*x > *y) ? *x : *y; + // log2(1 + z) = log(1 + z) / log(2) + return max_val + log1pl(exp2l(-abs_diff)) / M_LN2; +} + // comparison quad functions typedef npy_bool (*cmp_quad_def)(const Sleef_quad *, const Sleef_quad *); diff --git a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp index 5a24e678..fee2bc74 100644 --- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp @@ -243,5 +243,8 @@ init_quad_binary_ops(PyObject *numpy) if (create_quad_binary_ufunc(numpy, "logaddexp") < 0) { return -1; } + if (create_quad_binary_ufunc(numpy, "logaddexp2") < 0) { + return -1; + } return 0; } \ No newline at end of file diff --git a/quaddtype/release_tracker.md b/quaddtype/release_tracker.md index 1dc4de36..26fb871b 100644 --- a/quaddtype/release_tracker.md +++ b/quaddtype/release_tracker.md @@ -11,7 +11,7 @@ | matmul | ✅ | ✅ | | divide | ✅ | ✅ | | logaddexp | ✅ | ✅ | -| logaddexp2 | | | +| logaddexp2 | ✅ | ✅ | | true_divide | | | | floor_divide | | | | negative | ✅ | ✅ | diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index a669cdd4..f1741552 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -528,6 +528,113 @@ def test_logaddexp_special_properties(): np.testing.assert_allclose(float(result1), float(result2), rtol=1e-14) +@pytest.mark.parametrize("x", [ + # Regular values + "0.0", "1.0", "2.0", "-1.0", "-2.0", "0.5", "-0.5", + # Large values (test numerical stability) + "100.0", "1000.0", "-100.0", "-1000.0", + # Small values + "1e-10", "-1e-10", "1e-20", "-1e-20", + # Special values + "inf", "-inf", "nan", "-nan", "-0.0" +]) +@pytest.mark.parametrize("y", [ + # Regular values + "0.0", "1.0", "2.0", "-1.0", "-2.0", "0.5", "-0.5", + # Large values + "100.0", "1000.0", "-100.0", "-1000.0", + # Small values + "1e-10", "-1e-10", "1e-20", "-1e-20", + # Special values + "inf", "-inf", "nan", "-nan", "-0.0" +]) +def test_logaddexp2(x, y): + """Comprehensive test for logaddexp2 function: log2(2^x + 2^y)""" + quad_x = QuadPrecision(x) + quad_y = QuadPrecision(y) + float_x = float(x) + float_y = float(y) + + quad_result = np.logaddexp2(quad_x, quad_y) + float_result = np.logaddexp2(float_x, float_y) + + # Handle NaN cases + if np.isnan(float_result): + assert np.isnan(float(quad_result)), \ + f"Expected NaN for logaddexp2({x}, {y}), got {float(quad_result)}" + return + + # Handle infinity cases + if np.isinf(float_result): + assert np.isinf(float(quad_result)), \ + f"Expected inf for logaddexp2({x}, {y}), got {float(quad_result)}" + if not np.isnan(float_result): + assert np.sign(float_result) == np.sign(float(quad_result)), \ + f"Infinity sign mismatch for logaddexp2({x}, {y})" + return + + # For finite results, check with appropriate tolerance + # logaddexp2 is numerically sensitive, especially for large differences + if abs(float_x - float_y) > 50: + # When values differ greatly, result should be close to max(x, y) + rtol = 1e-10 + atol = 1e-10 + else: + rtol = 1e-13 + atol = 1e-15 + + np.testing.assert_allclose( + float(quad_result), float_result, + rtol=rtol, atol=atol, + err_msg=f"Value mismatch for logaddexp2({x}, {y})" + ) + + +def test_logaddexp2_special_properties(): + """Test special mathematical properties of logaddexp2""" + # logaddexp2(x, x) = x + 1 (since log2(2^x + 2^x) = log2(2 * 2^x) = log2(2) + log2(2^x) = 1 + x) + x = QuadPrecision("2.0") + result = np.logaddexp2(x, x) + expected = float(x) + 1.0 + np.testing.assert_allclose(float(result), expected, rtol=1e-14) + + # logaddexp2(x, -inf) = x + x = QuadPrecision("5.0") + result = np.logaddexp2(x, QuadPrecision("-inf")) + np.testing.assert_allclose(float(result), float(x), rtol=1e-14) + + # logaddexp2(-inf, x) = x + result = np.logaddexp2(QuadPrecision("-inf"), x) + np.testing.assert_allclose(float(result), float(x), rtol=1e-14) + + # logaddexp2(-inf, -inf) = -inf + result = np.logaddexp2(QuadPrecision("-inf"), QuadPrecision("-inf")) + assert np.isinf(float(result)) and float(result) < 0 + + # logaddexp2(inf, anything) = inf + result = np.logaddexp2(QuadPrecision("inf"), QuadPrecision("100.0")) + assert np.isinf(float(result)) and float(result) > 0 + + # logaddexp2(anything, inf) = inf + result = np.logaddexp2(QuadPrecision("100.0"), QuadPrecision("inf")) + assert np.isinf(float(result)) and float(result) > 0 + + # Commutativity: logaddexp2(x, y) = logaddexp2(y, x) + x = QuadPrecision("3.0") + y = QuadPrecision("5.0") + result1 = np.logaddexp2(x, y) + result2 = np.logaddexp2(y, x) + np.testing.assert_allclose(float(result1), float(result2), rtol=1e-14) + + # Relationship with logaddexp: logaddexp2(x, y) = logaddexp(x*ln2, y*ln2) / ln2 + x = QuadPrecision("2.0") + y = QuadPrecision("3.0") + result_logaddexp2 = np.logaddexp2(x, y) + ln2 = np.log(2.0) + result_logaddexp = np.logaddexp(float(x) * ln2, float(y) * ln2) / ln2 + np.testing.assert_allclose(float(result_logaddexp2), result_logaddexp, rtol=1e-13) + + def test_inf(): assert QuadPrecision("inf") > QuadPrecision("1e1000") assert np.signbit(QuadPrecision("inf")) == 0 From d4d5b4e3029de6237d1629f1fd27382deb7e7f72 Mon Sep 17 00:00:00 2001 From: swayaminsync Date: Thu, 16 Oct 2025 12:44:40 +0530 Subject: [PATCH 2/3] use log2l --- quaddtype/numpy_quaddtype/src/ops.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index 948a5434..11c56e3a 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/quaddtype/numpy_quaddtype/src/ops.hpp @@ -832,13 +832,12 @@ ld_logaddexp2(const long double *x, const long double *y) return *x; } - // Numerically stable computation // log2(2^x + 2^y) = max(x, y) + log2(1 + 2^(-abs(x - y))) long double diff = *x - *y; long double abs_diff = fabsl(diff); long double max_val = (*x > *y) ? *x : *y; - // log2(1 + z) = log(1 + z) / log(2) - return max_val + log1pl(exp2l(-abs_diff)) / M_LN2; + // Use native log2l function for base-2 logarithm + return max_val + log2l(1.0L + exp2l(-abs_diff)); } // comparison quad functions From eaf0812084d9647e472c5f6df517846c518ea330 Mon Sep 17 00:00:00 2001 From: swayaminsync Date: Thu, 16 Oct 2025 13:08:32 +0530 Subject: [PATCH 3/3] true_divide impl --- .../numpy_quaddtype/src/umath/binary_ops.cpp | 2 + quaddtype/release_tracker.md | 2 +- quaddtype/tests/test_quaddtype.py | 97 +++++++++++++++++++ 3 files changed, 100 insertions(+), 1 deletion(-) diff --git a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp index fee2bc74..b6b91d93 100644 --- a/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp @@ -216,6 +216,8 @@ init_quad_binary_ops(PyObject *numpy) if (create_quad_binary_ufunc(numpy, "divide") < 0) { return -1; } + // 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, "power") < 0) { return -1; } diff --git a/quaddtype/release_tracker.md b/quaddtype/release_tracker.md index 26fb871b..edcac8eb 100644 --- a/quaddtype/release_tracker.md +++ b/quaddtype/release_tracker.md @@ -12,7 +12,7 @@ | divide | ✅ | ✅ | | logaddexp | ✅ | ✅ | | logaddexp2 | ✅ | ✅ | -| true_divide | | | +| true_divide | ✅ | ✅ | | floor_divide | | | | negative | ✅ | ✅ | | positive | ✅ | ✅ | diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index f1741552..758a40df 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -635,6 +635,103 @@ def test_logaddexp2_special_properties(): np.testing.assert_allclose(float(result_logaddexp2), result_logaddexp, rtol=1e-13) +@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_true_divide(x_val, y_val): + """Test true_divide ufunc with comprehensive edge cases""" + x_quad = QuadPrecision(str(x_val)) + y_quad = QuadPrecision(str(y_val)) + + # Compute using QuadPrecision + result_quad = np.true_divide(x_quad, y_quad) + + # Compute using float64 for comparison + result_float64 = np.true_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 true_divide({x_val}, {y_val})" + elif np.isinf(result_float64): + assert np.isinf(float(result_quad)), f"Expected inf for true_divide({x_val}, {y_val})" + assert np.sign(float(result_quad)) == np.sign(result_float64), f"Sign mismatch for true_divide({x_val}, {y_val})" + else: + # For finite results, check relative tolerance + np.testing.assert_allclose( + float(result_quad), result_float64, rtol=1e-14, + err_msg=f"Mismatch for true_divide({x_val}, {y_val})" + ) + + +def test_true_divide_special_properties(): + """Test special mathematical properties of true_divide""" + # Division by 1 returns the original value + x = QuadPrecision("42.123456789") + result = np.true_divide(x, QuadPrecision("1.0")) + np.testing.assert_allclose(float(result), float(x), rtol=1e-30) + + # Division of 0 by any non-zero number is 0 + result = np.true_divide(QuadPrecision("0.0"), QuadPrecision("5.0")) + assert float(result) == 0.0 + + # Division by 0 gives inf (with appropriate sign) + result = np.true_divide(QuadPrecision("1.0"), QuadPrecision("0.0")) + assert np.isinf(float(result)) and float(result) > 0 + + result = np.true_divide(QuadPrecision("-1.0"), QuadPrecision("0.0")) + assert np.isinf(float(result)) and float(result) < 0 + + # 0 / 0 = NaN + result = np.true_divide(QuadPrecision("0.0"), QuadPrecision("0.0")) + assert np.isnan(float(result)) + + # inf / inf = NaN + result = np.true_divide(QuadPrecision("inf"), QuadPrecision("inf")) + assert np.isnan(float(result)) + + # inf / finite = inf + result = np.true_divide(QuadPrecision("inf"), QuadPrecision("100.0")) + assert np.isinf(float(result)) and float(result) > 0 + + # finite / inf = 0 + result = np.true_divide(QuadPrecision("100.0"), QuadPrecision("inf")) + assert float(result) == 0.0 + + # Self-division (x / x) = 1 for finite non-zero x + x = QuadPrecision("7.123456789") + result = np.true_divide(x, x) + np.testing.assert_allclose(float(result), 1.0, rtol=1e-30) + + # Sign preservation: (-x) / y = -(x / y) + x = QuadPrecision("5.5") + y = QuadPrecision("2.2") + result1 = np.true_divide(-x, y) + result2 = -np.true_divide(x, y) + np.testing.assert_allclose(float(result1), float(result2), rtol=1e-30) + + # Sign rule: negative / negative = positive + result = np.true_divide(QuadPrecision("-6.0"), QuadPrecision("-2.0")) + assert float(result) > 0 + np.testing.assert_allclose(float(result), 3.0, rtol=1e-30) + + def test_inf(): assert QuadPrecision("inf") > QuadPrecision("1e1000") assert np.signbit(QuadPrecision("inf")) == 0