pxfm/
asinf.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 6/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::acosf::poly12;
30use crate::common::{dd_fmlaf, f_fmla};
31
32#[cold]
33fn as_special(x: f32) -> f32 {
34    let t = x.to_bits();
35    let ax = t.wrapping_shl(1);
36    if ax > (0xffu32 << 24) {
37        return x + x;
38    } // nan
39    f32::NAN
40}
41
42/// Computes asin
43///
44/// Max found ULP 0.49999928
45#[inline]
46pub fn f_asinf(x: f32) -> f32 {
47    const PI2: f64 = f64::from_bits(0x3ff921fb54442d18);
48    let xs = x as f64;
49    let mut r;
50    let t = x.to_bits();
51    let ax = t.wrapping_shl(1);
52    if ax > 0x7f << 24 {
53        return as_special(x);
54    }
55    if ax < 0x7ec29000u32 {
56        // |x| < 1.49029
57        if ax < 115 << 24 {
58            // |x| < 0.000244141
59            return dd_fmlaf(x, f64::from_bits(0x3e60000000000000) as f32, x);
60        }
61        const B: [u64; 16] = [
62            0x3ff0000000000005,
63            0x3fc55557aeca105d,
64            0x3fb3314ec3db7d12,
65            0x3fa775738a5a6f92,
66            0x3f75d5f7ce1c8538,
67            0x3fd605c6d58740f0,
68            0xc005728b732d73c6,
69            0x402f152170f151eb,
70            0xc04f962ea3ca992e,
71            0x40671971e17375a0,
72            0xc07860512b4ba230,
73            0x40826a3b8d4bdb14,
74            0xc0836f2ea5698b51,
75            0x407b3d722aebfa2e,
76            0xc066cf89703b1289,
77            0x4041518af6a65e2d,
78        ];
79        let z = xs;
80        let z2 = z * z;
81        let w0 = f_fmla(z2, f64::from_bits(B[1]), f64::from_bits(B[0]));
82        let w1 = f_fmla(z2, f64::from_bits(B[3]), f64::from_bits(B[2]));
83        let w2 = f_fmla(z2, f64::from_bits(B[5]), f64::from_bits(B[4]));
84        let w3 = f_fmla(z2, f64::from_bits(B[7]), f64::from_bits(B[6]));
85        let w4 = f_fmla(z2, f64::from_bits(B[9]), f64::from_bits(B[8]));
86        let w5 = f_fmla(z2, f64::from_bits(B[11]), f64::from_bits(B[10]));
87        let w6 = f_fmla(z2, f64::from_bits(B[13]), f64::from_bits(B[12]));
88        let w7 = f_fmla(z2, f64::from_bits(B[15]), f64::from_bits(B[14]));
89
90        let z4 = z2 * z2;
91        let z8 = z4 * z4;
92        let z16 = z8 * z8;
93
94        r = z
95            * ((f_fmla(z4, w1, w0) + z8 * f_fmla(z4, w3, w2))
96                + z16 * (f_fmla(z4, w5, w4) + z8 * f_fmla(z4, w7, w6)));
97        let ub = r;
98        let lb = r - z * f64::from_bits(0x3e0efa8eb0000000);
99        // Ziv's accuracy test
100        if ub == lb {
101            return ub as f32;
102        }
103    }
104    if ax < (0x7eu32 << 24) {
105        const C: [u64; 12] = [
106            0x3fc555555555529c,
107            0x3fb333333337e0dd,
108            0x3fa6db6db3b4465e,
109            0x3f9f1c72e13ac306,
110            0x3f96e89cebe06bc4,
111            0x3f91c6dcf5289094,
112            0x3f8c6dbbcc7c6315,
113            0x3f88f8dc2615e996,
114            0x3f7a5833b7bf15e8,
115            0x3f943f44ace1665c,
116            0xbf90fb17df881c73,
117            0x3fa07520c026b2d6,
118        ];
119        let z = xs;
120        let z2 = z * z;
121        let c0 = poly12(z2, C);
122        r = z + (z * z2) * c0;
123    } else {
124        if ax == 0x7e55688au32 {
125            return f32::copysign(f64::from_bits(0x3fe75b8a20000000) as f32, x)
126                + f32::copysign(f64::from_bits(0x3e50000000000000) as f32, x);
127        }
128        if ax == 0x7e107434u32 {
129            return f32::copysign(f64::from_bits(0x3fe1f4b640000000) as f32, x)
130                + f32::copysign(f64::from_bits(0x3e50000000000000) as f32, x);
131        }
132        let bx = xs.abs();
133        let z = 1.0 - bx;
134        let s = z.sqrt();
135        const C: [u64; 12] = [
136            0x3ff6a09e667f3bcb,
137            0x3fbe2b7dddff2db9,
138            0x3f9b27247ab42dbc,
139            0x3f802995cc4e0744,
140            0x3f65ffb0276ec8ea,
141            0x3f5033885a928dec,
142            0x3f3911f2be23f8c7,
143            0x3f24c3c55d2437fd,
144            0x3f0af477e1d7b461,
145            0x3f0abd6bdff67dcb,
146            0xbef1717e86d0fa28,
147            0x3ef6ff526de46023,
148        ];
149        r = PI2 - s * poly12(z, C);
150        r = f64::copysign(r, xs);
151    }
152    r as f32
153}
154
155#[cfg(test)]
156mod tests {
157    use super::*;
158
159    #[test]
160    fn test_asinf() {
161        assert_eq!(f_asinf(-0.5), -std::f32::consts::FRAC_PI_6);
162        assert_eq!(f_asinf(0.5), std::f32::consts::FRAC_PI_6);
163        assert!(f_asinf(7.).is_nan());
164    }
165}
pxfm/asinf.rs

pxfm/
asinf.rs
