pxfm/sin_cosf/
sinpif.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 8/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::common::{f_fmla, is_integerf};
30use crate::polyeval::f_polyeval5;
31use crate::sin_cosf::argument_reduction_pi::ArgumentReducerPi;
32use crate::sin_cosf::sincosf_eval::{cospif_eval, sinpif_eval, sinpif_eval2};
33
34#[inline(always)]
35fn sinpif_gen_impl(x: f32) -> f32 {
36    let x_abs = x.to_bits() & 0x7fff_ffffu32;
37    let xd = x as f64;
38
39    // |x| <= 1/16
40    if x_abs <= 0x3d80_0000u32 {
41        // |x| < 0.0000009546391
42        if x_abs < 0x3580_2126u32 {
43            if x_abs == 0u32 {
44                // For signed zeros.
45                return x;
46            }
47            const PI: f64 = f64::from_bits(0x400921fb54442d18);
48            const MPI_E3_OVER_6: f64 = f64::from_bits(0xc014abbce625be53);
49
50            // Small values approximated with Taylor poly
51            // x = pi * x - pi^3*x^3/6
52            let x2 = xd * xd;
53            let p = f_fmla(x2, MPI_E3_OVER_6, PI);
54            return (xd * p) as f32;
55        }
56
57        let xsqr = xd * xd;
58
59        /*
60        Generated by Sollya:
61        d = [0, 1/16];
62        f_sin = sin(y*pi)/y;
63        Q = fpminimax(sin(y*pi)/y, [|0, 2, 4, 6, 8|], [|D...|], d, relative, floating);
64
65        See ./notes/sinpif.sollya
66         */
67        let p = f_polyeval5(
68            xsqr,
69            f64::from_bits(0x400921fb54442d18),
70            f64::from_bits(0xc014abbce625bbf2),
71            f64::from_bits(0x400466bc675e116a),
72            f64::from_bits(0xbfe32d2c0b62d41c),
73            f64::from_bits(0x3fb501ec4497cb7d),
74        );
75        return (xd * p) as f32;
76    }
77
78    // Numbers greater or equal to 2^23 are always integers or NaN
79    if x_abs >= 0x4b00_0000u32 || is_integerf(x) {
80        if x_abs >= 0x7f80_0000u32 {
81            return x + f32::NAN;
82        }
83        return if x.is_sign_negative() { -0. } else { 0. };
84    }
85
86    // We're computing sin(y) after argument reduction then return valid value
87    // based on quadrant
88
89    let reducer = ArgumentReducerPi { x: x as f64 };
90    let (y, k) = reducer.reduce_0p25();
91    // Decide based on quadrant what kernel function to use
92    (match k & 3 {
93        0 => sinpif_eval(y),
94        1 => cospif_eval(y),
95        2 => sinpif_eval(-y),
96        _ => -cospif_eval(y),
97    }) as f32
98}
99
100#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
101#[target_feature(enable = "avx", enable = "fma")]
102unsafe fn sinpif_fma_impl(x: f32) -> f32 {
103    let x_abs = x.to_bits() & 0x7fff_ffffu32;
104    let xd = x as f64;
105
106    // |x| <= 1/16
107    if x_abs <= 0x3d80_0000u32 {
108        // |x| < 0.0000009546391
109        if x_abs < 0x3580_2126u32 {
110            if x_abs == 0u32 {
111                // For signed zeros.
112                return x;
113            }
114            const PI: f64 = f64::from_bits(0x400921fb54442d18);
115            const MPI_E3_OVER_6: f64 = f64::from_bits(0xc014abbce625be53);
116
117            // Small values approximated with Taylor poly
118            // x = pi * x - pi^3*x^3/6
119            let x2 = xd * xd;
120            let p = f64::mul_add(x2, MPI_E3_OVER_6, PI);
121            return (xd * p) as f32;
122        }
123
124        let xsqr = xd * xd;
125
126        /*
127        Generated by Sollya:
128        d = [0, 1/16];
129        f_sin = sin(y*pi)/y;
130        Q = fpminimax(sin(y*pi)/y, [|0, 2, 4, 6, 8|], [|D...|], d, relative, floating);
131
132        See ./notes/sinpif.sollya
133         */
134        use crate::polyeval::d_polyeval5;
135        let p = d_polyeval5(
136            xsqr,
137            f64::from_bits(0x400921fb54442d18),
138            f64::from_bits(0xc014abbce625bbf2),
139            f64::from_bits(0x400466bc675e116a),
140            f64::from_bits(0xbfe32d2c0b62d41c),
141            f64::from_bits(0x3fb501ec4497cb7d),
142        );
143        return (xd * p) as f32;
144    }
145
146    // Numbers greater or equal to 2^23 are always integers or NaN
147    if x_abs >= 0x4b00_0000u32 || x.round_ties_even() == x {
148        if x_abs >= 0x7f80_0000u32 {
149            return x + f32::NAN;
150        }
151        return if x.is_sign_negative() { -0. } else { 0. };
152    }
153
154    // We're computing sin(y) after argument reduction then return valid value
155    // based on quadrant
156
157    let reducer = ArgumentReducerPi { x: x as f64 };
158    let (y, k) = reducer.reduce_0p25_fma();
159    // Decide based on quadrant what kernel function to use
160    use crate::sin_cosf::sincosf_eval::{cospif_eval_fma, sinpif_eval_fma};
161    (match k & 3 {
162        0 => sinpif_eval_fma(y),
163        1 => cospif_eval_fma(y),
164        2 => sinpif_eval_fma(-y),
165        _ => -cospif_eval_fma(y),
166    }) as f32
167}
168
169/// Computes sin(PI*x)
170///
171/// Max ULP 0.5
172#[inline]
173pub fn f_sinpif(x: f32) -> f32 {
174    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
175    {
176        sinpif_gen_impl(x)
177    }
178    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
179    {
180        use std::sync::OnceLock;
181        static EXECUTOR: OnceLock<unsafe fn(f32) -> f32> = OnceLock::new();
182        let q = EXECUTOR.get_or_init(|| {
183            if std::arch::is_x86_feature_detected!("avx")
184                && std::arch::is_x86_feature_detected!("fma")
185            {
186                sinpif_fma_impl
187            } else {
188                sinpif_gen_impl
189            }
190        });
191        unsafe { q(x) }
192    }
193}
194
195pub(crate) fn fast_sinpif(x: f32) -> f64 {
196    let x_abs = x.to_bits() & 0x7fff_ffffu32;
197    let xd = x as f64;
198
199    // |x| <= 1/16
200    if x_abs <= 0x3d80_0000u32 {
201        // |x| < 0.0000009546391
202        if x_abs < 0x3580_2126u32 {
203            if x_abs == 0u32 {
204                // For signed zeros.
205                return x as f64;
206            }
207            const PI: f64 = f64::from_bits(0x400921fb54442d18);
208            const MPI_E3_OVER_6: f64 = f64::from_bits(0xc014abbce625be53);
209
210            // Small values approximated with Taylor poly
211            // x = pi * x - pi^3*x^3/6
212            let x2 = xd * xd;
213            let p = f_fmla(x2, MPI_E3_OVER_6, PI);
214            return xd * p;
215        }
216
217        let xsqr = xd * xd;
218
219        /*
220        Generated by Sollya:
221        d = [0, 1/16];
222        f_sin = sin(y*pi)/y;
223        Q = fpminimax(sin(y*pi)/y, [|0, 2, 4, 6, 8|], [|D...|], d, relative, floating);
224
225        See ./notes/sinpif.sollya
226         */
227        let p = f_polyeval5(
228            xsqr,
229            f64::from_bits(0x400921fb54442d18),
230            f64::from_bits(0xc014abbce625bbf2),
231            f64::from_bits(0x400466bc675e116a),
232            f64::from_bits(0xbfe32d2c0b62d41c),
233            f64::from_bits(0x3fb501ec4497cb7d),
234        );
235        return xd * p;
236    }
237
238    let reducer = ArgumentReducerPi { x: x as f64 };
239    let (y, k) = reducer.reduce_0p25();
240    // Decide based on quadrant what kernel function to use
241    match k & 3 {
242        0 => sinpif_eval2(y),
243        1 => cospif_eval(y),
244        2 => sinpif_eval2(-y),
245        _ => -cospif_eval(y),
246    }
247}
248
249#[cfg(test)]
250mod tests {
251    use super::*;
252
253    #[test]
254    fn test_f_sinpif() {
255        assert_eq!(f_sinpif(3.), 0.);
256        assert_eq!(f_sinpif(1.12199515e-7), 3.524852e-7);
257        assert_eq!(f_sinpif(-0.31706), -0.83934295);
258        assert_eq!(f_sinpif(0.30706), 0.8218538);
259        assert_eq!(f_sinpif(115.30706), -0.82185423);
260        assert!(f_sinpif(f32::INFINITY).is_nan());
261        assert!(f_sinpif(f32::NEG_INFINITY).is_nan());
262        assert!(f_sinpif(f32::NAN).is_nan());
263    }
264}
pxfm/sin_cosf/sinpif.rs

pxfm/sin_cosf/
sinpif.rs
