core_math_rs/f16/
acos.rs

1/* Correctly-rounded 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
27#![expect(clippy::approx_constant)]
28
29use hexf::hexf32 as h;
30
31use super::utils::snanf16;
32
33// the following polynomials were generated using Sollya (cf acos.sollya)
34
35/* degree-4 minimax polynomial for acos(x) over [0,0.25], with relative error
36   bounded by 2^-22.943, manually optimized to reduce the number of exceptions
37*/
38const P0: [f32; 5] = [
39    h!("0x1.921fb4p0"),
40    h!("-0x1.fffb44p-1"),
41    h!("-0x1.25e6b8p-10"),
42    h!("-0x1.3cc114p-3"),
43    h!("-0x1.a85b22p-5"),
44];
45
46/* degree-4 minimax polynomial for acos(x) over [0.25,0.5], with relative error
47   bounded by 2^-20.789, manually optimized to reduce the number of exceptions
48*/
49const P1: [f32; 5] = [
50    h!("0x1.91b678p0"),
51    h!("-0x1.f515cap-1"),
52    h!("-0x1.bd043ap-4"),
53    h!("0x1.7e2d5ap-4"),
54    h!("-0x1.190806p-2"),
55];
56
57/* degree-4 minimax polynomial for acos(x)/sqrt(1-x) over [0.5,1],
58with relative error bounded by 2^-23.583 */
59const P2: [f32; 5] = [
60    h!("0x1.91fa1cp0"),
61    h!("-0x1.ae5c5ep-3"),
62    h!("0x1.31640cp-4"),
63    h!("-0x1.98038p-6"),
64    h!("0x1.251b5p-8"),
65];
66
67/// Correctly-rounded arc-cosine for binary16 value.
68pub fn cr_acosf16(x: f16) -> f16 {
69    let v = f32::from(x);
70    let u = v.to_bits();
71    let au = u & 0x7fffffff;
72
73    if au >= 0x3f800000 {
74        // NaN, Inf, or |x| >= 1
75        if au == 0x3f800000 {
76            if u == 0x3f800000 {
77                return 0.0;
78            }
79            return h!("0x1.921fb6p+1") as f16; // TODO
80        }
81
82        if (au >> 23) == 0x3ff && ((au & 0x7fffff) != 0) {
83            // qNaN or sNaN
84            return x;
85        }
86        return snanf16(); // will signal invalid and return sNaN
87    }
88
89    let mut t = v;
90    let tt = t * t;
91
92    if (u >> 31) != 0 {
93        // x < 0
94        t = -t;
95    }
96
97    let c1;
98    let c3;
99    let mut y;
100    if au < 0x3e800000 {
101        // |x| < 0.25
102        c1 = f32::mul_add(P0[2], t, P0[1]);
103        c3 = f32::mul_add(P0[4], t, P0[3]);
104        y = f32::mul_add(c3, tt, c1);
105        y = f32::mul_add(y, t, P0[0]);
106        if (u >> 31) != 0
107        // x < 0
108        {
109            y = h!("0x1.921fb4p+1") - y;
110            /* below we deal with a few exceptions that we were unable to remove
111            by tuning the coefficients of p0 */
112            if (0x36960000..=0x369a0000).contains(&au) {
113                y = h!("0x1.922002p+0"); // x = -0x1.2cp-18 or -0x1.3p-18 or -0x1.34p-18
114            }
115            if u == 0xbcf80000 {
116                y = h!("0x1.99e002p+0"); // x = -0x1.fp-6
117            }
118            if u == 0xbcfc0000 {
119                y = h!("0x1.9a0006p+0"); // -0x1.f8p-6
120            }
121        }
122    } else if au < 0x3f000000 {
123        // 0.25 <= |x| < 0.5
124        c1 = f32::mul_add(P1[2], t, P1[1]);
125        c3 = f32::mul_add(P1[4], t, P1[3]);
126        y = f32::mul_add(c3, tt, c1);
127        y = f32::mul_add(y, t, P1[0]);
128        if (u >> 31) != 0 {
129            // x < 0
130            y = h!("0x1.921fb4p+1") - y;
131        }
132    } else {
133        // 0.5 <= |x| <= 1
134        c1 = f32::mul_add(P2[2], t, P2[1]);
135        c3 = f32::mul_add(P2[4], t, P2[3]);
136        y = f32::mul_add(c3, tt, c1);
137        y = f32::mul_add(y, t, P2[0]);
138        y *= (1.0 - t).sqrt();
139        if (u >> 31) != 0 {
140            // x < 0
141            y = h!("0x1.921fb4p+1") - y;
142        }
143    }
144
145    y as f16
146}
147
148#[cfg(test)]
149mod tests {
150    #[test]
151    fn exhaustive() {
152        for b in 0..=u16::MAX {
153            let x = f16::from_bits(b);
154            let y1 = super::cr_acosf16(x);
155            let y2 = core_math::acosf16(x);
156            assert_eq!(
157                y1.to_bits(),
158                y2.to_bits(),
159                "acosf16({x} @ {b:#04x}) = ({y1} @ {y1b:#04x}) vs ({y2} @ {y2b:#04x})",
160                y1b = y1.to_bits(),
161                y2b = y2.to_bits()
162            );
163        }
164    }
165
166    #[test]
167    fn edge() {
168        assert_eq!(
169            super::cr_acosf16(-f16::NAN).to_bits(),
170            super::snanf16().to_bits() // FIXME: sign
171        );
172        assert_eq!(
173            super::cr_acosf16(-f16::INFINITY).to_bits(),
174            super::snanf16().to_bits() // FIXME: sign
175        );
176        assert_eq!(super::cr_acosf16(-1.0), std::f16::consts::PI);
177        assert_eq!(super::cr_acosf16(-0.0), std::f16::consts::FRAC_PI_2);
178        assert_eq!(super::cr_acosf16(0.0), std::f16::consts::FRAC_PI_2);
179        assert_eq!(super::cr_acosf16(1.0).to_bits(), (0.0_f16).to_bits());
180        assert_eq!(
181            super::cr_acosf16(f16::INFINITY).to_bits(),
182            super::snanf16().to_bits() // FIXME: sign
183        );
184        assert_eq!(
185            super::cr_acosf16(f16::NAN).to_bits(),
186            super::snanf16().to_bits() // FIXME: sign
187        );
188    }
189}
core_math_rs/f16/acos.rs

core_math_rs/f16/
acos.rs
