pxfm/hyperbolic/
cosh.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 7/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::common::{dd_fmla, f_fmla};
30use crate::double_double::DoubleDouble;
31use crate::exponents::{EXP_REDUCE_T0, EXP_REDUCE_T1};
32use crate::hyperbolic::acosh::lpoly_xd_generic;
33use crate::hyperbolic::sinh::hyperbolic_exp_accurate;
34
35#[cold]
36fn as_cosh_zero(x: f64) -> f64 {
37    static CH: [(u64, u64); 4] = [
38        (0xb90c7e8db669f624, 0x3fe0000000000000),
39        (0x3c45555555556135, 0x3fa5555555555555),
40        (0xbbef49f4a6e838f2, 0x3f56c16c16c16c17),
41        (0x3b3a4ffbe15316aa, 0x3efa01a01a01a01a),
42    ];
43    const CL: [u64; 4] = [
44        0x3e927e4fb7789f5c,
45        0x3e21eed8eff9089c,
46        0x3da939749ce13dad,
47        0x3d2ae9891efb6691,
48    ];
49
50    let dx2 = DoubleDouble::from_exact_mult(x, x);
51
52    let yw0 = f_fmla(dx2.hi, f64::from_bits(CL[3]), f64::from_bits(CL[2]));
53    let yw1 = f_fmla(dx2.hi, yw0, f64::from_bits(CL[1]));
54
55    let y2 = dx2.hi * f_fmla(dx2.hi, yw1, f64::from_bits(CL[0]));
56
57    let mut y1 = lpoly_xd_generic(dx2, CH, y2);
58    y1 = DoubleDouble::quick_mult(y1, dx2); // y2 = y1.l
59    let y0 = DoubleDouble::from_exact_add(1.0, y1.hi); // y0 = y0.hi
60    let mut p = DoubleDouble::from_exact_add(y0.lo, y1.lo);
61    let mut t = p.hi.to_bits();
62    if (t & 0x000fffffffffffff) == 0 {
63        let w = p.lo.to_bits();
64        if ((w ^ t) >> 63) != 0 {
65            t = t.wrapping_sub(1);
66        } else {
67            t = t.wrapping_add(1);
68        }
69        p.hi = f64::from_bits(t);
70    }
71    y0.hi + p.hi
72}
73
74/// Hyperbolic cosine function
75///
76/// Max ULP 0.5
77pub fn f_cosh(x: f64) -> f64 {
78    /*
79     The function sinh(x) is approximated by a minimax polynomial
80     cosh(x)~1+x^2*P(x^2) for |x|<0.125. For other arguments the
81     identity cosh(x)=(exp(|x|)+exp(-|x|))/2 is used. For |x|<5 both
82     exponents are calculated with slightly higher precision than
83     double. For 5<|x|<36.736801 the exp(-|x|) is rather small and is
84     calculated with double precision but exp(|x|) is calculated with
85     higher than double precision. For 36.736801<|x|<710.47586
86     exp(-|x|) becomes too small and only exp(|x|) is calculated.
87    */
88
89    const S: f64 = f64::from_bits(0x40b71547652b82fe);
90    let ax = x.abs();
91    let v0 = f_fmla(ax, S, f64::from_bits(0x4198000002000000));
92    let jt = v0.to_bits();
93    let mut v = v0.to_bits();
94    v &= 0xfffffffffc000000;
95    let t = f64::from_bits(v) - f64::from_bits(0x4198000000000000);
96    let ix = ax.to_bits();
97    let aix = ix;
98    if aix < 0x3fc0000000000000u64 {
99        if aix < 0x3e50000000000000u64 {
100            return f_fmla(ax, f64::from_bits(0x3c80000000000000), 1.);
101        }
102        const C: [u64; 5] = [
103            0x3fe0000000000000,
104            0x3fa555555555554e,
105            0x3f56c16c16c26737,
106            0x3efa019ffbbcdbda,
107            0x3e927ffe2df106cb,
108        ];
109        let x2 = x * x;
110        let x4 = x2 * x2;
111
112        let p0 = f_fmla(x2, f64::from_bits(C[3]), f64::from_bits(C[2]));
113        let p1 = f_fmla(x2, f64::from_bits(C[1]), f64::from_bits(C[0]));
114        let p2 = f_fmla(x4, f64::from_bits(C[4]), p0);
115
116        let p = x2 * f_fmla(x4, p2, p1);
117        let e = x2 * (4. * f64::from_bits(0x3ca0000000000000));
118        let lb = 1. + (p - e);
119        let ub = 1. + (p + e);
120        if lb == ub {
121            return lb;
122        }
123        return as_cosh_zero(x);
124    }
125
126    // treat large values apart to avoid a spurious invalid exception
127    if aix > 0x408633ce8fb9f87du64 {
128        // |x| > 0x1.633ce8fb9f87dp+9
129        if aix > 0x7ff0000000000000u64 {
130            return x + x;
131        } // nan
132        if aix == 0x7ff0000000000000u64 {
133            return x.abs();
134        } // inf
135        return f64::from_bits(0x7fe0000000000000) * 2.0;
136    }
137
138    let il: i64 = ((jt.wrapping_shl(14)) >> 40) as i64;
139    let jl: i64 = -il;
140    let i1 = il & 0x3f;
141    let i0 = (il >> 6) & 0x3f;
142    let ie = il >> 12;
143    let j1 = jl & 0x3f;
144    let j0 = (jl >> 6) & 0x3f;
145    let je = jl >> 12;
146    let mut sp = (1022i64.wrapping_add(ie) as u64).wrapping_shl(52);
147    let sm = (1022i64.wrapping_add(je) as u64).wrapping_shl(52);
148    let sn0 = EXP_REDUCE_T0[i0 as usize];
149    let sn1 = EXP_REDUCE_T1[i1 as usize];
150    let t0h = f64::from_bits(sn0.1);
151    let t0l = f64::from_bits(sn0.0);
152    let t1h = f64::from_bits(sn1.1);
153    let t1l = f64::from_bits(sn1.0);
154    let mut th = t0h * t1h;
155    let mut tl = f_fmla(t0h, t1l, t1h * t0l) + dd_fmla(t0h, t1h, -th);
156
157    const L2H: f64 = f64::from_bits(0x3f262e42ff000000);
158    const L2L: f64 = f64::from_bits(0x3d0718432a1b0e26);
159    let dx = f_fmla(L2L, t, f_fmla(-L2H, t, ax));
160    let dx2 = dx * dx;
161    let mx = -dx;
162    const CH: [u64; 4] = [
163        0x3ff0000000000000,
164        0x3fe0000000000000,
165        0x3fc5555555aaaaae,
166        0x3fa55555551c98c0,
167    ];
168
169    let pw0 = f_fmla(dx, f64::from_bits(CH[3]), f64::from_bits(CH[2]));
170    let pw1 = f_fmla(dx, f64::from_bits(CH[1]), f64::from_bits(CH[0]));
171
172    let pp = dx * f_fmla(dx2, pw0, pw1);
173    let (mut rh, mut rl);
174    if aix > 0x4014000000000000u64 {
175        // |x| > 5
176        if aix > 0x40425e4f7b2737fau64 {
177            // |x| >~ 36.736801
178            sp = (1021i64.wrapping_add(ie) as u64).wrapping_shl(52);
179            rh = th;
180            rl = f_fmla(th, pp, tl);
181            let e = 0.11e-18 * th;
182            let lb = rh + (rl - e);
183            let ub = rh + (rl + e);
184            if lb == ub {
185                return (lb * f64::from_bits(sp)) * 2.;
186            }
187
188            let mut tt = hyperbolic_exp_accurate(ax, t, DoubleDouble::new(tl, th));
189            tt = DoubleDouble::from_exact_add(tt.hi, tt.lo);
190            th = tt.hi;
191            tl = tt.lo;
192            th += tl;
193            th *= 2.;
194            th *= f64::from_bits(sp);
195            return th;
196        }
197
198        let q0h = f64::from_bits(EXP_REDUCE_T0[j0 as usize].1);
199        let q1h = f64::from_bits(EXP_REDUCE_T1[j1 as usize].1);
200        let mut qh = q0h * q1h;
201        th *= f64::from_bits(sp);
202        tl *= f64::from_bits(sp);
203        qh *= f64::from_bits(sm);
204
205        let pmw0 = f_fmla(mx, f64::from_bits(CH[3]), f64::from_bits(CH[2]));
206        let pmw1 = f_fmla(mx, f64::from_bits(CH[1]), f64::from_bits(CH[0]));
207
208        let pm = mx * f_fmla(dx2, pmw0, pmw1);
209        let em = f_fmla(qh, pm, qh);
210        rh = th;
211        rl = f_fmla(th, pp, tl + em);
212        let e = 0.09e-18 * rh;
213        let lb = rh + (rl - e);
214        let ub = rh + (rl + e);
215        if lb == ub {
216            return lb;
217        }
218        let tt = hyperbolic_exp_accurate(ax, t, DoubleDouble::new(tl, th));
219        th = tt.hi;
220        tl = tt.lo;
221        if aix > 0x403f666666666666u64 {
222            rh = th + qh;
223            rl = ((th - rh) + qh) + tl;
224        } else {
225            qh = q0h * q1h;
226            let q0l = f64::from_bits(EXP_REDUCE_T0[j0 as usize].0);
227            let q1l = f64::from_bits(EXP_REDUCE_T1[j1 as usize].0);
228            let mut ql = f_fmla(q0h, q1l, q1h * q0l) + dd_fmla(q0h, q1h, -qh);
229            qh *= f64::from_bits(sm);
230            ql *= f64::from_bits(sm);
231            let qq = hyperbolic_exp_accurate(-ax, -t, DoubleDouble::new(ql, qh));
232            qh = qq.hi;
233            ql = qq.lo;
234            rh = th + qh;
235            rl = (((th - rh) + qh) + ql) + tl;
236        }
237    } else {
238        let tq0 = EXP_REDUCE_T0[j0 as usize];
239        let tq1 = EXP_REDUCE_T1[j1 as usize];
240        let q0h = f64::from_bits(tq0.1);
241        let q0l = f64::from_bits(tq0.0);
242        let q1h = f64::from_bits(tq1.1);
243        let q1l = f64::from_bits(tq1.0);
244        let mut qh = q0h * q1h;
245        let mut ql = f_fmla(q0h, q1l, q1h * q0l) + dd_fmla(q0h, q1h, -qh);
246        th *= f64::from_bits(sp);
247        tl *= f64::from_bits(sp);
248        qh *= f64::from_bits(sm);
249        ql *= f64::from_bits(sm);
250
251        let pmw0 = f_fmla(mx, f64::from_bits(CH[3]), f64::from_bits(CH[2]));
252        let pmw1 = f_fmla(mx, f64::from_bits(CH[1]), f64::from_bits(CH[0]));
253
254        let pm = mx * f_fmla(dx2, pmw0, pmw1);
255        let fph = th;
256        let fpl = f_fmla(th, pp, tl);
257        let fmh = qh;
258        let fml = f_fmla(qh, pm, ql);
259
260        rh = fph + fmh;
261        rl = ((fph - rh) + fmh) + fml + fpl;
262        let e = 0.28e-18 * rh;
263        let lb = rh + (rl - e);
264        let ub = rh + (rl + e);
265        if lb == ub {
266            return lb;
267        }
268        let tt = hyperbolic_exp_accurate(ax, t, DoubleDouble::new(tl, th));
269        let qq = hyperbolic_exp_accurate(-ax, -t, DoubleDouble::new(ql, qh));
270        rh = tt.hi + qq.hi;
271        rl = ((tt.hi - rh) + qq.hi) + qq.lo + tt.lo;
272    }
273    let r = DoubleDouble::from_exact_add(rh, rl);
274    rh = r.hi;
275    rl = r.lo;
276    rh += rl;
277    rh
278}
279
280#[cfg(test)]
281mod tests {
282
283    use super::*;
284
285    #[test]
286    fn test_cosh() {
287        assert_eq!(f_cosh(1.), 1.5430806348152437);
288        assert_eq!(f_cosh(1.5454354343), 2.451616191647056);
289        assert_eq!(f_cosh(15.5454354343), 2820115.088877147);
290    }
291}
pxfm/hyperbolic/cosh.rs

pxfm/hyperbolic/
cosh.rs
