pxfm/hyperbolic/
sinh.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;
33
34#[cold]
35pub(crate) fn hyperbolic_exp_accurate(x: f64, t: f64, zt: DoubleDouble) -> DoubleDouble {
36    static CH: [(u64, u64); 3] = [
37        (0x3a16c16bd194535d, 0x3ff0000000000000),
38        (0xba28259d904fd34f, 0x3fe0000000000000),
39        (0x3c653e93e9f26e62, 0x3fc5555555555555),
40    ];
41    const L2H: f64 = f64::from_bits(0x3f262e42ff000000);
42    const L2L: f64 = f64::from_bits(0x3d0718432a1b0e26);
43    const L2LL: f64 = f64::from_bits(0x3999ff0342542fc3);
44    let dx = x - L2H * t;
45    let mut dxl = L2L * t;
46    let dxll = f_fmla(L2LL, t, dd_fmla(L2L, t, -dxl));
47    let dxh = dx + dxl;
48    dxl = ((dx - dxh) + dxl) + dxll;
49
50    let fl0 = f_fmla(
51        dxh,
52        f64::from_bits(0x3f56c16c169400a7),
53        f64::from_bits(0x3f811111113e93e9),
54    );
55
56    let fl = dxh * f_fmla(dxh, fl0, f64::from_bits(0x3fa5555555555555));
57    let mut f = lpoly_xd_generic(DoubleDouble::new(dxl, dxh), CH, fl);
58    f = DoubleDouble::quick_mult(DoubleDouble::new(dxl, dxh), f);
59    f = DoubleDouble::quick_mult(zt, f);
60    let zh = zt.hi + f.hi;
61    let zl = (zt.hi - zh) + f.hi;
62    let uh = zh + zt.lo;
63    let ul = ((zh - uh) + zt.lo) + zl;
64    let vh = uh + f.lo;
65    let vl = ((uh - vh) + f.lo) + ul;
66    DoubleDouble::new(vl, vh)
67}
68
69#[cold]
70fn as_sinh_zero(x: f64) -> f64 {
71    static CH: [(u64, u64); 5] = [
72        (0x3c6555555555552f, 0x3fc5555555555555),
73        (0x3c011111115cf00d, 0x3f81111111111111),
74        (0x3b6a0011c925b85c, 0x3f2a01a01a01a01a),
75        (0xbb6b4e2835532bcd, 0x3ec71de3a556c734),
76        (0xbaedefcf17a6ab79, 0x3e5ae64567f54482),
77    ];
78    let d2x = DoubleDouble::from_exact_mult(x, x);
79
80    let yw0 = f_fmla(
81        d2x.hi,
82        f64::from_bits(0x3ce95785063cd974),
83        f64::from_bits(0x3d6ae7f36beea815),
84    );
85
86    let y2 = d2x.hi * f_fmla(d2x.hi, yw0, f64::from_bits(0x3de6124613aef206));
87    let mut y1 = lpoly_xd_generic(d2x, CH, y2);
88    y1 = DoubleDouble::quick_mult_f64(y1, x);
89    y1 = DoubleDouble::quick_mult(y1, d2x); // y2 = y1.l
90    let y0 = DoubleDouble::from_exact_add(x, y1.hi); // y0 = y0.hi
91    let mut p = DoubleDouble::from_exact_add(y0.lo, y1.lo);
92    let mut t = p.hi.to_bits();
93    if (t & 0x000fffffffffffff) == 0 {
94        let w = p.lo.to_bits();
95        if ((w ^ t) >> 63) != 0 {
96            t = t.wrapping_sub(1);
97        } else {
98            t = t.wrapping_add(1);
99        }
100        p.hi = f64::from_bits(t);
101    }
102    y0.hi + p.hi
103}
104
105/// Hyperbolic sine function
106///
107/// Max ULP 0.5
108pub fn f_sinh(x: f64) -> f64 {
109    /*
110     The function sinh(x) is approximated by a minimax polynomial for
111     |x|<0.25. For other arguments the identity
112     sinh(x)=(exp(|x|)-exp(-|x|))/2*copysign(1,x) is used. For |x|<5
113     both exponents are calculated with slightly higher precision than
114     double. For 5<|x|<36.736801 the exp(-|x|) is small and is
115     calculated with double precision but exp(|x|) is calculated with
116     higher than double precision. For 36.736801<|x|<710.47586
117     exp(-|x|) becomes too small and only exp(|x|) is calculated.
118    */
119
120    const S: f64 = f64::from_bits(0x40b71547652b82fe);
121    let ax = x.abs();
122    let v0 = dd_fmla(ax, S, f64::from_bits(0x4198000002000000));
123    let jt = v0.to_bits();
124    let v = v0.to_bits() & 0xfffffffffc000000;
125    let t = f64::from_bits(v) - f64::from_bits(0x4198000000000000);
126    let ix = ax.to_bits();
127    let aix = ix;
128    if aix < 0x3fd0000000000000u64 {
129        // |x| < 0x1p-2
130        if aix < 0x3e57137449123ef7u64 {
131            // |x| < 0x1.7137449123ef7p-26
132            /* We have underflow exactly when 0 < |x| < 2^-1022:
133            for RNDU, sinh(2^-1022-2^-1074) would round to 2^-1022-2^-1075
134            with unbounded exponent range */
135            return dd_fmla(x, f64::from_bits(0x3c80000000000000), x);
136        }
137        const C: [u64; 5] = [
138            0x3fc5555555555555,
139            0x3f81111111111087,
140            0x3f2a01a01a12e1c3,
141            0x3ec71de2e415aa36,
142            0x3e5aed2bff4269e6,
143        ];
144        let x2 = x * x;
145        let x3 = x2 * x;
146        let x4 = x2 * x2;
147
148        let pw0 = f_fmla(x2, f64::from_bits(C[3]), f64::from_bits(C[2]));
149        let pw1 = f_fmla(x2, f64::from_bits(C[1]), f64::from_bits(C[0]));
150        let pw2 = f_fmla(x4, f64::from_bits(C[4]), pw0);
151
152        let p = x3 * f_fmla(x4, pw2, pw1);
153        let e = x3 * f64::from_bits(0x3ca9000000000000);
154        let lb = x + (p - e);
155        let ub = x + (p + e);
156        if lb == ub {
157            return lb;
158        }
159        return as_sinh_zero(x);
160    }
161
162    if aix > 0x408633ce8fb9f87du64 {
163        // |x| >~ 710.47586
164        if aix >= 0x7ff0000000000000u64 {
165            return x + x;
166        } // nan Inf
167        return f64::copysign(f64::from_bits(0x7fe0000000000000), x) * 2.0;
168    }
169    let il: i64 = ((jt.wrapping_shl(14)) >> 40) as i64;
170    let jl = -il;
171    let i1 = il & 0x3f;
172    let i0 = (il >> 6) & 0x3f;
173    let ie = il >> 12;
174    let j1 = jl & 0x3f;
175    let j0 = (jl >> 6) & 0x3f;
176    let je = jl >> 12;
177    let mut sp = (1022i64.wrapping_add(ie) as u64).wrapping_shl(52);
178    let sm = (1022i64.wrapping_add(je) as u64).wrapping_shl(52);
179
180    let sn0 = EXP_REDUCE_T0[i0 as usize];
181    let sn1 = EXP_REDUCE_T1[i1 as usize];
182    let t0h = f64::from_bits(sn0.1);
183    let t0l = f64::from_bits(sn0.0);
184    let t1h = f64::from_bits(sn1.1);
185    let t1l = f64::from_bits(sn1.0);
186    let mut th = t0h * t1h;
187    let mut tl = f_fmla(t0h, t1l, t1h * t0l) + dd_fmla(t0h, t1h, -th);
188    const L2H: f64 = f64::from_bits(0x3f262e42ff000000);
189    const L2L: f64 = f64::from_bits(0x3d0718432a1b0e26);
190    let dx = f_fmla(L2L, t, f_fmla(-L2H, t, ax));
191    let dx2 = dx * dx;
192    let mx = -dx;
193    const CH: [u64; 4] = [
194        0x3ff0000000000000,
195        0x3fe0000000000000,
196        0x3fc5555555aaaaae,
197        0x3fa55555551c98c0,
198    ];
199
200    let (mut rl, mut rh);
201
202    let pp0 = f_fmla(dx, f64::from_bits(CH[3]), f64::from_bits(CH[2]));
203    let pp1 = f_fmla(dx, f64::from_bits(CH[1]), f64::from_bits(CH[0]));
204
205    let pp = dx * f_fmla(dx2, pp0, pp1);
206    if aix > 0x4014000000000000u64 {
207        // |x| > 5
208        if aix > 0x40425e4f7b2737fau64 {
209            // |x| >~ 36.736801
210            sp = (1021i64.wrapping_add(ie) as u64).wrapping_shl(52);
211            let mut rh = th;
212            let mut rl = tl + th * pp;
213            rh *= f64::copysign(1., x);
214            rl *= f64::copysign(1., x);
215            let e = 0.11e-18 * th;
216            let lb = rh + (rl - e);
217            let ub = rh + (rl + e);
218            if lb == ub {
219                return (lb * f64::from_bits(sp)) * 2.;
220            }
221            let mut tt = hyperbolic_exp_accurate(ax, t, DoubleDouble::new(tl, th));
222            tt = DoubleDouble::from_exact_add(tt.hi, tt.lo);
223            th = tt.hi;
224            tl = tt.lo;
225            th *= f64::copysign(1., x);
226            tl *= f64::copysign(1., x);
227            th += tl;
228            th *= 2.;
229            th *= f64::from_bits(sp);
230            return th;
231        }
232
233        let q0h = f64::from_bits(EXP_REDUCE_T0[j0 as usize].1);
234        let q1h = f64::from_bits(EXP_REDUCE_T1[j1 as usize].1);
235        let mut qh = q0h * q1h;
236        th *= f64::from_bits(sp);
237        tl *= f64::from_bits(sp);
238        qh *= f64::from_bits(sm);
239
240        let pm0 = f_fmla(mx, f64::from_bits(CH[3]), f64::from_bits(CH[2]));
241        let pm1 = f_fmla(mx, f64::from_bits(CH[1]), f64::from_bits(CH[0]));
242
243        let pm = mx * f_fmla(dx2, pm0, pm1);
244        let em = f_fmla(qh, pm, qh);
245        rh = th;
246        rl = f_fmla(th, pp, tl - em);
247
248        rh *= f64::copysign(1., x);
249        rl *= f64::copysign(1., x);
250        let e = 0.09e-18 * rh;
251        let lb = rh + (rl - e);
252        let ub = rh + (rl + e);
253        if lb == ub {
254            return lb;
255        }
256
257        let tt = hyperbolic_exp_accurate(ax, t, DoubleDouble::new(tl, th));
258        th = tt.hi;
259        tl = tt.lo;
260        if aix > 0x403f666666666666u64 {
261            rh = th - qh;
262            rl = ((th - rh) - qh) + tl;
263        } else {
264            qh = q0h * q1h;
265            let q0l = f64::from_bits(EXP_REDUCE_T0[j0 as usize].0);
266            let q1l = f64::from_bits(EXP_REDUCE_T1[j1 as usize].0);
267            let mut ql = f_fmla(q0h, q1l, q1h * q0l) + dd_fmla(q0h, q1h, -qh);
268            qh *= f64::from_bits(sm);
269            ql *= f64::from_bits(sm);
270            let qq = hyperbolic_exp_accurate(-ax, -t, DoubleDouble::new(ql, qh));
271            rh = th - qq.hi;
272            rl = (((th - rh) - qq.hi) - qq.lo) + tl;
273        }
274    } else {
275        let tq0 = EXP_REDUCE_T0[j0 as usize];
276        let tq1 = EXP_REDUCE_T1[j1 as usize];
277        let q0h = f64::from_bits(tq0.1);
278        let q0l = f64::from_bits(tq0.0);
279        let q1h = f64::from_bits(tq1.1);
280        let q1l = f64::from_bits(tq1.0);
281        let mut qh = q0h * q1h;
282        let mut ql = f_fmla(q0h, q1l, q1h * q0l) + dd_fmla(q0h, q1h, -qh);
283        th *= f64::from_bits(sp);
284        tl *= f64::from_bits(sp);
285        qh *= f64::from_bits(sm);
286        ql *= f64::from_bits(sm);
287
288        let pm0 = f_fmla(mx, f64::from_bits(CH[3]), f64::from_bits(CH[2]));
289        let pm1 = f_fmla(mx, f64::from_bits(CH[1]), f64::from_bits(CH[0]));
290
291        let pm = mx * f_fmla(dx2, pm0, pm1);
292        let fph = th;
293        let fpl = f_fmla(th, pp, tl);
294        let fmh = qh;
295        let fml = f_fmla(qh, pm, ql);
296
297        rh = fph - fmh;
298        rl = ((fph - rh) - fmh) - fml + fpl;
299        rh *= f64::copysign(1., x);
300        rl *= f64::copysign(1., x);
301        let e = 0.28e-18 * rh;
302        let lb = rh + (rl - e);
303        let ub = rh + (rl + e);
304        if lb == ub {
305            return lb;
306        }
307        let tt = hyperbolic_exp_accurate(ax, t, DoubleDouble::new(tl, th));
308        let qq = hyperbolic_exp_accurate(-ax, -t, DoubleDouble::new(ql, qh));
309        rh = tt.hi - qq.hi;
310        rl = ((tt.hi - rh) - qq.hi) - qq.lo + tt.lo;
311    }
312    let r = DoubleDouble::from_exact_add(rh, rl);
313    rh = r.hi;
314    rl = r.lo;
315    rh *= f64::copysign(1., x);
316    rl *= f64::copysign(1., x);
317    rh += rl;
318    rh
319}
320
321#[cfg(test)]
322mod tests {
323    use super::*;
324
325    #[test]
326    fn test_f_sinh() {
327        assert_eq!(f_sinh(1.), 1.1752011936438014);
328    }
329}
pxfm/hyperbolic/sinh.rs

pxfm/hyperbolic/
sinh.rs
