core_math_rs/f16/
asinh.rs

1/* Correctly-rounded hyperbolic arc-sine for binary16 value.
2
3Copyright (c) 2025 Paul Zimmermann
4
5This file is ported from the CORE-MATH project
6(https://core-math.gitlabpages.inria.fr/).
7
8Permission is hereby granted, free of charge, to any person obtaining a copy
9of this software and associated documentation files (the "Software"), to deal
10in the Software without restriction, including without limitation the rights
11to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12copies of the Software, and to permit persons to whom the Software is
13furnished to do so, subject to the following conditions:
14
15The above copyright notice and this permission notice shall be included in all
16copies or substantial portions of the Software.
17
18THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24SOFTWARE.
25*/
26
27use hexf::hexf32 as h;
28
29/* The following polynomials were generated using Sollya (cf asinh.sollya).
30P[0] is a degree-7 polynomial with odd-degree coefficients only
31approximating asinh(x) over [0,2^-2], and for 1 <= i < 19,
32P[i] is a degree-6 polynomial approximating asinh(x) over [2^(i-3),2^(i-2)).
33They were afterwards optimized to decrease the number of exceptions. */
34
35const P: [[f32; 7]; 19] = [
36    [
37        h!("0x1.000002p+0"),
38        h!("-0x1.5554d6p-3"),
39        h!("0x1.32c4a4p-4"),
40        h!("-0x1.5336f2p-5"),
41        0.0,
42        0.0,
43        0.0,
44    ], // [0,2^-2]
45    [
46        h!("0x1.707d7ap-17"),
47        h!("0x1.ffe5ap-1"),
48        h!("0x1.86a284p-10"),
49        h!("-0x1.608874p-3"),
50        h!("0x1.1d36dp-7"),
51        h!("0x1.435f08p-4"),
52        h!("-0x1.2861d6p-5"),
53    ], // [2^-2,2^-1)
54    [
55        h!("0x1.2c296cp-10"),
56        h!("0x1.f9d57p-1"),
57        h!("0x1.b8aab8p-5"),
58        h!("-0x1.30b6a4p-2"),
59        h!("0x1.73af8ap-3"),
60        h!("-0x1.a1564ep-5"),
61        h!("0x1.6908p-8"),
62    ], // [2^-1,2^0)
63    [
64        h!("-0x1.171628p-6"),
65        h!("0x1.1382e6p+0"),
66        h!("-0x1.f8ac18p-4"),
67        h!("-0x1.c6795cp-4"),
68        h!("0x1.2d2dfep-4"),
69        h!("-0x1.38f626p-6"),
70        h!("0x1.f9d716p-10"),
71    ], // [2^0,2^1)
72    [
73        h!("-0x1.29739cp-4"),
74        h!("0x1.43decap+0"),
75        h!("-0x1.91761cp-2"),
76        h!("0x1.82261ap-4"),
77        h!("-0x1.eee226p-7"),
78        h!("0x1.746546p-10"),
79        h!("-0x1.f02e6ap-15"),
80    ], // [2^1,2^2)
81    [
82        h!("0x1.6722aep-3"),
83        h!("0x1.d6b2e4p-1"),
84        h!("-0x1.7dbe34p-3"),
85        h!("0x1.c10cbep-6"),
86        h!("-0x1.502cc2p-9"),
87        h!("0x1.1e2eb2p-13"),
88        h!("-0x1.a56434p-19"),
89    ], // [2^2,2^3)
90    [
91        h!("0x1.6f7ecp-1"),
92        h!("0x1.05fe1p-1"),
93        h!("-0x1.c290eap-5"),
94        h!("0x1.12d0bep-8"),
95        h!("-0x1.a5b282p-13"),
96        h!("0x1.6d57cp-18"),
97        h!("-0x1.108528p-24"),
98    ], // [2^3,2^4)
99    [
100        h!("0x1.5ed2bep+0"),
101        h!("0x1.0d3bp-2"),
102        h!("-0x1.d55ab8p-7"),
103        h!("0x1.20902cp-11"),
104        h!("-0x1.bd0e6ap-17"),
105        h!("0x1.82e9d2p-23"),
106        h!("-0x1.214f5cp-30"),
107    ], // [2^4,2^5)
108    [
109        h!("0x1.06cee4p+1"),
110        h!("0x1.0f17c2p-3"),
111        h!("-0x1.da366ep-9"),
112        h!("0x1.24229ap-14"),
113        h!("-0x1.c327d6p-21"),
114        h!("0x1.889372p-28"),
115        h!("-0x1.25be4cp-36"),
116    ], // [2^5,2^6)
117    [
118        h!("0x1.5f39aep+1"),
119        h!("0x1.0f8076p-4"),
120        h!("-0x1.db39e6p-11"),
121        h!("0x1.24d726p-17"),
122        h!("-0x1.c44a96p-25"),
123        h!("0x1.899114p-33"),
124        h!("-0x1.26784ap-42"),
125    ], // [2^6,2^7)
126    [
127        h!("0x1.b7e21p+1"),
128        h!("0x1.0f9492p-5"),
129        h!("-0x1.db65a6p-13"),
130        h!("0x1.24f0f8p-20"),
131        h!("-0x1.c46c0ep-29"),
132        h!("0x1.89a692p-38"),
133        h!("-0x1.2681d4p-48"),
134    ], // [2^7,2^8)
135    [
136        h!("0x1.084d12p+2"),
137        h!("0x1.0f923p-6"),
138        h!("-0x1.db5616p-15"),
139        h!("0x1.24df02p-23"),
140        h!("-0x1.c4422p-33"),
141        h!("0x1.89751cp-43"),
142        h!("-0x1.2653cp-54"),
143    ], // [2^8,2^9)
144    [
145        h!("0x1.34afcep+2"),
146        h!("0x1.0f76f6p-7"),
147        h!("-0x1.daf3c8p-17"),
148        h!("0x1.24818cp-26"),
149        h!("-0x1.c37d7ep-37"),
150        h!("0x1.889b7ep-48"),
151        h!("-0x1.258ca6p-60"),
152    ], // [2^9,2^10)
153    [
154        h!("0x1.610b6ep+2"),
155        h!("0x1.0f7ad8p-8"),
156        h!("-0x1.db022p-19"),
157        h!("0x1.248fa8p-29"),
158        h!("-0x1.c39c7p-41"),
159        h!("0x1.88bf6p-53"),
160        h!("-0x1.25aeeep-66"),
161    ], // [2^10,2^11)
162    [
163        h!("0x1.8d6a98p+2"),
164        h!("0x1.0f6f2ep-9"),
165        h!("-0x1.dad874p-21"),
166        h!("0x1.24689p-32"),
167        h!("-0x1.c34ab6p-45"),
168        h!("0x1.886532p-58"),
169        h!("-0x1.255d5ep-72"),
170    ], // [2^11,2^12)
171    [
172        h!("0x1.b9c986p+2"),
173        h!("0x1.0f651cp-10"),
174        h!("-0x1.dab52cp-23"),
175        h!("0x1.2447e4p-35"),
176        h!("-0x1.c3079ep-49"),
177        h!("0x1.881ca6p-63"),
178        h!("-0x1.251c46p-78"),
179    ], // [2^12,2^13)
180    [
181        h!("0x1.e6269ap+2"),
182        h!("0x1.0f62aep-11"),
183        h!("-0x1.daac6ep-25"),
184        h!("0x1.243f9ep-38"),
185        h!("-0x1.c2f63cp-53"),
186        h!("0x1.88096cp-68"),
187        h!("-0x1.250ac6p-84"),
188    ], // [2^13,2^14)
189    [
190        h!("0x1.09425p+3"),
191        h!("0x1.0f5c6ap-12"),
192        h!("-0x1.da96ep-27"),
193        h!("0x1.242c12p-41"),
194        h!("-0x1.c2ceccp-57"),
195        h!("0x1.87df72p-73"),
196        h!("-0x1.24e5eep-90"),
197    ], // [2^14,2^15)
198    [
199        h!("0x1.1f707ep+3"),
200        h!("0x1.0f5d12p-13"),
201        h!("-0x1.da990ap-29"),
202        h!("0x1.242df2p-44"),
203        h!("-0x1.c2d272p-61"),
204        h!("0x1.87e32ep-78"),
205        h!("-0x1.24e916p-96"),
206    ], // [2^15,2^16)
207];
208
209/// Correctly-rounded hyperbolic arc-sine for binary16 value.
210pub fn cr_asinhf16(x: f16) -> f16 {
211    let v = f32::from(x);
212    let u = v.to_bits();
213    let au = u & 0x7fffffff;
214
215    if au >= 0x7f800000 {
216        // NaN or Inf
217        // asinh(+Inf) = +Inf, otherwise we get qNaN
218        if (u == 0x7f800000) || ((au & 0x7fffff) != 0) {
219            // +Inf or NaN
220            return x + x;
221        }
222        return x; // -Inf
223    }
224
225    let i = ((au >> 23) as i32) - 124; // 2^(i-3) <= x < 2^(i-2)
226    let mut t = v;
227    let tt = t * t;
228    let mut y;
229    let p;
230    if i <= 0 {
231        // |x| < 2^-2
232        /* For |x| <= 0x1.714p-5, asinh(x) rounds to x to nearest,
233        we deal with that case separately, so that for x subnormal
234        and a power of two, we get an underflow. */
235        if au <= 0x3d38a000 {
236            if au == 0 {
237                return x;
238            }
239            return if au == u {
240                v - h!("0x1p-26")
241            } else {
242                v + h!("0x1p-26")
243            } as f16;
244        }
245        p = P[0];
246        let c5 = f32::mul_add(p[3], tt, p[2]);
247        let mut c1 = f32::mul_add(p[1], tt, p[0]);
248        c1 = f32::mul_add(c5, tt * tt, c1);
249        y = t * c1;
250    } else {
251        // |x| >= 2^-2
252        const S: [f32; 2] = [1.0, -1.0];
253        /* we make t positive since the polynomials for |x| >= 2^-2 are not odd,
254        thus only work for x > 0 */
255        t *= S[(u >> 31) as usize];
256        p = P[i as usize];
257        let mut c4 = f32::mul_add(p[5], t, p[4]);
258        let c2 = f32::mul_add(p[3], t, p[2]);
259        let mut c0 = f32::mul_add(p[1], t, p[0]);
260        c4 = f32::mul_add(p[6], tt, c4);
261        c0 = f32::mul_add(c2, tt, c0);
262        y = f32::mul_add(c4, tt * tt, c0);
263
264        // deal with exceptions
265        match (i, au) {
266            // i is the exponent of x, plus 3
267            (7, 0x41936000) => {
268                y = h!("0x1.cdbffep+1"); // |x|=0x1.26cp+4
269            }
270            (10, 0x436cc000) => {
271                y = h!("0x1.8a4002p+2"); // |x|=0x1.d98p+7
272            }
273            (11, 0x43f10000) => {
274                y = h!("0x1.b7bffap+2"); // |x|=0x1.e2p+8
275            }
276            (12, 0x44588000) => {
277                y = h!("0x1.dd4004p+2"); // |x|=0x1.b1p+9
278            }
279            _ => (),
280        }
281
282        y *= S[(u >> 31) as usize]; // restore sign
283    }
284
285    y as f16
286}
287
288#[cfg(test)]
289mod tests {
290    #[test]
291    fn exhaustive() {
292        for b in 0..=u16::MAX {
293            let x = f16::from_bits(b);
294            let y1 = super::cr_asinhf16(x);
295            let y2 = core_math::asinhf16(x);
296            assert_eq!(
297                y1.to_bits(),
298                y2.to_bits(),
299                "asinhf16({x} @ {b:#04x}) = ({y1} @ {y1b:#04x}) vs ({y2} @ {y2b:#04x})",
300                y1b = y1.to_bits(),
301                y2b = y2.to_bits()
302            );
303        }
304    }
305
306    #[test]
307    fn edge() {
308        assert_eq!(
309            super::cr_asinhf16(-f16::NAN).to_bits(),
310            (-f16::NAN).to_bits()
311        );
312        assert_eq!(super::cr_asinhf16(-f16::INFINITY), -f16::INFINITY);
313        assert_eq!(super::cr_asinhf16(-0.0).to_bits(), (-0.0_f16).to_bits());
314        assert_eq!(super::cr_asinhf16(0.0).to_bits(), (0.0_f16).to_bits());
315        assert_eq!(super::cr_asinhf16(f16::INFINITY), f16::INFINITY);
316        assert_eq!(super::cr_asinhf16(f16::NAN).to_bits(), (f16::NAN).to_bits());
317    }
318}
core_math_rs/f16/asinh.rs

core_math_rs/f16/
asinh.rs
