core_math_rs/f16/
atan.rs

1/* Correctly-rounded arc-tangent 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 atan.sollya)
30
31/* Degree-7 minimax polynomial for atan(x) over [0,0.25], with relative error
32   bounded by 2^-25.419, with coefficients of odd degree only, and degree-1
33   coefficient forced to 1. Coefficients were later optimized to reduce the
34   number of exceptions.
35*/
36const P0: [f32; 4] = [
37    h!("0x1.fffffcp-1"),
38    h!("-0x1.55546cp-2"),
39    h!("0x1.98d0ep-3"),
40    h!("-0x1.0c7c54p-3"),
41];
42
43/* degree-4 minimax polynomial for atan(x) over [0.25,0.5], with relative error
44bounded by 2^-22.573 */
45const P1: [f32; 5] = [
46    h!("0x1.411612p-14"),
47    h!("0x1.ff076cp-1"),
48    h!("0x1.1ee64cp-6"),
49    h!("-0x1.a6fc96p-2"),
50    h!("0x1.81dd3ep-3"),
51];
52
53/* degree-4 minimax polynomial for atan(x) over [0.5,0.75], with relative error
54bounded by 2^-21.757 */
55const P2: [f32; 5] = [
56    h!("-0x1.95964cp-8"),
57    h!("0x1.0bad76p+0"),
58    h!("-0x1.e652e2p-4"),
59    h!("-0x1.e70ab4p-3"),
60    h!("0x1.a5eb88p-4"),
61];
62
63/* degree-4 minimax polynomial for atan(x) over [0.75,1], with relative error
64bounded by 2^-23.027 */
65const P3: [f32; 5] = [
66    h!("-0x1.f75b16p-6"),
67    h!("0x1.2d3d1p+0"),
68    h!("-0x1.882376p-2"),
69    h!("0x1.d18c96p-13"),
70    h!("0x1.6aa268p-6"),
71];
72
73/// Correctly-rounded arc-tangent for binary16 value.
74pub fn cr_atanf16(x: f16) -> f16 {
75    let v = f32::from(x);
76    let u = v.to_bits();
77    let au = u & 0x7fffffff;
78
79    const HALF_PI: f32 = h!("0x1.921fb6p+0");
80
81    if au >= 0x7f800000 {
82        // NaN or Inf
83        if au == 0x7f800000 {
84            // +/-Inf
85            return if u == 0x7f800000 { HALF_PI } else { -HALF_PI } as f16;
86        }
87        return x + x; // will signal invalid for sNaN and return qNaN
88    }
89
90    let mut t = v;
91
92    // for x < 0 we use atan(-x) = -atan(x)
93    let neg = (u >> 31) != 0;
94    const NEG: [f32; 2] = [1.0, -1.0];
95    let s = NEG[usize::from(neg)];
96    t *= s;
97
98    // now t >= 0
99
100    // for x > 1 we use atan(x) = pi/2 - atan(1/x)
101    let reduce = au > 0x3f800000;
102    if reduce {
103        t = 1.0 / t;
104    }
105
106    // now 0 <= t <= 1
107
108    let tt = t * t;
109    let mut y;
110    if t <= 0.25 {
111        // for |x| < 0x1.d14p-6, atan(x) rounds to x to nearest
112        if !reduce && (t <= h!("0x1.d14p-6")) {
113            if au == 0 {
114                return x; // x = 0
115            }
116            t *= s;
117            return f32::mul_add(t, h!("-0x1p-23"), t) as f16;
118        }
119
120        // deal with exceptional cases
121        if au == 0x3e56a000 {
122            // |x| = 0x1.ad4p-3
123            return if au == u {
124                h!("0x1.a72002p-3")
125            } else {
126                h!("-0x1.a72002p-3")
127            } as f16;
128        }
129        if au == 0x4115c000 {
130            // |x| = 0x1.2b8p+3
131            return if au == u {
132                h!("0x1.76dffep+0")
133            } else {
134                h!("-0x1.76dffep+0")
135            } as f16;
136        }
137        if au == 0x42c32000 {
138            // |x| = 0x1.864p+6
139            return if au == u {
140                h!("0x1.8f7ffep+0")
141            } else {
142                h!("-0x1.8f7ffep+0")
143            } as f16;
144        }
145
146        let p = P0;
147        let c5 = f32::mul_add(p[3], tt, p[2]);
148        let c1 = f32::mul_add(p[1], tt, p[0]);
149        let c1 = f32::mul_add(c5, tt * tt, c1);
150        y = t * c1;
151    } else {
152        let p = if t <= 0.5 {
153            P1
154        } else if t <= 0.75 {
155            P2
156        } else {
157            P3
158        };
159        let c3 = f32::mul_add(p[4], t, p[3]);
160        let c2 = f32::mul_add(c3, t, p[2]);
161        y = f32::mul_add(p[1], t, p[0]);
162        y = f32::mul_add(c2, tt, y);
163    }
164
165    if reduce {
166        y = HALF_PI - y; // argument reconstruction
167    }
168
169    if neg {
170        y = -y;
171    }
172
173    y as f16
174}
175
176#[cfg(test)]
177mod tests {
178    #[test]
179    fn exhaustive() {
180        for b in 0..=u16::MAX {
181            let x = f16::from_bits(b);
182            let y1 = super::cr_atanf16(x);
183            let y2 = core_math::atanf16(x);
184            assert_eq!(
185                y1.to_bits(),
186                y2.to_bits(),
187                "atanf16({x} @ {b:#04x}) = ({y1} @ {y1b:#04x}) vs ({y2} @ {y2b:#04x})",
188                y1b = y1.to_bits(),
189                y2b = y2.to_bits()
190            );
191        }
192    }
193
194    #[test]
195    fn edge() {
196        assert_eq!(
197            super::cr_atanf16(-f16::NAN).to_bits(),
198            (-f16::NAN).to_bits()
199        );
200        assert_eq!(
201            super::cr_atanf16(-f16::INFINITY),
202            -std::f16::consts::FRAC_PI_2,
203        );
204        assert_eq!(super::cr_atanf16(-0.0).to_bits(), (-0.0_f16).to_bits());
205        assert_eq!(super::cr_atanf16(0.0).to_bits(), (0.0_f16).to_bits());
206        assert_eq!(
207            super::cr_atanf16(f16::INFINITY),
208            std::f16::consts::FRAC_PI_2
209        );
210        assert_eq!(super::cr_atanf16(f16::NAN).to_bits(), (f16::NAN).to_bits());
211    }
212}
core_math_rs/f16/atan.rs

core_math_rs/f16/
atan.rs
