pxfm/sin_cosf/sincospif.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, is_odd_integerf};
30use crate::polyeval::f_polyeval5;
31use crate::sin_cosf::sincosf_eval::sincospif_eval;
32
33#[inline(always)]
34fn sincospif_gen_impl(x: f32) -> (f32, f32) {
35 let x_abs = x.to_bits() & 0x7fff_ffffu32;
36 let xd = x as f64;
37
38 // |x| <= 1/16
39 if x_abs <= 0x3d80_0000u32 {
40 // |x| < 0.00000009546391
41 if x_abs < 0x38a2_f984u32 {
42 const PI: f64 = f64::from_bits(0x400921fb54442d18);
43 const MPI_E3_OVER_6: f64 = f64::from_bits(0xc014abbce625be53);
44
45 // Small values approximated with Taylor poly for sine
46 // x = pi * x - pi^3*x^3/6
47 let x2 = xd * xd;
48 let p = f_fmla(x2, MPI_E3_OVER_6, PI);
49 let sf = (xd * p) as f32;
50 #[cfg(any(
51 all(
52 any(target_arch = "x86", target_arch = "x86_64"),
53 target_feature = "fma"
54 ),
55 target_arch = "aarch64"
56 ))]
57 {
58 use crate::common::f_fmlaf;
59 let cs = f_fmlaf(x, f32::from_bits(0xb3000000), 1.);
60 return (sf, cs);
61 }
62 #[cfg(not(any(
63 all(
64 any(target_arch = "x86", target_arch = "x86_64"),
65 target_feature = "fma"
66 ),
67 target_arch = "aarch64"
68 )))]
69 {
70 let cs = f_fmla(xd, f64::from_bits(0xbe60000000000000), 1.) as f32;
71 return (sf, cs);
72 }
73 }
74
75 // Cos(x*PI)
76 // Generated poly by Sollya:
77 // d = [0, 1/16];
78 // f_cos = cos(y*pi);
79 // Q = fpminimax(f_cos, [|0, 2, 4, 6, 8|], [|D...|], d, relative, floating);
80 //
81 // See ./notes/cospif.sollya
82
83 let x2 = xd * xd;
84 let cs = f_polyeval5(
85 x2,
86 f64::from_bits(0x3ff0000000000000),
87 f64::from_bits(0xc013bd3cc9be43f7),
88 f64::from_bits(0x40103c1f08091fe0),
89 f64::from_bits(0xbff55d3ba3d94835),
90 f64::from_bits(0x3fce173c2a00e74e),
91 ) as f32;
92 /*
93 Sin(x*PI)
94 Generated by Sollya:
95 d = [0, 1/16];
96 f_sin = sin(y*pi)/y;
97 Q = fpminimax(sin(y*pi)/y, [|0, 2, 4, 6, 8|], [|D...|], d, relative, floating);
98
99 See ./notes/sinpif.sollya
100 */
101 let p = f_polyeval5(
102 x2,
103 f64::from_bits(0x400921fb54442d18),
104 f64::from_bits(0xc014abbce625bbf2),
105 f64::from_bits(0x400466bc675e116a),
106 f64::from_bits(0xbfe32d2c0b62d41c),
107 f64::from_bits(0x3fb501ec4497cb7d),
108 );
109 let sf = (xd * p) as f32;
110
111 return (sf, cs);
112 }
113
114 // Numbers greater or equal to 2^23 are always integers or NaN
115 if x_abs >= 0x4b00_0000u32 || is_integerf(x) {
116 if x_abs >= 0x7f80_0000u32 {
117 return (x + f32::NAN, x + f32::NAN);
118 }
119 static SF: [f32; 2] = [0., -0.];
120 let sf = SF[x.is_sign_negative() as usize];
121 if x_abs < 0x4b80_0000u32 {
122 static CF: [f32; 2] = [1., -1.];
123 let cs = CF[is_odd_integerf(x) as usize];
124 return (sf, cs);
125 }
126 return (sf, 1.);
127 }
128
129 // Formula:
130 // cos(x) = cos((k + y)*pi/32)
131 // = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
132 // The values of sin(k*pi/32) and cos(k*pi/32) for k = 0..63 are precomputed
133 // and stored using a vector of 32 doubles. Sin(y*pi/32) and cos(y*pi/32) are
134 // computed using degree-7 and degree-6 minimax polynomials generated by
135 // Sollya respectively.
136 // Combine the results with the sine of sum formula:
137 // cos(x) = cos((k + y)*pi/32)
138 // = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
139 // = cosm1_y * cos_k + sin_y * sin_k
140 // = (cosm1_y * cos_k + cos_k) + sin_y * sin_k
141
142 // sin(x) = sin((k + y)*pi/32)
143 // = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
144
145 let rs = sincospif_eval(xd);
146 let cs = f_fmla(rs.sin_y, -rs.sin_k, f_fmla(rs.cosm1_y, rs.cos_k, rs.cos_k)) as f32;
147 let sf = f_fmla(rs.sin_y, rs.cos_k, f_fmla(rs.cosm1_y, rs.sin_k, rs.sin_k)) as f32;
148 (sf, cs)
149}
150
151#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
152#[target_feature(enable = "avx", enable = "fma")]
153unsafe fn sincospif_fma_impl(x: f32) -> (f32, f32) {
154 let x_abs = x.to_bits() & 0x7fff_ffffu32;
155 let xd = x as f64;
156
157 // |x| <= 1/16
158 if x_abs <= 0x3d80_0000u32 {
159 // |x| < 0.00000009546391
160 if x_abs < 0x38a2_f984u32 {
161 const PI: f64 = f64::from_bits(0x400921fb54442d18);
162 const MPI_E3_OVER_6: f64 = f64::from_bits(0xc014abbce625be53);
163
164 // Small values approximated with Taylor poly for sine
165 // x = pi * x - pi^3*x^3/6
166 let x2 = xd * xd;
167 let p = f64::mul_add(x2, MPI_E3_OVER_6, PI);
168 let sf = (xd * p) as f32;
169 let cs = f32::mul_add(x, f32::from_bits(0xb3000000), 1.);
170 return (sf, cs);
171 }
172
173 use crate::polyeval::d_polyeval5;
174
175 // Cos(x*PI)
176 // Generated poly by Sollya:
177 // d = [0, 1/16];
178 // f_cos = cos(y*pi);
179 // Q = fpminimax(f_cos, [|0, 2, 4, 6, 8|], [|D...|], d, relative, floating);
180 //
181 // See ./notes/cospif.sollya
182
183 let x2 = xd * xd;
184 let cs = d_polyeval5(
185 x2,
186 f64::from_bits(0x3ff0000000000000),
187 f64::from_bits(0xc013bd3cc9be43f7),
188 f64::from_bits(0x40103c1f08091fe0),
189 f64::from_bits(0xbff55d3ba3d94835),
190 f64::from_bits(0x3fce173c2a00e74e),
191 ) as f32;
192 /*
193 Sin(x*PI)
194 Generated by Sollya:
195 d = [0, 1/16];
196 f_sin = sin(y*pi)/y;
197 Q = fpminimax(sin(y*pi)/y, [|0, 2, 4, 6, 8|], [|D...|], d, relative, floating);
198
199 See ./notes/sinpif.sollya
200 */
201 let p = d_polyeval5(
202 x2,
203 f64::from_bits(0x400921fb54442d18),
204 f64::from_bits(0xc014abbce625bbf2),
205 f64::from_bits(0x400466bc675e116a),
206 f64::from_bits(0xbfe32d2c0b62d41c),
207 f64::from_bits(0x3fb501ec4497cb7d),
208 );
209 let sf = (xd * p) as f32;
210
211 return (sf, cs);
212 }
213
214 // Numbers greater or equal to 2^23 are always integers or NaN
215 if x_abs >= 0x4b00_0000u32 || x.round_ties_even() == x {
216 if x_abs >= 0x7f80_0000u32 {
217 return (x + f32::NAN, x + f32::NAN);
218 }
219 static SF: [f32; 2] = [0., -0.];
220 let sf = SF[x.is_sign_negative() as usize];
221 if x_abs < 0x4b80_0000u32 {
222 static CF: [f32; 2] = [1., -1.];
223 let is_odd_integer = unsafe { (x.to_int_unchecked::<i32>() & 1) != 0 };
224 let cs = CF[is_odd_integer as usize];
225 return (sf, cs);
226 }
227 return (sf, 1.);
228 }
229
230 // Formula:
231 // cos(x) = cos((k + y)*pi/32)
232 // = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
233 // The values of sin(k*pi/32) and cos(k*pi/32) for k = 0..63 are precomputed
234 // and stored using a vector of 32 doubles. Sin(y*pi/32) and cos(y*pi/32) are
235 // computed using degree-7 and degree-6 minimax polynomials generated by
236 // Sollya respectively.
237 // Combine the results with the sine of sum formula:
238 // cos(x) = cos((k + y)*pi/32)
239 // = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
240 // = cosm1_y * cos_k + sin_y * sin_k
241 // = (cosm1_y * cos_k + cos_k) + sin_y * sin_k
242
243 // sin(x) = sin((k + y)*pi/32)
244 // = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
245 use crate::sin_cosf::sincosf_eval::sincospif_eval_fma;
246 let rs = sincospif_eval_fma(xd);
247 let cs = f64::mul_add(
248 rs.sin_y,
249 -rs.sin_k,
250 f64::mul_add(rs.cosm1_y, rs.cos_k, rs.cos_k),
251 ) as f32;
252 let sf = f64::mul_add(
253 rs.sin_y,
254 rs.cos_k,
255 f64::mul_add(rs.cosm1_y, rs.sin_k, rs.sin_k),
256 ) as f32;
257 (sf, cs)
258}
259
260/// Computes sin(x) and cos(x) at the same time
261///
262/// Max ULP 0.5
263#[inline]
264pub fn f_sincospif(x: f32) -> (f32, f32) {
265 #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
266 {
267 sincospif_gen_impl(x)
268 }
269 #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
270 {
271 use std::sync::OnceLock;
272 static EXECUTOR: OnceLock<unsafe fn(f32) -> (f32, f32)> = OnceLock::new();
273 let q = EXECUTOR.get_or_init(|| {
274 if std::arch::is_x86_feature_detected!("avx")
275 && std::arch::is_x86_feature_detected!("fma")
276 {
277 sincospif_fma_impl
278 } else {
279 sincospif_gen_impl
280 }
281 });
282 unsafe { q(x) }
283 }
284}
285
286#[cfg(test)]
287mod tests {
288 use super::*;
289 use crate::{f_cospif, f_sinpif};
290
291 #[test]
292 fn test_sincospif() {
293 let v0 = f_sincospif(-5.);
294 assert_eq!(v0.0, f_sinpif(-5.));
295 assert_eq!(v0.1, f_cospif(-5.));
296
297 let v0 = f_sincospif(-4.);
298 assert_eq!(v0.0, f_sinpif(-4.));
299 assert_eq!(v0.1, f_cospif(-4.));
300
301 let v0 = f_sincospif(4.);
302 assert_eq!(v0.0, f_sinpif(4.));
303 assert_eq!(v0.1, f_cospif(4.));
304
305 let v0 = f_sincospif(-8489897.0);
306 assert_eq!(v0.0, f_sinpif(-8489897.0));
307 assert_eq!(v0.1, f_cospif(-8489897.0));
308
309 let v1 = f_sincospif(3.23);
310 assert_eq!(v1.0, f_sinpif(3.23));
311 assert_eq!(v1.1, f_cospif(3.23));
312 }
313}