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/// Computes sin(PI*x)
35///
36/// Max ULP 0.5
37#[inline]
38pub fn f_sinpif(x: f32) -> f32 {
39    let x_abs = x.to_bits() & 0x7fff_ffffu32;
40    let xd = x as f64;
41
42    // |x| <= 1/16
43    if x_abs <= 0x3d80_0000u32 {
44        // |x| < 0.0000009546391
45        if x_abs < 0x3580_2126u32 {
46            if x_abs == 0u32 {
47                // For signed zeros.
48                return x;
49            }
50            const PI: f64 = f64::from_bits(0x400921fb54442d18);
51            const MPI_E3_OVER_6: f64 = f64::from_bits(0xc014abbce625be53);
52
53            // Small values approximated with Taylor poly
54            // x = pi * x - pi^3*x^3/6
55            let x2 = xd * xd;
56            let p = f_fmla(x2, MPI_E3_OVER_6, PI);
57            return (xd * p) as f32;
58        }
59
60        let xsqr = xd * xd;
61
62        /*
63        Generated by Sollya:
64        d = [0, 1/16];
65        f_sin = sin(y*pi)/y;
66        Q = fpminimax(sin(y*pi)/y, [|0, 2, 4, 6, 8|], [|D...|], d, relative, floating);
67
68        See ./notes/sinpif.sollya
69         */
70        let p = f_polyeval5(
71            xsqr,
72            f64::from_bits(0x400921fb54442d18),
73            f64::from_bits(0xc014abbce625bbf2),
74            f64::from_bits(0x400466bc675e116a),
75            f64::from_bits(0xbfe32d2c0b62d41c),
76            f64::from_bits(0x3fb501ec4497cb7d),
77        );
78        return (xd * p) as f32;
79    }
80
81    // Numbers greater or equal to 2^23 are always integers or NaN
82    if x_abs >= 0x4b00_0000u32 || is_integerf(x) {
83        if x_abs >= 0x7f80_0000u32 {
84            return x + f32::NAN;
85        }
86        return if x.is_sign_negative() { -0. } else { 0. };
87    }
88
89    // We're computing sin(y) after argument reduction then return valid value
90    // based on quadrant
91
92    let reducer = ArgumentReducerPi { x: x as f64 };
93    let (y, k) = reducer.reduce_0p25();
94    // Decide based on quadrant what kernel function to use
95    (match k & 3 {
96        0 => sinpif_eval(y),
97        1 => cospif_eval(y),
98        2 => sinpif_eval(-y),
99        _ => -cospif_eval(y),
100    }) as f32
101}
102
103pub(crate) fn fast_sinpif(x: f32) -> f64 {
104    let x_abs = x.to_bits() & 0x7fff_ffffu32;
105    let xd = x as f64;
106
107    // |x| <= 1/16
108    if x_abs <= 0x3d80_0000u32 {
109        // |x| < 0.0000009546391
110        if x_abs < 0x3580_2126u32 {
111            if x_abs == 0u32 {
112                // For signed zeros.
113                return x as f64;
114            }
115            const PI: f64 = f64::from_bits(0x400921fb54442d18);
116            const MPI_E3_OVER_6: f64 = f64::from_bits(0xc014abbce625be53);
117
118            // Small values approximated with Taylor poly
119            // x = pi * x - pi^3*x^3/6
120            let x2 = xd * xd;
121            let p = f_fmla(x2, MPI_E3_OVER_6, PI);
122            return xd * p;
123        }
124
125        let xsqr = xd * xd;
126
127        /*
128        Generated by Sollya:
129        d = [0, 1/16];
130        f_sin = sin(y*pi)/y;
131        Q = fpminimax(sin(y*pi)/y, [|0, 2, 4, 6, 8|], [|D...|], d, relative, floating);
132
133        See ./notes/sinpif.sollya
134         */
135        let p = f_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;
144    }
145
146    let reducer = ArgumentReducerPi { x: x as f64 };
147    let (y, k) = reducer.reduce_0p25();
148    // Decide based on quadrant what kernel function to use
149    match k & 3 {
150        0 => sinpif_eval2(y),
151        1 => cospif_eval(y),
152        2 => sinpif_eval2(-y),
153        _ => -cospif_eval(y),
154    }
155}
156
157#[cfg(test)]
158mod tests {
159    use super::*;
160
161    #[test]
162    fn test_f_sinpif() {
163        assert_eq!(f_sinpif(3.), 0.);
164        assert_eq!(f_sinpif(1.12199515e-7), 3.524852e-7);
165        assert_eq!(f_sinpif(-0.31706), -0.83934295);
166        assert_eq!(f_sinpif(0.30706), 0.8218538);
167        assert_eq!(f_sinpif(115.30706), -0.82185423);
168        assert!(f_sinpif(f32::INFINITY).is_nan());
169        assert!(f_sinpif(f32::NEG_INFINITY).is_nan());
170        assert!(f_sinpif(f32::NAN).is_nan());
171    }
172}
pxfm/sin_cosf/sinpif.rs

pxfm/sin_cosf/
sinpif.rs
