core_math_rs/f16/
acosh.rs

1/* Correctly-rounded hyperbolic arc-cosine 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
29use super::utils::snanf16;
30
31/* The following polynomials were generated using Sollya (cf acosh.sollya).
32For 0 <= i < 16, P[i] is a degree-5 polynomial approximating
33acosh(x)/sqrt(x-1) for 0 <= i < 5, and acosh(x) for 5 <= i < 16.
34They were afterwards optimized to decrease the number of exceptions. */
35
36// ([+-]?0x[0-9a-f]+(\.[0-9a-f]+)?p[+-][0-9]+) -> h!("$1")
37const P: [[f32; 6]; 16] = [
38    [
39        h!("0x1.91940ap+0"),
40        h!("-0x1.a1ee52p-3"),
41        h!("0x1.0dc64ap-4"),
42        h!("-0x1.453108p-6"),
43        h!("0x1.040a46p-8"),
44        h!("-0x1.803198p-12"),
45    ], /* [2^0,2^1) */
46    [
47        h!("0x1.8eb5fap+0"),
48        h!("-0x1.6c52a2p-3"),
49        h!("0x1.4dddcp-5"),
50        h!("-0x1.e9e95cp-8"),
51        h!("0x1.b697b2p-11"),
52        h!("-0x1.5b2df4p-15"),
53    ], /* [2^1,2^2) */
54    [
55        h!("0x1.854f5ap+0"),
56        h!("-0x1.11e128p-3"),
57        h!("0x1.362606p-6"),
58        h!("-0x1.fb23bp-10"),
59        h!("0x1.e2d03p-14"),
60        h!("-0x1.8d21eap-19"),
61    ], /* [2^2,2^3) */
62    [
63        h!("0x1.717322p+0"),
64        h!("-0x1.608d74p-4"),
65        h!("0x1.c3ceb2p-8"),
66        h!("-0x1.8859d4p-12"),
67        h!("0x1.828e08p-17"),
68        h!("-0x1.44eab2p-23"),
69    ], /* [2^3,2^4) */
70    [
71        h!("0x1.52b1b4p+0"),
72        h!("-0x1.8e60acp-5"),
73        h!("0x1.12e502p-9"),
74        h!("-0x1.ee8866p-15"),
75        h!("0x1.f10d92p-21"),
76        h!("-0x1.a71f3ep-28"),
77    ], /* [2^4,2^5) */
78    [
79        h!("0x1.1b4176p+1"),
80        h!("0x1.c58fdp-4"),
81        h!("-0x1.3d491cp-9"),
82        h!("0x1.244e58p-15"),
83        h!("-0x1.2b59e6p-22"),
84        h!("0x1.029682p-30"),
85    ], /* [2^5,2^6) */
86    [
87        h!("0x1.744536p+1"),
88        h!("0x1.c4d402p-5"),
89        h!("-0x1.3c787ap-11"),
90        h!("0x1.235b58p-18"),
91        h!("-0x1.2a375cp-26"),
92        h!("0x1.017ceap-35"),
93    ], /* [2^6,2^7) */
94    // 4 exceptions
95    [
96        h!("0x1.cd16eap+1"),
97        h!("0x1.c490cp-6"),
98        h!("-0x1.3c27dep-13"),
99        h!("0x1.22f70cp-21"),
100        h!("-0x1.29b90ep-30"),
101        h!("0x1.00fd26p-40"),
102    ], /* [2^7,2^8) */
103    // 1 exception
104    [
105        h!("0x1.12ec74p+2"),
106        h!("0x1.c475e4p-7"),
107        h!("-0x1.3c05dap-15"),
108        h!("0x1.22cb2cp-24"),
109        h!("-0x1.29804ap-34"),
110        h!("0x1.00c29ap-45"),
111    ], /* [2^8,2^9) */
112    // 6 exceptions
113    [
114        h!("0x1.3f4c7ep+2"),
115        h!("0x1.c45ddap-8"),
116        h!("-0x1.3be5p-17"),
117        h!("0x1.229e72p-27"),
118        h!("-0x1.2943fp-38"),
119        h!("0x1.008224p-50"),
120    ], /* [2^9,2^10) */
121    // // 1 exception
122    [
123        h!("0x1.6baa5p+2"),
124        h!("0x1.c45576p-9"),
125        h!("-0x1.3bd988p-19"),
126        h!("0x1.228f6ep-30"),
127        h!("-0x1.29308ep-42"),
128        h!("0x1.006ddp-55"),
129    ], /* [2^10,2^11) */
130    // no exception
131    [
132        h!("0x1.980882p+2"),
133        h!("0x1.c449a8p-10"),
134        h!("-0x1.3bc944p-21"),
135        h!("0x1.2278ep-33"),
136        h!("-0x1.2911bap-46"),
137        h!("0x1.004d3p-60"),
138    ], /* [2^11,2^12) */
139    // no exception
140    [
141        h!("0x1.c46664p+2"),
142        h!("0x1.c44006p-11"),
143        h!("-0x1.3bbbdp-23"),
144        h!("0x1.22665ep-36"),
145        h!("-0x1.28f8a4p-50"),
146        h!("0x1.00326p-65"),
147    ], /* [2^12,2^13) */
148    // 1 exception
149    [
150        h!("0x1.f0c4d6p+2"),
151        h!("0x1.c43246p-12"),
152        h!("-0x1.3ba86ap-25"),
153        h!("0x1.224b62p-39"),
154        h!("-0x1.28d39ep-54"),
155        h!("0x1.000a46p-70"),
156    ], /* [2^13,2^14) */
157    // 1 exception
158    [
159        h!("0x1.0e912cp+3"),
160        h!("0x1.c42b4ep-13"),
161        h!("-0x1.3b9eb4p-27"),
162        h!("0x1.223e1p-42"),
163        h!("-0x1.28c1ap-58"),
164        h!("0x1.ffee3cp-76"),
165    ], /* [2^14,2^15) */
166    // no exception
167    [
168        h!("0x1.24bfccp+3"),
169        h!("0x1.c4263p-14"),
170        h!("-0x1.3b97b8p-29"),
171        h!("0x1.22347ap-45"),
172        h!("-0x1.28b4c4p-62"),
173        h!("0x1.ffd2ecp-81"),
174    ], /* [2^15,2^16) */ // no exception
175];
176
177/// Correctly-rounded hyperbolic arc-cosine for binary16 value.
178pub fn cr_acoshf16(x: f16) -> f16 {
179    let v = f32::from(x);
180    let u = v.to_bits();
181    let au = u & 0x7fffffff;
182
183    if au >= 0x7f800000 {
184        // NaN or Inf
185        // acosh(+Inf) = +Inf, otherwise we get qNaN
186        if (u == 0x7f800000) || ((au & 0x7fffff) != 0) {
187            // +Inf or NaN
188            return x + x;
189        }
190        return snanf16(); // -Inf
191    }
192
193    if (u >= 0x80000000) || (au <= 0x3f800000) {
194        // x <= 1
195        if u == 0x3f800000 {
196            // x = 1
197            return 0.0; // for x=1 the code below yields -0
198        }
199        // for x=1 the code below yields -0
200        return snanf16();
201    }
202
203    let i = (au >> 23) - 127; // 2^i <= x < 2^(i+1)
204    let t = v;
205    let tt = t * t;
206    let p = P[i as usize];
207    let c4 = f32::mul_add(p[5], t, p[4]);
208    let mut c2 = f32::mul_add(p[3], t, p[2]);
209    c2 = f32::mul_add(c4, tt, c2);
210    let c0 = f32::mul_add(p[1], t, p[0]);
211    let mut y = f32::mul_add(c2, tt, c0);
212    if i < 5 {
213        y *= (t - 1.0).sqrt();
214    }
215
216    // deal with exceptions
217    match (i, u) {
218        (1, 0x40422000) => {
219            y = h!("0x1.c63ffcp+0"); // x = 0x1.844p+1
220        }
221        (3, 0x411c0000) => {
222            y = h!("0x1.7be006p+1"); // x = 0x1.38p+3
223        }
224        (3, 0x415f6000) => {
225            y = h!("0x1.aa0002p+1"); // x = 0x1.becp+3
226        }
227        (4, 0x41bfe000) => {
228            y = h!("0x1.ef5fecp+1"); // x = 0x1.7fcp+4
229        }
230        (4, 0x41c04000) => {
231            y = h!("0x1.ef9ff6p+1"); // x = 0x1.808p+4
232        }
233        (4, 0x41dbc000) => {
234            y = h!("0x1.006016p+2"); // x = 0x1.b78p+4
235        }
236        (4, 0x41d3c000) => {
237            y = h!("0x1.fc002ap+1"); // x = 0x1.a78p+4
238        }
239        (5, 0x42374000) => {
240            y = h!("0x1.21201cp+2"); // x = 0x1.6e8p+5
241        }
242        (5, 0x4245c000) => {
243            y = h!("0x1.260012p+2"); // x = 0x1.8b8p+5
244        }
245        (6, 0x42890000) => {
246            y = h!("0x1.3ae018p+2"); // x = 0x1.12p+6
247        }
248        (6, 0x429d6000) => {
249            y = h!("0x1.43bff4p+2"); // x = 0x1.3acp+6
250        }
251        (7, 0x433d2000) => {
252            y = h!("0x1.7be006p+2"); // x = 0x1.7a4p+7
253        }
254        (8, 0x43884000) => {
255            y = h!("0x1.934004p+2"); // x = 0x1.108p+8
256        }
257        (8, 0x43f9a000) => {
258            y = h!("0x1.ba000ep+2"); // x = 0x1.f34p+8
259        }
260        (9, 0x443b0000) => {
261            y = h!("0x1.d3e00cp+2"); // x = 0x1.76p+9
262        }
263        (12, 0x45b78000) => {
264            y = h!("0x1.2be008p+3"); // x = 0x1.6fp+12
265        }
266        (13, 0x4673a000) => {
267            y = h!("0x1.4b2008p+3"); // x = 0x1.e74p+13
268        }
269        _ => (),
270    }
271
272    y as f16
273}
274
275#[cfg(test)]
276mod tests {
277    #[test]
278    fn exhaustive() {
279        for b in 0..=u16::MAX {
280            let x = f16::from_bits(b);
281            let y1 = super::cr_acoshf16(x);
282            let y2 = core_math::acoshf16(x);
283            assert_eq!(
284                y1.to_bits(),
285                y2.to_bits(),
286                "acoshf16({x} @ {b:#04x}) = ({y1} @ {y1b:#04x}) vs ({y2} @ {y2b:#04x})",
287                y1b = y1.to_bits(),
288                y2b = y2.to_bits()
289            );
290        }
291    }
292
293    #[test]
294    fn edge() {
295        assert_eq!(
296            super::cr_acoshf16(-f16::NAN).to_bits(),
297            (-f16::NAN).to_bits()
298        );
299        assert_eq!(
300            super::cr_acoshf16(-f16::INFINITY).to_bits(),
301            super::snanf16().to_bits() // FIXME: sign
302        );
303        assert_eq!(
304            super::cr_acoshf16(-0.0).to_bits(),
305            super::snanf16().to_bits() // FIXME: sign
306        );
307        assert_eq!(
308            super::cr_acoshf16(0.0).to_bits(),
309            super::snanf16().to_bits() // FIXME: sign
310        );
311        assert_eq!(super::cr_acoshf16(1.0).to_bits(), (0.0_f16).to_bits());
312        assert_eq!(super::cr_acoshf16(f16::INFINITY), f16::INFINITY);
313        assert_eq!(super::cr_acoshf16(f16::NAN).to_bits(), (f16::NAN).to_bits());
314    }
315}
core_math_rs/f16/acosh.rs

core_math_rs/f16/
acosh.rs
