1 /* 2 * Double-precision cosh(x) function. 3 * 4 * Copyright (c) 2022-2023, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 8 #include "math_config.h" 9 #include "pl_sig.h" 10 #include "pl_test.h" 11 12 #define AbsMask 0x7fffffffffffffff 13 #define SpecialBound \ 14 0x40861da04cbafe44 /* 0x1.61da04cbafe44p+9, above which exp overflows. */ 15 16 double 17 __exp_dd (double, double); 18 19 static double 20 specialcase (double x, uint64_t iax) 21 { 22 if (iax == 0x7ff0000000000000) 23 return INFINITY; 24 if (iax > 0x7ff0000000000000) 25 return __math_invalid (x); 26 /* exp overflows above SpecialBound. At this magnitude cosh(x) is dominated by 27 exp(x), so we can approximate cosh(x) by (exp(|x|/2)) ^ 2 / 2. */ 28 double t = __exp_dd (asdouble (iax) / 2, 0); 29 return (0.5 * t) * t; 30 } 31 32 /* Approximation for double-precision cosh(x). 33 cosh(x) = (exp(x) + exp(-x)) / 2. 34 The greatest observed error is in the special region, 1.93 ULP: 35 cosh(0x1.628af341989dap+9) got 0x1.fdf28623ef921p+1021 36 want 0x1.fdf28623ef923p+1021. 37 38 The greatest observed error in the non-special region is 1.03 ULP: 39 cosh(0x1.502cd8e56ab3bp+0) got 0x1.fe54962842d0ep+0 40 want 0x1.fe54962842d0fp+0. */ 41 double 42 cosh (double x) 43 { 44 uint64_t ix = asuint64 (x); 45 uint64_t iax = ix & AbsMask; 46 47 /* exp overflows a little bit before cosh, so use special-case handler for the 48 gap, as well as special values. */ 49 if (unlikely (iax >= SpecialBound)) 50 return specialcase (x, iax); 51 52 double ax = asdouble (iax); 53 /* Use double-precision exp helper to calculate exp(x), then: 54 cosh(x) = exp(|x|) / 2 + 1 / (exp(|x| * 2). */ 55 double t = __exp_dd (ax, 0); 56 return 0.5 * t + 0.5 / t; 57 } 58 59 PL_SIG (S, D, 1, cosh, -10.0, 10.0) 60 PL_TEST_ULP (cosh, 1.43) 61 PL_TEST_SYM_INTERVAL (cosh, 0, 0x1.61da04cbafe44p+9, 100000) 62 PL_TEST_SYM_INTERVAL (cosh, 0x1.61da04cbafe44p+9, 0x1p10, 1000) 63 PL_TEST_SYM_INTERVAL (cosh, 0x1p10, inf, 100) 64