logaddexp2 ufunc

Author: SwayamInSync

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 *);
 

quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp

@@ -243,5 +243,8 @@ init_quad_binary_ops(PyObject *numpy)
     if (create_quad_binary_ufunc<quad_logaddexp, ld_logaddexp>(numpy, "logaddexp") < 0) {
         return -1;
     }
+    if (create_quad_binary_ufunc<quad_logaddexp2, ld_logaddexp2>(numpy, "logaddexp2") < 0) {
+        return -1;
+    }
     return 0;
 }

quaddtype/release_tracker.md

@@ -11,7 +11,7 @@
 | matmul        | ✅    | ✅                                                                   |
 | divide        | ✅    | ✅                                                                   |
 | logaddexp     | ✅    | ✅                                                                   |
-| logaddexp2    |       |                                                                      |
+| logaddexp2    | ✅    | ✅                                                                   |
 | true_divide   |       |                                                                      |
 | floor_divide  |       |                                                                      |
 | negative      | ✅    | ✅                                                                   |

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