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/// Computes sin(x) and cos(x) at the same time
34///
35/// Max ULP 0.5
36#[inline]
37pub fn f_sincospif(x: f32) -> (f32, f32) {
38 let x_abs = x.to_bits() & 0x7fff_ffffu32;
39 let xd = x as f64;
40
41 // |x| <= 1/16
42 if x_abs <= 0x3d80_0000u32 {
43 // |x| < 0.00000009546391
44 if x_abs < 0x38a2_f984u32 {
45 const PI: f64 = f64::from_bits(0x400921fb54442d18);
46 const MPI_E3_OVER_6: f64 = f64::from_bits(0xc014abbce625be53);
47
48 // Small values approximated with Taylor poly for sine
49 // x = pi * x - pi^3*x^3/6
50 let x2 = xd * xd;
51 let p = f_fmla(x2, MPI_E3_OVER_6, PI);
52 let sf = (xd * p) as f32;
53 #[cfg(any(
54 all(
55 any(target_arch = "x86", target_arch = "x86_64"),
56 target_feature = "fma"
57 ),
58 all(target_arch = "aarch64", target_feature = "neon")
59 ))]
60 {
61 use crate::common::f_fmlaf;
62 let cs = f_fmlaf(x, f32::from_bits(0xb3000000), 1.);
63 return (sf, cs);
64 }
65 #[cfg(not(any(
66 all(
67 any(target_arch = "x86", target_arch = "x86_64"),
68 target_feature = "fma"
69 ),
70 all(target_arch = "aarch64", target_feature = "neon")
71 )))]
72 {
73 let cs = f_fmla(xd, f64::from_bits(0xbe60000000000000), 1.) as f32;
74 return (sf, cs);
75 }
76 }
77
78 // Cos(x*PI)
79 // Generated poly by Sollya:
80 // d = [0, 1/16];
81 // f_cos = cos(y*pi);
82 // Q = fpminimax(f_cos, [|0, 2, 4, 6, 8|], [|D...|], d, relative, floating);
83 //
84 // See ./notes/cospif.sollya
85
86 let x2 = xd * xd;
87 let cs = f_polyeval5(
88 x2,
89 f64::from_bits(0x3ff0000000000000),
90 f64::from_bits(0xc013bd3cc9be43f7),
91 f64::from_bits(0x40103c1f08091fe0),
92 f64::from_bits(0xbff55d3ba3d94835),
93 f64::from_bits(0x3fce173c2a00e74e),
94 ) as f32;
95 /*
96 Sin(x*PI)
97 Generated by Sollya:
98 d = [0, 1/16];
99 f_sin = sin(y*pi)/y;
100 Q = fpminimax(sin(y*pi)/y, [|0, 2, 4, 6, 8|], [|D...|], d, relative, floating);
101
102 See ./notes/sinpif.sollya
103 */
104 let p = f_polyeval5(
105 x2,
106 f64::from_bits(0x400921fb54442d18),
107 f64::from_bits(0xc014abbce625bbf2),
108 f64::from_bits(0x400466bc675e116a),
109 f64::from_bits(0xbfe32d2c0b62d41c),
110 f64::from_bits(0x3fb501ec4497cb7d),
111 );
112 let sf = (xd * p) as f32;
113
114 return (sf, cs);
115 }
116
117 // Numbers greater or equal to 2^23 are always integers or NaN
118 if x_abs >= 0x4b00_0000u32 || is_integerf(x) {
119 if x_abs >= 0x7f80_0000u32 {
120 return (x + f32::NAN, x + f32::NAN);
121 }
122 static SF: [f32; 2] = [0., -0.];
123 let sf = SF[x.is_sign_negative() as usize];
124 if x_abs < 0x4b80_0000u32 {
125 static CF: [f32; 2] = [1., -1.];
126 let cs = CF[is_odd_integerf(x) as usize];
127 return (sf, cs);
128 }
129 return (sf, 1.);
130 }
131
132 // Formula:
133 // cos(x) = cos((k + y)*pi/32)
134 // = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
135 // The values of sin(k*pi/32) and cos(k*pi/32) for k = 0..63 are precomputed
136 // and stored using a vector of 32 doubles. Sin(y*pi/32) and cos(y*pi/32) are
137 // computed using degree-7 and degree-6 minimax polynomials generated by
138 // Sollya respectively.
139 // Combine the results with the sine of sum formula:
140 // cos(x) = cos((k + y)*pi/32)
141 // = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
142 // = cosm1_y * cos_k + sin_y * sin_k
143 // = (cosm1_y * cos_k + cos_k) + sin_y * sin_k
144
145 // sin(x) = sin((k + y)*pi/32)
146 // = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
147
148 let rs = sincospif_eval(xd);
149 let cs = f_fmla(rs.sin_y, -rs.sin_k, f_fmla(rs.cosm1_y, rs.cos_k, rs.cos_k)) as f32;
150 let sf = f_fmla(rs.sin_y, rs.cos_k, f_fmla(rs.cosm1_y, rs.sin_k, rs.sin_k)) as f32;
151 (sf, cs)
152}
153
154#[cfg(test)]
155mod tests {
156 use super::*;
157 use crate::{f_cospif, f_sinpif};
158
159 #[test]
160 fn test_sincospif() {
161 let v0 = f_sincospif(-5.);
162 assert_eq!(v0.0, f_sinpif(-5.));
163 assert_eq!(v0.1, f_cospif(-5.));
164
165 let v0 = f_sincospif(-4.);
166 assert_eq!(v0.0, f_sinpif(-4.));
167 assert_eq!(v0.1, f_cospif(-4.));
168
169 let v0 = f_sincospif(4.);
170 assert_eq!(v0.0, f_sinpif(4.));
171 assert_eq!(v0.1, f_cospif(4.));
172
173 let v0 = f_sincospif(-8489897.0);
174 assert_eq!(v0.0, f_sinpif(-8489897.0));
175 assert_eq!(v0.1, f_cospif(-8489897.0));
176
177 let v1 = f_sincospif(3.23);
178 assert_eq!(v1.0, f_sinpif(3.23));
179 assert_eq!(v1.1, f_cospif(3.23));
180 }
181}