core_math_rs/f16/
asin.rs

1/* Correctly-rounded 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
29use super::utils::snanf16;
30
31// the following polynomials were generated using Sollya (cf asin.sollya)
32
33/* degree-5 minimax polynomial for asin(x) over [0,0.25], with relative error
34   bounded by 2^-20.817, with coefficients of odd degree only, and degree-1
35   coefficient forced to 1
36*/
37const P0: [f32; 3] = [1.0, h!("0x1.552e8ep-3"), h!("0x1.43696ep-4")];
38
39/* degree-7 minimax polynomial for asin(x) over [0.25,0.5], with relative error
40   bounded by 2^-19.202, with coefficients of odd degree only, and degree-1
41   coefficient forced to 1
42*/
43const P1: [f32; 4] = [
44    1.0,
45    h!("0x1.55a7ep-3"),
46    h!("0x1.258e54p-4"),
47    h!("0x1.08e44p-4"),
48];
49
50/// Correctly-rounded arc-sine for binary16 value.
51pub fn cr_asinf16(x: f16) -> f16 {
52    let v = f32::from(x);
53    let u = v.to_bits();
54    let mut au = u & 0x7fffffff;
55
56    const HALF_PI: f32 = h!("0x1.921fb6p+0");
57
58    if au >= 0x3f800000 {
59        // NaN, Inf, or |x| >= 1
60        if au == 0x3f800000 {
61            // |x| = 1
62            return if u == 0x3f800000 { HALF_PI } else { -HALF_PI } as f16;
63        }
64        if (au >> 23) == 0x3ff && ((au & 0x7fffff) != 0) {
65            // qNaN or sNaN
66            return x + x;
67        }
68        return snanf16(); // will signal invalid and return sNaN
69    }
70
71    let mut t = v;
72
73    /* for |x| < 0.5 we don't reduce t into -t, since the polynomials
74    P0 and P1 are valid also on the negative side */
75
76    let reduce = au >= 0x3f000000; // |x| >= 0.5
77    if reduce {
78        if (u >> 31) != 0 {
79            // x < 0
80            t = -t;
81        }
82        // argument reduction: asin(x) = pi/2 - 2*asin(sqrt((1-x)/2))
83        t = ((1.0 - t) * 0.5).sqrt();
84        au = t.to_bits() & 0x7fffffff;
85    }
86
87    // now 0 <= t <= 0.5
88    let tt = t * t;
89
90    let mut y;
91    if au < 0x3e800000 {
92        // 0 <= t < 0.25
93        /* For |x| <= 0x1.71p-5, asin(x) rounds to x to nearest,
94        we deal with that case separately, so that for x subnormal
95        and a power of two, we get an underflow. */
96        if (au <= 0x3d388000) && !reduce {
97            if au == 0 {
98                return x;
99            }
100            return if au == u {
101                v + h!("0x1p-26")
102            } else {
103                v - h!("0x1p-26")
104            } as f16;
105        }
106        if au == 0x3dd30000 {
107            // |x| = 0x1.a6p-4
108            return if au == u {
109                h!("0x1.a6c00ap-4")
110            } else {
111                h!("-0x1.a6c00ap-4")
112            } as f16;
113        }
114        if au == 0x3d688000 {
115            // |x| = 0x1.d1p-5
116            return if au == u {
117                h!("0x1.d14004p-5")
118            } else {
119                h!("-0x1.d14004p-5")
120            } as f16;
121        }
122        if au == 0x3dfa0000 {
123            // |x| = 0x1.f4p-4
124            return if au == u {
125                h!("0x1.f5400ap-4")
126            } else {
127                h!("-0x1.f5400ap-4")
128            } as f16;
129        }
130        /* Warning for rounding toward -Inf: let P0(t) = t*q(t). If we first
131        compute q(t) and then multiply by t, for tiny t and rounding we will
132        get q(t)=1, and then t, whereas the correct result is nextbelow(t). */
133        let c1 = f32::mul_add(P0[2], tt, P0[1]);
134        y = f32::mul_add(c1, tt * t, t);
135    } else {
136        // 0.25 <= t <= 0.5
137        if au == 0x3eb24000 {
138            // |x| = 0x1.648p-2
139            return if au == u {
140                h!("0x1.6c2012p-2")
141            } else {
142                h!("-0x1.6c2012p-2")
143            } as f16;
144        }
145        if au == 0x3ed96000 {
146            // |x| = 0x1.b2cp-2
147            return if au == u {
148                h!("0x1.c0fffap-2")
149            } else {
150                h!("-0x1.c0fffap-2")
151            } as f16;
152        }
153        if au == 0x3ef0a000 {
154            // |x| = 0x1.e14p-2
155            return if au == u {
156                h!("0x1.f4fffp-2")
157            } else {
158                h!("-0x1.f4fffp-2")
159            } as f16;
160        }
161        let c5 = f32::mul_add(P1[3], tt, P1[2]);
162        let c1 = f32::mul_add(P1[1], tt, P1[0]);
163        y = f32::mul_add(c5, tt * tt, c1);
164        y *= t;
165    }
166
167    if reduce {
168        // argument reconstruction
169        y = HALF_PI - 2.0 * y;
170        if (u >> 31) != 0 {
171            // x < 0
172            y = -y;
173        }
174    }
175
176    y as f16
177}
178
179#[cfg(test)]
180mod tests {
181    #[test]
182    fn exhaustive() {
183        for b in 0..=u16::MAX {
184            let x = f16::from_bits(b);
185            let y1 = super::cr_asinf16(x);
186            let y2 = core_math::asinf16(x);
187            assert_eq!(
188                y1.to_bits(),
189                y2.to_bits(),
190                "asinf16({x} @ {b:#04x}) = ({y1} @ {y1b:#04x}) vs ({y2} @ {y2b:#04x})",
191                y1b = y1.to_bits(),
192                y2b = y2.to_bits()
193            );
194        }
195    }
196
197    #[test]
198    fn edge() {
199        assert_eq!(
200            super::cr_asinf16(-f16::NAN).to_bits(),
201            super::snanf16().to_bits() // FIXME: sign
202        );
203        assert_eq!(
204            super::cr_asinf16(-f16::INFINITY).to_bits(),
205            super::snanf16().to_bits() // FIXME: sign
206        );
207        assert_eq!(super::cr_asinf16(-1.0), -std::f16::consts::FRAC_PI_2);
208        assert_eq!(super::cr_asinf16(-0.0).to_bits(), (-0.0_f16).to_bits());
209        assert_eq!(super::cr_asinf16(0.0).to_bits(), (0.0_f16).to_bits());
210        assert_eq!(super::cr_asinf16(1.0), std::f16::consts::FRAC_PI_2);
211        assert_eq!(
212            super::cr_asinf16(f16::INFINITY).to_bits(),
213            super::snanf16().to_bits() // FIXME: sign
214        );
215        assert_eq!(
216            super::cr_asinf16(f16::NAN).to_bits(),
217            super::snanf16().to_bits() // FIXME: sign
218        );
219    }
220}
core_math_rs/f16/asin.rs

core_math_rs/f16/
asin.rs
