Skip to main content

pxfm/sin_cosf/
sincosf_eval.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::polyeval::{f_estrin_polyeval7, f_polyeval3, f_polyeval4, f_polyeval6};
30use crate::sin_cosf::argument_reduction::ArgumentReducer;
31use crate::sin_cosf::argument_reduction_pi::ArgumentReducerPi;
32
33/**
34Generated in SageMath:
35
36```text
37print("[")
38for k in range(64):
39    k = RealField(150)(k) * RealField(150).pi() / RealField(150)(32)
40    print(double_to_hex(k.sin()) + ",")
41print("];")
42```
43**/
44static SIN_K_PI_OVER32: [u64; 64] = [
45    0x0000000000000000,
46    0x3fb917a6bc29b42c,
47    0x3fc8f8b83c69a60b,
48    0x3fd294062ed59f06,
49    0x3fd87de2a6aea963,
50    0x3fde2b5d3806f63b,
51    0x3fe1c73b39ae68c8,
52    0x3fe44cf325091dd6,
53    0x3fe6a09e667f3bcd,
54    0x3fe8bc806b151741,
55    0x3fea9b66290ea1a3,
56    0x3fec38b2f180bdb1,
57    0x3fed906bcf328d46,
58    0x3fee9f4156c62dda,
59    0x3fef6297cff75cb0,
60    0x3fefd88da3d12526,
61    0x3ff0000000000000,
62    0x3fefd88da3d12526,
63    0x3fef6297cff75cb0,
64    0x3fee9f4156c62dda,
65    0x3fed906bcf328d46,
66    0x3fec38b2f180bdb1,
67    0x3fea9b66290ea1a3,
68    0x3fe8bc806b151741,
69    0x3fe6a09e667f3bcd,
70    0x3fe44cf325091dd6,
71    0x3fe1c73b39ae68c8,
72    0x3fde2b5d3806f63b,
73    0x3fd87de2a6aea963,
74    0x3fd294062ed59f06,
75    0x3fc8f8b83c69a60b,
76    0x3fb917a6bc29b42c,
77    0xb69f77598338bfdf,
78    0xbfb917a6bc29b42c,
79    0xbfc8f8b83c69a60b,
80    0xbfd294062ed59f06,
81    0xbfd87de2a6aea963,
82    0xbfde2b5d3806f63b,
83    0xbfe1c73b39ae68c8,
84    0xbfe44cf325091dd6,
85    0xbfe6a09e667f3bcd,
86    0xbfe8bc806b151741,
87    0xbfea9b66290ea1a3,
88    0xbfec38b2f180bdb1,
89    0xbfed906bcf328d46,
90    0xbfee9f4156c62dda,
91    0xbfef6297cff75cb0,
92    0xbfefd88da3d12526,
93    0xbff0000000000000,
94    0xbfefd88da3d12526,
95    0xbfef6297cff75cb0,
96    0xbfee9f4156c62dda,
97    0xbfed906bcf328d46,
98    0xbfec38b2f180bdb1,
99    0xbfea9b66290ea1a3,
100    0xbfe8bc806b151741,
101    0xbfe6a09e667f3bcd,
102    0xbfe44cf325091dd6,
103    0xbfe1c73b39ae68c8,
104    0xbfde2b5d3806f63b,
105    0xbfd87de2a6aea963,
106    0xbfd294062ed59f06,
107    0xbfc8f8b83c69a60b,
108    0xbfb917a6bc29b42c,
109];
110
111pub(crate) struct SinCosf {
112    pub(crate) sin_k: f64,
113    pub(crate) cos_k: f64,
114    pub(crate) sin_y: f64,
115    pub(crate) cosm1_y: f64,
116}
117
118#[inline]
119pub(crate) fn sincosf_eval(x: f64, x_abs: u32) -> SinCosf {
120    let argument_reduction = ArgumentReducer { x, x_abs };
121
122    let (y, k) = argument_reduction.reduce();
123
124    let y_sqr = y * y;
125
126    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
127    // So k is an integer and -0.5 <= y <= 0.5.
128    // Then sin(x) = sin((k + y)*pi/32)
129    //             = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
130
131    let sin_k = f64::from_bits(SIN_K_PI_OVER32[((k as u64) & 63) as usize]);
132    // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).
133    // cos_k = cos(k * pi/32)
134    let cos_k = f64::from_bits(SIN_K_PI_OVER32[((k.wrapping_add(16) as u64) & 63) as usize]);
135
136    // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
137    // with:
138    // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
139    //
140    // See ./notes/sincosf.sollya
141    let sin_y = y * f_polyeval4(
142        y_sqr,
143        f64::from_bits(0x3fb921fb54442d18),
144        f64::from_bits(0xbf24abbce625abb1),
145        f64::from_bits(0x3e7466bc624f2776),
146        f64::from_bits(0xbdb32c3a619d4a7e),
147    );
148    // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
149    // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);
150    // Note that cosm1_y = cos(y*pi/32) - 1.
151    //
152    // See ./notes/sincosf.sollya
153    let cosm1_y = y_sqr
154        * f_polyeval3(
155            y_sqr,
156            f64::from_bits(0xbf73bd3cc9be430b),
157            f64::from_bits(0x3ed03c1f070c2e27),
158            f64::from_bits(0xbe155cc84bd94200),
159        );
160
161    SinCosf {
162        sin_k,
163        cos_k,
164        sin_y,
165        cosm1_y,
166    }
167}
168
169#[inline(always)]
170#[allow(unused)]
171pub(crate) fn sincosf_eval_fma(x: f64, x_abs: u32) -> SinCosf {
172    let argument_reduction = ArgumentReducer { x, x_abs };
173
174    let (y, k) = argument_reduction.reduce_fma();
175
176    let y_sqr = y * y;
177
178    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
179    // So k is an integer and -0.5 <= y <= 0.5.
180    // Then sin(x) = sin((k + y)*pi/32)
181    //             = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
182
183    let sin_k = f64::from_bits(SIN_K_PI_OVER32[((k as u64) & 63) as usize]);
184    // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).
185    // cos_k = cos(k * pi/32)
186    let cos_k = f64::from_bits(SIN_K_PI_OVER32[((k.wrapping_add(16) as u64) & 63) as usize]);
187
188    // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
189    // with:
190    // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
191    //
192    // See ./notes/sincosf.sollya
193    use crate::polyeval::d_polyeval4;
194    let sin_y = y * d_polyeval4(
195        y_sqr,
196        f64::from_bits(0x3fb921fb54442d18),
197        f64::from_bits(0xbf24abbce625abb1),
198        f64::from_bits(0x3e7466bc624f2776),
199        f64::from_bits(0xbdb32c3a619d4a7e),
200    );
201    // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
202    // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);
203    // Note that cosm1_y = cos(y*pi/32) - 1.
204    //
205    // See ./notes/sincosf.sollya
206    use crate::polyeval::d_polyeval3;
207    let cosm1_y = y_sqr
208        * d_polyeval3(
209            y_sqr,
210            f64::from_bits(0xbf73bd3cc9be430b),
211            f64::from_bits(0x3ed03c1f070c2e27),
212            f64::from_bits(0xbe155cc84bd94200),
213        );
214
215    SinCosf {
216        sin_k,
217        cos_k,
218        sin_y,
219        cosm1_y,
220    }
221}
222
223#[inline(always)]
224pub(crate) fn sincospif_eval(x: f64) -> SinCosf {
225    let argument_reduction = ArgumentReducerPi { x };
226
227    let (y, k) = argument_reduction.reduce();
228
229    let y_sqr = y * y;
230
231    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
232    // So k is an integer and -0.5 <= y <= 0.5.
233    // Then sin(x) = sin((k + y)*pi/32)
234    //             = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
235
236    let sin_k = f64::from_bits(SIN_K_PI_OVER32[((k as u64) & 63) as usize]);
237    // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).
238    // cos_k = cos(k * pi/32)
239    let cos_k = f64::from_bits(SIN_K_PI_OVER32[((k.wrapping_add(16) as u64) & 63) as usize]);
240
241    // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
242    // with:
243    // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
244    //
245    // See ./notes/sincosf.sollya
246    let sin_y = y * f_polyeval4(
247        y_sqr,
248        f64::from_bits(0x3fb921fb54442d18),
249        f64::from_bits(0xbf24abbce625abb1),
250        f64::from_bits(0x3e7466bc624f2776),
251        f64::from_bits(0xbdb32c3a619d4a7e),
252    );
253    // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
254    // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);
255    // Note that cosm1_y = cos(y*pi/32) - 1.
256    //
257    // See ./notes/sincosf.sollya
258    let cosm1_y = y_sqr
259        * f_polyeval3(
260            y_sqr,
261            f64::from_bits(0xbf73bd3cc9be430b),
262            f64::from_bits(0x3ed03c1f070c2e27),
263            f64::from_bits(0xbe155cc84bd94200),
264        );
265
266    SinCosf {
267        sin_k,
268        cos_k,
269        sin_y,
270        cosm1_y,
271    }
272}
273
274#[inline(always)]
275#[allow(unused)]
276pub(crate) fn sincospif_eval_fma(x: f64) -> SinCosf {
277    let argument_reduction = ArgumentReducerPi { x };
278
279    let (y, k) = argument_reduction.reduce_fma();
280
281    let y_sqr = y * y;
282
283    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
284    // So k is an integer and -0.5 <= y <= 0.5.
285    // Then sin(x) = sin((k + y)*pi/32)
286    //             = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
287
288    let sin_k = f64::from_bits(SIN_K_PI_OVER32[((k as u64) & 63) as usize]);
289    // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).
290    // cos_k = cos(k * pi/32)
291    let cos_k = f64::from_bits(SIN_K_PI_OVER32[((k.wrapping_add(16) as u64) & 63) as usize]);
292
293    use crate::polyeval::{d_polyeval3, d_polyeval4};
294
295    // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
296    // with:
297    // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
298    //
299    // See ./notes/sincosf.sollya
300    let sin_y = y * d_polyeval4(
301        y_sqr,
302        f64::from_bits(0x3fb921fb54442d18),
303        f64::from_bits(0xbf24abbce625abb1),
304        f64::from_bits(0x3e7466bc624f2776),
305        f64::from_bits(0xbdb32c3a619d4a7e),
306    );
307    // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
308    // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);
309    // Note that cosm1_y = cos(y*pi/32) - 1.
310    //
311    // See ./notes/sincosf.sollya
312    let cosm1_y = y_sqr
313        * d_polyeval3(
314            y_sqr,
315            f64::from_bits(0xbf73bd3cc9be430b),
316            f64::from_bits(0x3ed03c1f070c2e27),
317            f64::from_bits(0xbe155cc84bd94200),
318        );
319
320    SinCosf {
321        sin_k,
322        cos_k,
323        sin_y,
324        cosm1_y,
325    }
326}
327
328/**
329Poly for sin(pi*y) on [-0.25; 0.25].
330Generated by Sollya:
331```text
332d = [0,0.25];
333f_sin = sin(y*pi)/y;
334Q = fpminimax(f_sin, [|0, 2, 4, 6, 8, 10|], [|D...|], d, relative, floating);
335```
336See ./notes/sinpif_0p25.sollya
337**/
338#[inline(always)]
339pub(crate) fn sinpif_eval(y: f64) -> f64 {
340    let y2 = y * y;
341    f_polyeval6(
342        y2,
343        f64::from_bits(0x400921fb54442cf8),
344        f64::from_bits(0xc014abbce6257778),
345        f64::from_bits(0x400466bc670fa464),
346        f64::from_bits(0xbfe32d2c627f47da),
347        f64::from_bits(0x3fb5071be0a2d3da),
348        f64::from_bits(0xbf7dd4e0ae5b9fcd),
349    ) * y
350}
351
352/**
353Poly for sin(pi*y) on [-0.25; 0.25].
354Generated by Sollya:
355```text
356d = [0,0.25];
357f_sin = sin(y*pi)/y;
358Q = fpminimax(f_sin, [|0, 2, 4, 6, 8, 10|], [|D...|], d, relative, floating);
359```
360See ./notes/sinpif_0p25.sollya
361**/
362#[inline(always)]
363#[allow(unused)]
364pub(crate) fn sinpif_eval_fma(y: f64) -> f64 {
365    let y2 = y * y;
366    use crate::polyeval::d_polyeval6;
367    d_polyeval6(
368        y2,
369        f64::from_bits(0x400921fb54442cf8),
370        f64::from_bits(0xc014abbce6257778),
371        f64::from_bits(0x400466bc670fa464),
372        f64::from_bits(0xbfe32d2c627f47da),
373        f64::from_bits(0x3fb5071be0a2d3da),
374        f64::from_bits(0xbf7dd4e0ae5b9fcd),
375    ) * y
376}
377
378/**
379Poly for sin(pi*y) on [-0.25; 0.25].
380Generated by Sollya:
381```text
382d = [0,0.25];
383f_sin = sin(y*pi)/y;
384Q = fpminimax(f_sin, [|0, 2, 4, 6, 8, 10, 12|], [|D...|], d, relative, floating);
385```
386See ./notes/sinpif_0p25_2.sollya
387**/
388#[inline]
389pub(crate) fn sinpif_eval2(y: f64) -> f64 {
390    let y2 = y * y;
391    f_estrin_polyeval7(
392        y2,
393        f64::from_bits(0x400921fb54442d18),
394        f64::from_bits(0xc014abbce625be09),
395        f64::from_bits(0x400466bc67754fff),
396        f64::from_bits(0xbfe32d2ccdfe9424),
397        f64::from_bits(0x3fb50782d5f14825),
398        f64::from_bits(0xbf7e2fe76fdffd2b),
399        f64::from_bits(0x3f3e357ef99eb0bb),
400    ) * y
401}
402
403/**
404Poly for cos(pi*y) on [-0.25; 0.25].
405Generated by Sollya:
406```text
407d = [0, 0.25];
408f_cos = cos(y*pi);
409Q = fpminimax(f_cos, [|0, 2, 4, 6, 8, 10|], [|1, D...|], d, relative, floating);
410```
411See ./notes/cospif_0p25.sollya
412**/
413#[inline(always)]
414pub(crate) fn cospif_eval(y: f64) -> f64 {
415    let y2 = y * y;
416    f_estrin_polyeval7(
417        y2,
418        f64::from_bits(0x3ff0000000000000),
419        f64::from_bits(0xc013bd3cc9be459b),
420        f64::from_bits(0x40103c1f081b0c98),
421        f64::from_bits(0xbff55d3c7dc25308),
422        f64::from_bits(0x3fce1f4fb6e12563),
423        f64::from_bits(0xbf9a6c9c101c1182),
424        f64::from_bits(0x3f5f3d9faffda250),
425    )
426}
427
428/**
429Poly for cos(pi*y) on [-0.25; 0.25].
430Generated by Sollya:
431```text
432d = [0, 0.25];
433f_cos = cos(y*pi);
434Q = fpminimax(f_cos, [|0, 2, 4, 6, 8, 10|], [|1, D...|], d, relative, floating);
435```
436See ./notes/cospif_0p25.sollya
437**/
438#[inline(always)]
439#[allow(unused)]
440pub(crate) fn cospif_eval_fma(y: f64) -> f64 {
441    let y2 = y * y;
442    use crate::polyeval::d_estrin_polyeval7;
443    d_estrin_polyeval7(
444        y2,
445        f64::from_bits(0x3ff0000000000000),
446        f64::from_bits(0xc013bd3cc9be459b),
447        f64::from_bits(0x40103c1f081b0c98),
448        f64::from_bits(0xbff55d3c7dc25308),
449        f64::from_bits(0x3fce1f4fb6e12563),
450        f64::from_bits(0xbf9a6c9c101c1182),
451        f64::from_bits(0x3f5f3d9faffda250),
452    )
453}