pxfm/hyperbolic/
asinh.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::hyperbolic::acosh::{
32    ACOSH_ASINH_B, ACOSH_ASINH_L1, ACOSH_ASINH_L2, ACOSH_ASINH_LL, ACOSH_ASINH_R1, ACOSH_ASINH_R2,
33    ACOSH_ASINH_REFINE_T2, ACOSH_ASINH_REFINE_T4, ACOSH_SINH_REFINE_T1, ACOSH_SINH_REFINE_T3,
34    lpoly_xd_generic,
35};
36
37#[cold]
38fn asinh_refine(x: f64, a: f64, z: DoubleDouble) -> f64 {
39    let mut t = z.hi.to_bits();
40    let ex: i32 = (t >> 52) as i32;
41    let e = ex.wrapping_sub(0x3ff) + if z.lo == 0.0 { 1i32 } else { 0i32 };
42    t &= 0x000fffffffffffff;
43    t |= 0x3ffu64 << 52;
44    let ed = e as f64;
45    let v = (a - ed + f64::from_bits(0x3ff0000800000000)).to_bits();
46    let i = (v.wrapping_sub(0x3ffu64 << 52)) >> (52 - 16);
47    let i1 = (i >> 12) & 0x1f;
48    let i2 = (i >> 8) & 0xf;
49    let i3 = (i >> 4) & 0xf;
50    let i4 = i & 0xf;
51    const L20: f64 = f64::from_bits(0x3fd62e42fefa3800);
52    const L21: f64 = f64::from_bits(0x3d1ef35793c76800);
53    const L22: f64 = f64::from_bits(0xba49ff0342542fc3);
54    let el2 = L22 * ed;
55    let el1 = L21 * ed;
56    let el0 = L20 * ed;
57    let mut dl0: f64;
58
59    let ll0i1 = ACOSH_ASINH_LL[0][i1 as usize];
60    let ll1i2 = ACOSH_ASINH_LL[1][i2 as usize];
61    let ll2i3 = ACOSH_ASINH_LL[2][i3 as usize];
62    let ll3i4 = ACOSH_ASINH_LL[3][i4 as usize];
63
64    dl0 = f64::from_bits(ll0i1.0)
65        + f64::from_bits(ll1i2.0)
66        + (f64::from_bits(ll2i3.0) + f64::from_bits(ll3i4.0));
67    let dl1 = f64::from_bits(ll0i1.1)
68        + f64::from_bits(ll1i2.1)
69        + (f64::from_bits(ll2i3.1) + f64::from_bits(ll3i4.1));
70    let dl2 = f64::from_bits(ll0i1.2)
71        + f64::from_bits(ll1i2.2)
72        + (f64::from_bits(ll2i3.2) + f64::from_bits(ll3i4.2));
73    dl0 += el0;
74    let t12 = f64::from_bits(ACOSH_SINH_REFINE_T1[i1 as usize])
75        * f64::from_bits(ACOSH_ASINH_REFINE_T2[i2 as usize]);
76    let t34 = f64::from_bits(ACOSH_SINH_REFINE_T3[i3 as usize])
77        * f64::from_bits(ACOSH_ASINH_REFINE_T4[i4 as usize]);
78    let th = t12 * t34;
79    let tl = dd_fmla(t12, t34, -th);
80    let dh = th * f64::from_bits(t);
81    let dl = dd_fmla(th, f64::from_bits(t), -dh);
82    let sh = tl * f64::from_bits(t);
83    let sl = dd_fmla(tl, f64::from_bits(t), -sh);
84
85    let mut dx = DoubleDouble::from_exact_add(dh - 1., dl);
86    if z.lo != 0.0 {
87        t = z.lo.to_bits();
88        t = t.wrapping_sub((e as i64).wrapping_shl(52) as u64);
89        dx.lo += th * f64::from_bits(t);
90    }
91    dx = DoubleDouble::add(dx, DoubleDouble::new(sl, sh));
92    const CL: [u64; 3] = [0xbfc0000000000000, 0x3fb9999999a0754f, 0xbfb55555555c3157];
93
94    let sl0 = f_fmla(dx.hi, f64::from_bits(CL[2]), f64::from_bits(CL[1]));
95    let sl1 = f_fmla(dx.hi, sl0, f64::from_bits(CL[0]));
96
97    let sl = dx.hi * sl1;
98    const CH: [(u64, u64); 3] = [
99        (0x39024b67ee516e3b, 0x3fe0000000000000),
100        (0xb91932ce43199a8d, 0xbfd0000000000000),
101        (0x3c655540c15cf91f, 0x3fc5555555555555),
102    ];
103    let mut s = lpoly_xd_generic(dx, CH, sl);
104    s = DoubleDouble::quick_mult(dx, s);
105    s = DoubleDouble::add(s, DoubleDouble::new(el2, el1));
106    s = DoubleDouble::add(s, DoubleDouble::new(dl2, dl1));
107    let mut v02 = DoubleDouble::from_exact_add(dl0, s.hi);
108    let mut v12 = DoubleDouble::from_exact_add(v02.lo, s.lo);
109    let scale = f64::copysign(2., x);
110    v02.hi *= scale;
111    v12.hi *= scale;
112    v12.lo *= scale;
113    t = v12.hi.to_bits();
114    if (t & 0x000fffffffffffff) == 0 {
115        let w = v12.lo.to_bits();
116        if ((w ^ t) >> 63) != 0 {
117            t = t.wrapping_sub(1);
118        } else {
119            t = t.wrapping_add(1);
120        }
121        v12.hi = f64::from_bits(t);
122    }
123    v02.hi + v12.hi
124}
125
126#[cold]
127fn as_asinh_zero(x: f64, x2h: f64, x2l: f64) -> f64 {
128    static CH: [(u64, u64); 12] = [
129        (0xbc65555555555555, 0xbfc5555555555555),
130        (0x3c499999999949df, 0x3fb3333333333333),
131        (0x3c32492496091b0c, 0xbfa6db6db6db6db7),
132        (0x3c1c71a35cfa0671, 0x3f9f1c71c71c71c7),
133        (0x3c317f937248cf81, 0xbf96e8ba2e8ba2e9),
134        (0xbc374e3c1dfd4c3d, 0x3f91c4ec4ec4ec4f),
135        (0xbc238e7a467ecc55, 0xbf8c999999999977),
136        (0x3c2a83c7bace55eb, 0x3f87a87878786c7e),
137        (0xbc2d024df7fa0542, 0xbf83fde50d764083),
138        (0xbc2ba9c13deb261f, 0x3f812ef3ceae4d12),
139        (0xbc1546da9bc5b32a, 0xbf7df3bd104aa267),
140        (0x3c140d284a1d67f9, 0x3f7a685fc5de7a04),
141    ];
142
143    const CL: [u64; 5] = [
144        0xbf77828d553ec800,
145        0x3f751712f7bee368,
146        0xbf72e6d98527bcc6,
147        0x3f70095da47b392c,
148        0xbf63b92d6368192c,
149    ];
150
151    let yw0 = f_fmla(x2h, f64::from_bits(CL[4]), f64::from_bits(CL[3]));
152    let yw1 = f_fmla(x2h, yw0, f64::from_bits(CL[2]));
153    let yw2 = f_fmla(x2h, yw1, f64::from_bits(CL[1]));
154
155    let y2 = x2h * f_fmla(x2h, yw2, f64::from_bits(CL[0]));
156    let mut y1 = lpoly_xd_generic(DoubleDouble::new(x2l, x2h), CH, y2);
157    y1 = DoubleDouble::quick_mult(y1, DoubleDouble::new(x2l, x2h));
158    y1 = DoubleDouble::quick_mult_f64(y1, x);
159
160    let y0 = DoubleDouble::from_exact_add(x, y1.hi);
161    let mut p = DoubleDouble::from_exact_add(y0.lo, y1.lo);
162
163    let mut t = p.hi.to_bits();
164    if (t & 0x000fffffffffffff) == 0 {
165        let w = p.lo.to_bits();
166        if ((w ^ t) >> 63) != 0 {
167            t = t.wrapping_sub(1);
168        } else {
169            t = t.wrapping_add(1);
170        }
171        p.hi = f64::from_bits(t);
172    }
173    y0.hi + p.hi
174}
175
176/// Huperbolic sine function
177///
178/// Max ULP 0.5
179pub fn f_asinh(x: f64) -> f64 {
180    let ax = x.abs();
181    let ix = ax.to_bits();
182    let u = ix;
183    if u < 0x3fbb000000000000u64 {
184        // |x| < 0x1.bp-4
185        // for |x| < 0x1.7137449123ef7p-26, asinh(x) rounds to x to nearest
186        // for |x| < 0x1p-1022 we have underflow but not for 0x1p-1022 (to nearest)
187        if u < 0x3e57137449123ef7u64 {
188            // |x| < 0x1.7137449123ef7p-26
189            if u == 0 {
190                return x;
191            }
192            let res = f_fmla(f64::from_bits(0xbc30000000000000), x, x);
193            return res;
194        }
195        let x2h = x * x;
196        let x2l = dd_fmla(x, x, -x2h);
197        let x3h = x2h * x;
198        let sl;
199        if u < 0x3f93000000000000u64 {
200            // |x| < 0x1.3p-6
201            if u < 0x3f30000000000000u64 {
202                // |x| < 0x1p-12
203                if u < 0x3e5a000000000000u64 {
204                    // |x| < 0x1.ap-26
205                    sl = x3h * f64::from_bits(0xbfc5555555555555);
206                } else {
207                    const CL: [u64; 2] = [0xbfc5555555555555, 0x3fb3333327c57c60];
208                    let sl0 = f_fmla(x2h, f64::from_bits(CL[1]), f64::from_bits(CL[0]));
209                    sl = x3h * sl0;
210                }
211            } else {
212                const CL: [u64; 4] = [
213                    0xbfc5555555555555,
214                    0x3fb333333332f2ff,
215                    0xbfa6db6d9a665159,
216                    0x3f9f186866d775f0,
217                ];
218                let sl0 = f_fmla(x2h, f64::from_bits(CL[3]), f64::from_bits(CL[2]));
219                let sl1 = f_fmla(x2h, sl0, f64::from_bits(CL[1]));
220                sl = x3h * f_fmla(x2h, sl1, f64::from_bits(CL[0]));
221            }
222        } else {
223            const CL: [u64; 7] = [
224                0xbfc5555555555555,
225                0x3fb3333333333310,
226                0xbfa6db6db6da466c,
227                0x3f9f1c71c2ea7be4,
228                0xbf96e8b651b09d72,
229                0x3f91c309fc0e69c2,
230                0xbf8bab7833c1e000,
231            ];
232            let c1 = f_fmla(x2h, f64::from_bits(CL[2]), f64::from_bits(CL[1]));
233            let c3 = f_fmla(x2h, f64::from_bits(CL[4]), f64::from_bits(CL[3]));
234            let c5 = f_fmla(x2h, f64::from_bits(CL[6]), f64::from_bits(CL[5]));
235            let x4 = x2h * x2h;
236
237            let sl0 = f_fmla(x4, c5, c3);
238            let sl1 = f_fmla(x4, sl0, c1);
239
240            sl = x3h * f_fmla(x2h, sl1, f64::from_bits(CL[0]));
241        }
242        let eps = f64::from_bits(0x3ca6000000000000) * x3h;
243        let lb = x + (sl - eps);
244        let ub = x + (sl + eps);
245        if lb == ub {
246            return lb;
247        }
248        return as_asinh_zero(x, x2h, x2l);
249    }
250
251    // |x| >= 0x1.bp-4
252    let mut x2h: f64 = 0.;
253    let mut x2l: f64 = 0.;
254    let mut off: i32 = 0x3ff;
255
256    let va: DoubleDouble = if u < 0x4190000000000000 {
257        // x < 0x1p+26
258        x2h = x * x;
259        x2l = dd_fmla(x, x, -x2h);
260        let mut dt = if u < 0x3ff0000000000000 {
261            DoubleDouble::from_exact_add(1., x2h)
262        } else {
263            DoubleDouble::from_exact_add(x2h, 1.)
264        };
265        dt.lo += x2l;
266
267        let ah = dt.hi.sqrt();
268        let rs = 0.5 / dt.hi;
269        let al = (dt.lo - dd_fmla(ah, ah, -dt.hi)) * (rs * ah);
270        let mut ma = DoubleDouble::from_exact_add(ah, ax);
271        ma.lo += al;
272        ma
273    } else if u < 0x4330000000000000 {
274        DoubleDouble::new(0.5 / ax, 2. * ax)
275    } else {
276        if u >= 0x7ff0000000000000u64 {
277            return x + x;
278        } // +-inf or nan
279        off = 0x3fe;
280        DoubleDouble::new(0., ax)
281    };
282
283    let mut t = va.hi.to_bits();
284    let ex = (t >> 52) as i32;
285    let e = ex.wrapping_sub(off);
286    t &= 0x000fffffffffffff;
287    let ed = e as f64;
288    let i = t >> (52 - 5);
289    let d = (t & 0x00007fffffffffff) as i64;
290    let b_i = ACOSH_ASINH_B[i as usize];
291    let j: u64 = t
292        .wrapping_add((b_i[0] as u64).wrapping_shl(33))
293        .wrapping_add((b_i[1] as i64).wrapping_mul(d >> 16) as u64)
294        >> (52 - 10);
295    t |= 0x3ffu64 << 52;
296    let i1 = (j >> 5) as i32;
297    let i2 = j & 0x1f;
298    let r =
299        f64::from_bits(ACOSH_ASINH_R1[i1 as usize]) * f64::from_bits(ACOSH_ASINH_R2[i2 as usize]);
300    let dx = dd_fmla(r, f64::from_bits(t), -1.);
301    let dx2 = dx * dx;
302    const C: [u64; 5] = [
303        0xbfe0000000000000,
304        0x3fd5555555555530,
305        0xbfcfffffffffffa0,
306        0x3fc99999e33a6366,
307        0xbfc555559ef9525f,
308    ];
309
310    let fw0 = f_fmla(dx, f64::from_bits(C[3]), f64::from_bits(C[2]));
311    let fw1 = f_fmla(dx, f64::from_bits(C[1]), f64::from_bits(C[0]));
312    let fw2 = f_fmla(dx2, f64::from_bits(C[4]), fw0);
313
314    let f = dx2 * f_fmla(dx2, fw2, fw1);
315    const L2H: f64 = f64::from_bits(0x3fe62e42fefa3800);
316    const L2L: f64 = f64::from_bits(0x3d2ef35793c76730);
317    let l1i1 = ACOSH_ASINH_L1[i1 as usize];
318    let l1i2 = ACOSH_ASINH_L2[i2 as usize];
319    let mut lh = f_fmla(L2H, ed, f64::from_bits(l1i1.1) + f64::from_bits(l1i2.1));
320    let mut ll = f_fmla(
321        L2L,
322        ed,
323        f64::from_bits(l1i1.0) + f64::from_bits(l1i2.0) + va.lo / va.hi + f,
324    );
325    ll += dx;
326    lh *= f64::copysign(1., x);
327    ll *= f64::copysign(1., x);
328    let eps = 1.63e-19;
329    let lb = lh + (ll - eps);
330    let ub = lh + (ll + eps);
331    if lb == ub {
332        return lb;
333    }
334    if ax < f64::from_bits(0x3fd0000000000000) {
335        return as_asinh_zero(x, x2h, x2l);
336    }
337    asinh_refine(x, f64::from_bits(0x3ff71547652b82fe) * lb.abs(), va)
338}
339
340#[cfg(test)]
341mod tests {
342    use super::*;
343
344    #[test]
345    fn test_asinh() {
346        assert_eq!(f_asinh(-0.05268859863273256), -0.05266425100170862);
347        assert_eq!(f_asinh(1.05268859863273256), 0.9181436936151385);
348    }
349}
pxfm/hyperbolic/asinh.rs

pxfm/hyperbolic/
asinh.rs
