From e0f48f48dbf586b7590a5d677063f2fee6be6410 Mon Sep 17 00:00:00 2001 From: Juniper Tyree <50025784+juntyr@users.noreply.github.com> Date: Fri, 8 Aug 2025 10:33:49 +0000 Subject: [PATCH] Implement hyperbolic ufuncs Add hyperbolic ufunc tests Add more tests for -0.0 Guarantee that min/max are zero-sign-sensitive --- quaddtype/numpy_quaddtype/src/ops.hpp | 116 ++++++++++++++++-- .../numpy_quaddtype/src/umath/unary_ops.cpp | 18 +++ quaddtype/release_tracker.md | 12 +- quaddtype/tests/test_quaddtype.py | 115 +++++++++++++++-- 4 files changed, 235 insertions(+), 26 deletions(-) diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp index ae38dbed..e103c8d8 100644 --- a/quaddtype/numpy_quaddtype/src/ops.hpp +++ b/quaddtype/numpy_quaddtype/src/ops.hpp @@ -151,6 +151,42 @@ quad_atan(const Sleef_quad *op) return Sleef_atanq1_u10(*op); } +static inline Sleef_quad +quad_sinh(const Sleef_quad *op) +{ + return Sleef_sinhq1_u10(*op); +} + +static inline Sleef_quad +quad_cosh(const Sleef_quad *op) +{ + return Sleef_coshq1_u10(*op); +} + +static inline Sleef_quad +quad_tanh(const Sleef_quad *op) +{ + return Sleef_tanhq1_u10(*op); +} + +static inline Sleef_quad +quad_asinh(const Sleef_quad *op) +{ + return Sleef_asinhq1_u10(*op); +} + +static inline Sleef_quad +quad_acosh(const Sleef_quad *op) +{ + return Sleef_acoshq1_u10(*op); +} + +static inline Sleef_quad +quad_atanh(const Sleef_quad *op) +{ + return Sleef_atanhq1_u10(*op); +} + // Unary long double operations typedef long double (*unary_op_longdouble_def)(const long double *); @@ -299,6 +335,42 @@ ld_atan(const long double *op) return atanl(*op); } +static inline long double +ld_sinh(const long double *op) +{ + return sinhl(*op); +} + +static inline long double +ld_cosh(const long double *op) +{ + return coshl(*op); +} + +static inline long double +ld_tanh(const long double *op) +{ + return tanhl(*op); +} + +static inline long double +ld_asinh(const long double *op) +{ + return asinhl(*op); +} + +static inline long double +ld_acosh(const long double *op) +{ + return acoshl(*op); +} + +static inline long double +ld_atanh(const long double *op) +{ + return atanhl(*op); +} + // Unary Quad properties typedef npy_bool (*unary_prop_quad_def)(const Sleef_quad *); @@ -442,33 +514,53 @@ quad_mod(const Sleef_quad *a, const Sleef_quad *b) static inline Sleef_quad quad_minimum(const Sleef_quad *in1, const Sleef_quad *in2) { - return Sleef_iunordq1(*in1, *in2) ? (Sleef_iunordq1(*in1, *in1) ? *in1 : *in2) - : Sleef_icmpleq1(*in1, *in2) ? *in1 - : *in2; + if (Sleef_iunordq1(*in1, *in2)) { + return Sleef_iunordq1(*in1, *in1) ? *in1 : *in2; + } + // minimum(-0.0, +0.0) = -0.0 + if (Sleef_icmpeqq1(*in1, QUAD_ZERO) && Sleef_icmpeqq1(*in2, QUAD_ZERO)) { + return Sleef_icmpleq1(Sleef_copysignq1(QUAD_ONE, *in1), Sleef_copysignq1(QUAD_ONE, *in2)) ? *in1 : *in2; + } + return Sleef_fminq1(*in1, *in2); } static inline Sleef_quad quad_maximum(const Sleef_quad *in1, const Sleef_quad *in2) { - return Sleef_iunordq1(*in1, *in2) ? (Sleef_iunordq1(*in1, *in1) ? *in1 : *in2) - : Sleef_icmpgeq1(*in1, *in2) ? *in1 - : *in2; + if (Sleef_iunordq1(*in1, *in2)) { + return Sleef_iunordq1(*in1, *in1) ? *in1 : *in2; + } + // maximum(-0.0, +0.0) = +0.0 + if (Sleef_icmpeqq1(*in1, QUAD_ZERO) && Sleef_icmpeqq1(*in2, QUAD_ZERO)) { + return Sleef_icmpgeq1(Sleef_copysignq1(QUAD_ONE, *in1), Sleef_copysignq1(QUAD_ONE, *in2)) ? *in1 : *in2; + } + return Sleef_fmaxq1(*in1, *in2); } static inline Sleef_quad quad_fmin(const Sleef_quad *in1, const Sleef_quad *in2) { - return Sleef_iunordq1(*in1, *in2) ? (Sleef_iunordq1(*in2, *in2) ? *in1 : *in2) - : Sleef_icmpleq1(*in1, *in2) ? *in1 - : *in2; + if (Sleef_iunordq1(*in1, *in2)) { + return Sleef_iunordq1(*in2, *in2) ? *in1 : *in2; + } + // fmin(-0.0, +0.0) = -0.0 + if (Sleef_icmpeqq1(*in1, QUAD_ZERO) && Sleef_icmpeqq1(*in2, QUAD_ZERO)) { + return Sleef_icmpleq1(Sleef_copysignq1(QUAD_ONE, *in1), Sleef_copysignq1(QUAD_ONE, *in2)) ? *in1 : *in2; + } + return Sleef_fminq1(*in1, *in2); } static inline Sleef_quad quad_fmax(const Sleef_quad *in1, const Sleef_quad *in2) { - return Sleef_iunordq1(*in1, *in2) ? (Sleef_iunordq1(*in2, *in2) ? *in1 : *in2) - : Sleef_icmpgeq1(*in1, *in2) ? *in1 - : *in2; + if (Sleef_iunordq1(*in1, *in2)) { + return Sleef_iunordq1(*in2, *in2) ? *in1 : *in2; + } + // maximum(-0.0, +0.0) = +0.0 + if (Sleef_icmpeqq1(*in1, QUAD_ZERO) && Sleef_icmpeqq1(*in2, QUAD_ZERO)) { + return Sleef_icmpgeq1(Sleef_copysignq1(QUAD_ONE, *in1), Sleef_copysignq1(QUAD_ONE, *in2)) ? *in1 : *in2; + } + return Sleef_fmaxq1(*in1, *in2); } static inline Sleef_quad diff --git a/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp b/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp index 09e45e97..d60591e8 100644 --- a/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/unary_ops.cpp @@ -216,5 +216,23 @@ init_quad_unary_ops(PyObject *numpy) if (create_quad_unary_ufunc(numpy, "arctan") < 0) { return -1; } + if (create_quad_unary_ufunc(numpy, "sinh") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "cosh") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "tanh") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "arcsinh") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "arccosh") < 0) { + return -1; + } + if (create_quad_unary_ufunc(numpy, "arctanh") < 0) { + return -1; + } return 0; } \ No newline at end of file diff --git a/quaddtype/release_tracker.md b/quaddtype/release_tracker.md index 9cff305d..76186f6c 100644 --- a/quaddtype/release_tracker.md +++ b/quaddtype/release_tracker.md @@ -50,12 +50,12 @@ | arctan | ✅ | ❌ _Need: basic tests + edge cases (NaN/inf/0/asymptotes)_ | | arctan2 | ✅ | ❌ _Need: basic tests + edge cases (NaN/inf/0/quadrant coverage)_ | | hypot | | | -| sinh | | | -| cosh | | | -| tanh | | | -| arcsinh | | | -| arccosh | | | -| arctanh | | | +| sinh | ✅ | ✅ | +| cosh | ✅ | ✅ | +| tanh | ✅ | ✅ | +| arcsinh | ✅ | ✅ | +| arccosh | ✅ | ✅ | +| arctanh | ✅ | ✅ | | degrees | | | | radians | | | | deg2rad | | | diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index 1b4591cb..25550db8 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -35,22 +35,28 @@ def test_basic_equality(): @pytest.mark.parametrize("op", ["add", "sub", "mul", "truediv", "pow", "copysign"]) -@pytest.mark.parametrize("other", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"]) -def test_binary_ops(op, other): - if op == "truediv" and float(other) == 0: +@pytest.mark.parametrize("a", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"]) +@pytest.mark.parametrize("b", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"]) +def test_binary_ops(op, a, b): + if op == "truediv" and float(b) == 0: pytest.xfail("float division by zero") op_func = getattr(operator, op, None) or getattr(np, op) - quad_a = QuadPrecision("12.5") - quad_b = QuadPrecision(other) - float_a = 12.5 - float_b = float(other) + quad_a = QuadPrecision(a) + quad_b = QuadPrecision(b) + float_a = float(a) + float_b = float(b) quad_result = op_func(quad_a, quad_b) float_result = op_func(float_a, float_b) np.testing.assert_allclose(np.float64(quad_result), float_result, atol=1e-10, rtol=0, equal_nan=True) + # Check sign for zero results + if float_result == 0.0: + assert np.signbit(float_result) == np.signbit( + quad_result), f"Zero sign mismatch for {op}({a}, {b})" + @pytest.mark.parametrize("op", ["eq", "ne", "le", "lt", "ge", "gt"]) @pytest.mark.parametrize("a", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"]) @@ -91,8 +97,20 @@ def test_array_minmax(op, a, b): quad_res = op_func(quad_a, quad_b) float_res = op_func(float_a, float_b) + # native implementation may not be sensitive to zero signs + # but we want to enforce it for the quad dtype + # e.g. min(+0.0, -0.0) = -0.0 + if float_a == 0.0 and float_b == 0.0: + assert float_res == 0.0 + float_res = np.copysign(0.0, op_func(np.copysign(1.0, float_a), np.copysign(1.0, float_b))) + np.testing.assert_array_equal(quad_res.astype(float), float_res) + # Check sign for zero results + if float_res == 0.0: + assert np.signbit(float_res) == np.signbit( + quad_res), f"Zero sign mismatch for {op}({a}, {b})" + @pytest.mark.parametrize("op", ["amin", "amax", "nanmin", "nanmax"]) @pytest.mark.parametrize("a", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"]) @@ -105,8 +123,20 @@ def test_array_aminmax(op, a, b): quad_res = op_func(quad_ab) float_res = op_func(float_ab) + # native implementation may not be sensitive to zero signs + # but we want to enforce it for the quad dtype + # e.g. min(+0.0, -0.0) = -0.0 + if float(a) == 0.0 and float(b) == 0.0: + assert float_res == 0.0 + float_res = np.copysign(0.0, op_func(np.array([np.copysign(1.0, float(a)), np.copysign(1.0, float(b))]))) + np.testing.assert_array_equal(np.array(quad_res).astype(float), float_res) + # Check sign for zero results + if float_res == 0.0: + assert np.signbit(float_res) == np.signbit( + quad_res), f"Zero sign mismatch for {op}({a}, {b})" + @pytest.mark.parametrize("op", ["negative", "positive", "absolute", "sign", "signbit", "isfinite", "isinf", "isnan", "sqrt", "square", "reciprocal"]) @pytest.mark.parametrize("val", ["3.0", "-3.0", "12.5", "100.0", "1e100", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"]) @@ -126,7 +156,7 @@ def test_unary_ops(op, val): np.testing.assert_array_equal(np.array(quad_result).astype(float), float_result) - if op in ["negative", "positive", "absolute", "sign"]: + if (float_result == 0.0) and (op not in ["signbit", "isfinite", "isinf", "isnan"]): assert np.signbit(float_result) == np.signbit(quad_result) @@ -290,6 +320,11 @@ def test_logarithmic_functions(op, val): np.testing.assert_allclose(float(quad_result), float_result, rtol=rtol, atol=atol, err_msg=f"Value mismatch for {op}({val})") + # Check sign for zero results + if float_result == 0.0: + assert np.signbit(float_result) == np.signbit( + quad_result), f"Zero sign mismatch for {op}({a}, {b})" + @pytest.mark.parametrize("val", [ # Basic cases around -1 (critical point for log1p) @@ -304,6 +339,8 @@ def test_logarithmic_functions(op, val): "-1.1", "-2.0", "-10.0", # Large positive values "1e10", "1e15", "1e100", + # Edge cases + "0.0", "-0.0", # Special values "inf", "-inf", "nan", "-nan" ]) @@ -341,9 +378,16 @@ def test_log1p(val): np.testing.assert_allclose(float(quad_result), float_result, rtol=rtol, atol=atol, err_msg=f"Value mismatch for log1p({val})") + # Check sign for zero results + if float_result == 0.0: + assert np.signbit(float_result) == np.signbit( + quad_result), f"Zero sign mismatch for {op}({val})" + def test_inf(): assert QuadPrecision("inf") > QuadPrecision("1e1000") + assert np.signbit(QuadPrecision("inf")) == 0 assert QuadPrecision("-inf") < QuadPrecision("-1e1000") + assert np.signbit(QuadPrecision("-inf")) == 1 def test_dtype_creation(): @@ -448,3 +492,58 @@ def test_mod(a, b, backend, op): numpy_negative = numpy_result < 0 assert result_negative == numpy_negative, f"Sign mismatch for {a} % {b}: quad={result_negative}, numpy={numpy_negative}" + + +@pytest.mark.parametrize("op", ["sinh", "cosh", "tanh", "arcsinh", "arccosh", "arctanh"]) +@pytest.mark.parametrize("val", [ + # Basic cases + "0.0", "-0.0", "1.0", "-1.0", "2.0", "-2.0", + # Small values + "1e-10", "-1e-10", "1e-15", "-1e-15", + # Values near one + "0.9", "-0.9", "0.9999", "-0.9999", + "1.1", "-1.1", "1.0001", "-1.0001", + # Medium values + "10.0", "-10.0", "20.0", "-20.0", + # Large values + "100.0", "200.0", "700.0", "1000.0", "1e100", "1e308", + "-100.0", "-200.0", "-700.0", "-1000.0", "-1e100", "-1e308", + # Fractional values + "0.5", "-0.5", "1.5", "-1.5", "2.5", "-2.5", + # Special values + "inf", "-inf", "nan", "-nan" +]) +def test_hyperbolic_functions(op, val): + """Comprehensive test for hyperbolic functions: sinh, cosh, tanh, arcsinh, arccosh, arctanh""" + op_func = getattr(np, op) + + quad_val = QuadPrecision(val) + float_val = float(val) + + quad_result = op_func(quad_val) + float_result = op_func(float_val) + + # Handle NaN cases + if np.isnan(float_result): + assert np.isnan( + float(quad_result)), f"Expected NaN for {op}({val}), got {float(quad_result)}" + return + + # Handle infinity cases + if np.isinf(float_result): + assert np.isinf( + float(quad_result)), f"Expected inf for {op}({val}), got {float(quad_result)}" + assert np.sign(float_result) == np.sign( + float(quad_result)), f"Infinity sign mismatch for {op}({val})" + return + + # For finite non-zero results + # Use relative tolerance for exponential functions due to their rapid growth + rtol = 1e-13 if abs(float_result) < 1e100 else 1e-10 + np.testing.assert_allclose(float(quad_result), float_result, rtol=rtol, atol=1e-15, + err_msg=f"Value mismatch for {op}({val})") + + # Check sign for zero results + if float_result == 0.0: + assert np.signbit(float_result) == np.signbit( + quad_result), f"Zero sign mismatch for {op}({val})"