pxfm/
sincospi.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 6/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::{dd_fmla, dyad_fmla, f_fmla, is_odd_integer};
30use crate::double_double::DoubleDouble;
31use crate::polyeval::{f_polyeval3, f_polyeval4};
32use crate::rounding::CpuRound;
33use crate::sin::SinCos;
34use crate::sincospi_tables::SINPI_K_PI_OVER_64;
35
36/**
37Cospi(x) on [0; 0.000244140625]
38
39Generated by Sollya:
40```text
41d = [0, 0.000244140625];
42f_cos = cos(y*pi);
43Q = fpminimax(f_cos, [|0, 2, 4, 6, 8, 10|], [|107...|], d, relative, floating);
44```
45
46See ./notes/cospi_zero_dd.sollya
47**/
48#[cold]
49#[inline(always)]
50fn as_cospi_zero<B: SinCosPiBackend>(x: f64, backend: &B) -> f64 {
51    const C: [(u64, u64); 5] = [
52        (0xbcb692b71366cc04, 0xc013bd3cc9be45de),
53        (0xbcb32b33fb803bd5, 0x40103c1f081b5ac4),
54        (0xbc9f5b752e98b088, 0xbff55d3c7e3cbff9),
55        (0x3c30023d540b9350, 0x3fce1f506446cb66),
56        (0x3c1a5d47937787d2, 0xbf8a9b062a36ba1c),
57    ];
58    let x2 = backend.exact_mult(x, x);
59    let mut p = backend.quick_mul_add(
60        x2,
61        DoubleDouble::from_bit_pair(C[3]),
62        DoubleDouble::from_bit_pair(C[3]),
63    );
64    p = backend.quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[2]));
65    p = backend.quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[1]));
66    p = backend.quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[0]));
67    p = backend.mul_add_f64(x2, p, 1.);
68    p.to_f64()
69}
70
71/**
72Sinpi on range [0.0, 0.03515625]
73
74Generated poly by Sollya:
75```text
76d = [0, 0.03515625];
77
78f_sin = sin(y*pi)/y;
79Q = fpminimax(f_sin, [|0, 2, 4, 6, 8, 10|], [|107...|], d, relative, floating);
80```
81See ./notes/sinpi_zero_dd.sollya
82**/
83#[cold]
84#[inline(always)]
85fn as_sinpi_zero<B: SinCosPiBackend>(x: f64, backend: &B) -> f64 {
86    const C: [(u64, u64); 6] = [
87        (0x3ca1a626311d9056, 0x400921fb54442d18),
88        (0x3cb055f12c462211, 0xc014abbce625be53),
89        (0xbc9789ea63534250, 0x400466bc6775aae1),
90        (0xbc78b86de6962184, 0xbfe32d2cce62874e),
91        (0x3c4eddf7cd887302, 0x3fb507833e2b781f),
92        (0x3bf180c9d4af2894, 0xbf7e2ea4e143707e),
93    ];
94    let x2 = backend.exact_mult(x, x);
95    let mut p = backend.quick_mul_add(
96        x2,
97        DoubleDouble::from_bit_pair(C[5]),
98        DoubleDouble::from_bit_pair(C[4]),
99    );
100    p = backend.quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[3]));
101    p = backend.quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[2]));
102    p = backend.quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[1]));
103    p = backend.quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[0]));
104    p = backend.quick_mult_f64(p, x);
105    p.to_f64()
106}
107
108// Return k and y, where
109// k = round(x * 64 / pi) and y = (x * 64 / pi) - k.
110#[inline]
111pub(crate) fn reduce_pi_64(x: f64) -> (f64, i64) {
112    let kd = (x * 64.).cpu_round();
113    let y = dd_fmla(kd, -1. / 64., x);
114    (y, unsafe {
115        kd.to_int_unchecked::<i64>() // indeterminate values is always filtered out before this call, as well only lowest bits are used
116    })
117}
118
119// Return k and y, where
120// k = round(x * 64 / pi) and y = (x * 64 / pi) - k.
121#[inline(always)]
122#[allow(unused)]
123pub(crate) fn reduce_pi_64_fma(x: f64) -> (f64, i64) {
124    let kd = (x * 64.).round();
125    let y = f64::mul_add(kd, -1. / 64., x);
126    (y, unsafe {
127        kd.to_int_unchecked::<i64>() // indeterminate values is always filtered out before this call, as well only lowest bits are used
128    })
129}
130
131pub(crate) trait SinCosPiBackend {
132    fn fma(&self, x: f64, y: f64, z: f64) -> f64;
133    fn dd_fma(&self, x: f64, y: f64, z: f64) -> f64;
134    fn dyad_fma(&self, x: f64, y: f64, z: f64) -> f64;
135    fn polyeval3(&self, x: f64, a0: f64, a1: f64, a2: f64) -> f64;
136    fn arg_reduce_pi_64(&self, x: f64) -> (f64, i64);
137    fn quick_mult_f64(&self, x: DoubleDouble, y: f64) -> DoubleDouble;
138    fn quick_mult(&self, x: DoubleDouble, y: DoubleDouble) -> DoubleDouble;
139    fn odd_integer(&self, x: f64) -> bool;
140    fn div(&self, x: DoubleDouble, y: DoubleDouble) -> DoubleDouble;
141    fn mul_add_f64(&self, a: DoubleDouble, b: DoubleDouble, c: f64) -> DoubleDouble;
142    fn quick_mul_add(&self, a: DoubleDouble, b: DoubleDouble, c: DoubleDouble) -> DoubleDouble;
143    fn mul_add(&self, a: DoubleDouble, b: DoubleDouble, c: DoubleDouble) -> DoubleDouble;
144    fn exact_mult(&self, x: f64, y: f64) -> DoubleDouble;
145}
146
147pub(crate) struct GenSinCosPiBackend {}
148
149impl SinCosPiBackend for GenSinCosPiBackend {
150    #[inline(always)]
151    fn fma(&self, x: f64, y: f64, z: f64) -> f64 {
152        f_fmla(x, y, z)
153    }
154    #[inline(always)]
155    fn dd_fma(&self, x: f64, y: f64, z: f64) -> f64 {
156        dd_fmla(x, y, z)
157    }
158    #[inline(always)]
159    fn dyad_fma(&self, x: f64, y: f64, z: f64) -> f64 {
160        dyad_fmla(x, y, z)
161    }
162    #[inline(always)]
163    fn polyeval3(&self, x: f64, a0: f64, a1: f64, a2: f64) -> f64 {
164        use crate::polyeval::f_polyeval3;
165        f_polyeval3(x, a0, a1, a2)
166    }
167    #[inline(always)]
168    fn arg_reduce_pi_64(&self, x: f64) -> (f64, i64) {
169        reduce_pi_64(x)
170    }
171    #[inline(always)]
172    fn quick_mult_f64(&self, x: DoubleDouble, y: f64) -> DoubleDouble {
173        DoubleDouble::quick_mult_f64(x, y)
174    }
175    #[inline(always)]
176    fn quick_mult(&self, x: DoubleDouble, y: DoubleDouble) -> DoubleDouble {
177        DoubleDouble::quick_mult(x, y)
178    }
179
180    #[inline(always)]
181    fn odd_integer(&self, x: f64) -> bool {
182        is_odd_integer(x)
183    }
184
185    #[inline(always)]
186    fn div(&self, x: DoubleDouble, y: DoubleDouble) -> DoubleDouble {
187        DoubleDouble::div(x, y)
188    }
189
190    #[inline(always)]
191    fn mul_add_f64(&self, a: DoubleDouble, b: DoubleDouble, c: f64) -> DoubleDouble {
192        DoubleDouble::mul_add_f64(a, b, c)
193    }
194
195    #[inline(always)]
196    fn quick_mul_add(&self, a: DoubleDouble, b: DoubleDouble, c: DoubleDouble) -> DoubleDouble {
197        DoubleDouble::quick_mul_add(a, b, c)
198    }
199
200    #[inline(always)]
201    fn mul_add(&self, a: DoubleDouble, b: DoubleDouble, c: DoubleDouble) -> DoubleDouble {
202        DoubleDouble::mul_add(a, b, c)
203    }
204
205    #[inline(always)]
206    fn exact_mult(&self, x: f64, y: f64) -> DoubleDouble {
207        DoubleDouble::from_exact_mult(x, y)
208    }
209}
210
211#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
212pub(crate) struct FmaSinCosPiBackend {}
213
214#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
215impl SinCosPiBackend for FmaSinCosPiBackend {
216    #[inline(always)]
217    fn fma(&self, x: f64, y: f64, z: f64) -> f64 {
218        f64::mul_add(x, y, z)
219    }
220    #[inline(always)]
221    fn dd_fma(&self, x: f64, y: f64, z: f64) -> f64 {
222        f64::mul_add(x, y, z)
223    }
224    #[inline(always)]
225    fn dyad_fma(&self, x: f64, y: f64, z: f64) -> f64 {
226        f64::mul_add(x, y, z)
227    }
228    #[inline(always)]
229    fn polyeval3(&self, x: f64, a0: f64, a1: f64, a2: f64) -> f64 {
230        use crate::polyeval::d_polyeval3;
231        d_polyeval3(x, a0, a1, a2)
232    }
233    #[inline(always)]
234    fn arg_reduce_pi_64(&self, x: f64) -> (f64, i64) {
235        reduce_pi_64_fma(x)
236    }
237    #[inline(always)]
238    fn quick_mult_f64(&self, x: DoubleDouble, y: f64) -> DoubleDouble {
239        DoubleDouble::quick_mult_f64_fma(x, y)
240    }
241    #[inline(always)]
242    fn quick_mult(&self, x: DoubleDouble, y: DoubleDouble) -> DoubleDouble {
243        DoubleDouble::quick_mult_fma(x, y)
244    }
245
246    #[inline(always)]
247    fn odd_integer(&self, x: f64) -> bool {
248        is_odd_integer(x)
249    }
250
251    #[inline(always)]
252    fn div(&self, x: DoubleDouble, y: DoubleDouble) -> DoubleDouble {
253        DoubleDouble::div_fma(x, y)
254    }
255
256    #[inline(always)]
257    fn mul_add_f64(&self, a: DoubleDouble, b: DoubleDouble, c: f64) -> DoubleDouble {
258        DoubleDouble::mul_add_f64_fma(a, b, c)
259    }
260
261    #[inline(always)]
262    fn quick_mul_add(&self, a: DoubleDouble, b: DoubleDouble, c: DoubleDouble) -> DoubleDouble {
263        DoubleDouble::quick_mul_add_fma(a, b, c)
264    }
265
266    #[inline(always)]
267    fn mul_add(&self, a: DoubleDouble, b: DoubleDouble, c: DoubleDouble) -> DoubleDouble {
268        DoubleDouble::mul_add_fma(a, b, c)
269    }
270
271    #[inline(always)]
272    fn exact_mult(&self, x: f64, y: f64) -> DoubleDouble {
273        DoubleDouble::from_exact_mult_fma(x, y)
274    }
275}
276
277#[inline(always)]
278pub(crate) fn sincospi_eval<B: SinCosPiBackend>(x: f64, backend: &B) -> SinCos {
279    let x2 = x * x;
280    /*
281        sinpi(pi*x) poly generated by Sollya:
282        d = [0, 0.0078128];
283        f_sin = sin(y*pi)/y;
284        Q = fpminimax(f_sin, [|0, 2, 4, 6|], [|107, D...|], d, relative, floating);
285        See ./notes/sinpi.sollya
286    */
287    let sin_lop = backend.polyeval3(
288        x2,
289        f64::from_bits(0xc014abbce625be4d),
290        f64::from_bits(0x400466bc6767f259),
291        f64::from_bits(0xbfe32d176b0b3baf),
292    ) * x2;
293    // We're splitting polynomial in two parts, since first term dominates
294    // we compute: (a0_lo + a0_hi) * x + x * (a1 * x^2 + a2 + x^4) ...
295    let sin_lo = backend.dd_fma(f64::from_bits(0x3ca1a5c04563817a), x, sin_lop * x);
296    let sin_hi = x * f64::from_bits(0x400921fb54442d18);
297
298    /*
299       cospi(pi*x) poly generated by Sollya:
300       d = [0, 0.015625];
301       f_cos = cos(y*pi);
302       Q = fpminimax(f_cos, [|0, 2, 4, 6, 8|], [|107, D...|], d, relative, floating);
303       See ./notes/cospi.sollya
304    */
305    let p = backend.polyeval3(
306        x2,
307        f64::from_bits(0xc013bd3cc9be45cf),
308        f64::from_bits(0x40103c1f08085ad1),
309        f64::from_bits(0xbff55d1e43463fc3),
310    );
311
312    // We're splitting polynomial in two parts, since first term dominates
313    // we compute: (a0_lo + a0_hi) + (a1 * x^2 + a2 + x^4)...
314    let cos_lo = backend.dd_fma(p, x2, f64::from_bits(0xbbdf72adefec0800));
315    let cos_hi = f64::from_bits(0x3ff0000000000000);
316
317    let err = backend.fma(
318        x2,
319        f64::from_bits(0x3cb0000000000000), // 2^-52
320        f64::from_bits(0x3c40000000000000), // 2^-59
321    );
322    SinCos {
323        v_sin: DoubleDouble::from_exact_add(sin_hi, sin_lo),
324        v_cos: DoubleDouble::from_exact_add(cos_hi, cos_lo),
325        err,
326    }
327}
328
329#[inline(always)]
330pub(crate) fn sincospi_eval_dd<B: SinCosPiBackend>(x: f64, backend: &B) -> SinCos {
331    let x2 = backend.exact_mult(x, x);
332    // Sin coeffs
333    // poly sin(pi*x) generated by Sollya:
334    // d = [0, 0.0078128];
335    // f_sin = sin(y*pi)/y;
336    // Q = fpminimax(f_sin, [|0, 2, 4, 6, 8|], [|107...|], d, relative, floating);
337    // see ./notes/sinpi_dd.sollya
338    const SC: [(u64, u64); 5] = [
339        (0x3ca1a626330ccf19, 0x400921fb54442d18),
340        (0x3cb05540f6323de9, 0xc014abbce625be53),
341        (0xbc9050fdd1229756, 0x400466bc6775aadf),
342        (0xbc780d406f3472e8, 0xbfe32d2cce5a7bf1),
343        (0x3c4cfcf8b6b817f2, 0x3fb5077069d8a182),
344    ];
345
346    let mut sin_y = backend.quick_mul_add(
347        x2,
348        DoubleDouble::from_bit_pair(SC[4]),
349        DoubleDouble::from_bit_pair(SC[3]),
350    );
351    sin_y = backend.quick_mul_add(x2, sin_y, DoubleDouble::from_bit_pair(SC[2]));
352    sin_y = backend.quick_mul_add(x2, sin_y, DoubleDouble::from_bit_pair(SC[1]));
353    sin_y = backend.quick_mul_add(x2, sin_y, DoubleDouble::from_bit_pair(SC[0]));
354    sin_y = backend.quick_mult_f64(sin_y, x);
355
356    // Cos coeffs
357    // d = [0, 0.0078128];
358    // f_cos = cos(y*pi);
359    // Q = fpminimax(f_cos, [|0, 2, 4, 6, 8|], [|107...|], d, relative, floating);
360    // See ./notes/cospi_dd.sollya
361    const CC: [(u64, u64); 5] = [
362        (0xbaaa70a580000000, 0x3ff0000000000000),
363        (0xbcb69211d8dd1237, 0xc013bd3cc9be45de),
364        (0xbcbd96cfd637eeb7, 0x40103c1f081b5abf),
365        (0x3c994d75c577f029, 0xbff55d3c7e2e4ba5),
366        (0xbc5c542d998a4e48, 0x3fce1f2f5f747411),
367    ];
368
369    let mut cos_y = backend.quick_mul_add(
370        x2,
371        DoubleDouble::from_bit_pair(CC[4]),
372        DoubleDouble::from_bit_pair(CC[3]),
373    );
374    cos_y = backend.quick_mul_add(x2, cos_y, DoubleDouble::from_bit_pair(CC[2]));
375    cos_y = backend.quick_mul_add(x2, cos_y, DoubleDouble::from_bit_pair(CC[1]));
376    cos_y = backend.quick_mul_add(x2, cos_y, DoubleDouble::from_bit_pair(CC[0]));
377    SinCos {
378        v_sin: sin_y,
379        v_cos: cos_y,
380        err: 0.,
381    }
382}
383
384#[cold]
385#[inline(always)]
386fn sinpi_dd<B: SinCosPiBackend>(
387    x: f64,
388    sin_k: DoubleDouble,
389    cos_k: DoubleDouble,
390    backend: &B,
391) -> f64 {
392    let r_sincos = sincospi_eval_dd(x, backend);
393    let cos_k_sin_y = backend.quick_mult(cos_k, r_sincos.v_sin);
394    let rr = backend.mul_add(sin_k, r_sincos.v_cos, cos_k_sin_y);
395    rr.to_f64()
396}
397
398#[cold]
399#[inline(always)]
400fn sincospi_dd<B: SinCosPiBackend>(
401    x: f64,
402    sin_sin_k: DoubleDouble,
403    sin_cos_k: DoubleDouble,
404    cos_sin_k: DoubleDouble,
405    cos_cos_k: DoubleDouble,
406    backend: &B,
407) -> (f64, f64) {
408    let r_sincos = sincospi_eval_dd(x, backend);
409
410    let cos_k_sin_y = backend.quick_mult(sin_cos_k, r_sincos.v_sin);
411    let rr_sin = backend.mul_add(sin_sin_k, r_sincos.v_cos, cos_k_sin_y);
412
413    let cos_k_sin_y = backend.quick_mult(cos_cos_k, r_sincos.v_sin);
414    let rr_cos = backend.mul_add(cos_sin_k, r_sincos.v_cos, cos_k_sin_y);
415
416    (rr_sin.to_f64(), rr_cos.to_f64())
417}
418
419// [sincospi_eval] gives precision around 2^-66 what is not enough for DD case this method gives 2^-84
420#[inline]
421fn sincospi_eval_extended(x: f64) -> SinCos {
422    let x2 = DoubleDouble::from_exact_mult(x, x);
423    /*
424        sinpi(pi*x) poly generated by Sollya:
425        d = [0, 0.0078128];
426        f_sin = sin(y*pi)/y;
427        Q = fpminimax(f_sin, [|0, 2, 4, 6, 8|], [|107, 107, D...|], d, relative, floating);
428        See ./notes/sinpi.sollya
429    */
430    let sin_lop = f_polyeval3(
431        x2.hi,
432        f64::from_bits(0x400466bc67763662),
433        f64::from_bits(0xbfe32d2cce5aad86),
434        f64::from_bits(0x3fb5077099a1f35b),
435    );
436    let mut v_sin = DoubleDouble::mul_f64_add(
437        x2,
438        sin_lop,
439        DoubleDouble::from_bit_pair((0x3cb0553d6ee5e8ec, 0xc014abbce625be53)),
440    );
441    v_sin = DoubleDouble::mul_add(
442        x2,
443        v_sin,
444        DoubleDouble::from_bit_pair((0x3ca1a626330dd130, 0x400921fb54442d18)),
445    );
446    v_sin = DoubleDouble::quick_mult_f64(v_sin, x);
447
448    /*
449       cospi(pi*x) poly generated by Sollya:
450       d = [0, 0.015625];
451       f_cos = cos(y*pi);
452       Q = fpminimax(f_cos, [|0, 2, 4, 6, 8|], [|107, 107, D...|], d, relative, floating);
453       See ./notes/cospi_fast_dd.sollya
454    */
455    let p = f_polyeval3(
456        x2.hi,
457        f64::from_bits(0x40103c1f081b5abf),
458        f64::from_bits(0xbff55d3c7e2edd89),
459        f64::from_bits(0x3fce1f2fd9d79484),
460    );
461
462    let mut v_cos = DoubleDouble::mul_f64_add(
463        x2,
464        p,
465        DoubleDouble::from_bit_pair((0xbcb69236a9b3ed73, 0xc013bd3cc9be45de)),
466    );
467    v_cos = DoubleDouble::mul_add_f64(x2, v_cos, f64::from_bits(0x3ff0000000000000));
468
469    SinCos {
470        v_sin: DoubleDouble::from_exact_add(v_sin.hi, v_sin.lo),
471        v_cos: DoubleDouble::from_exact_add(v_cos.hi, v_cos.lo),
472        err: 0.,
473    }
474}
475
476pub(crate) fn f_fast_sinpi_dd(x: f64) -> DoubleDouble {
477    let ix = x.to_bits();
478    let ax = ix & 0x7fff_ffff_ffff_ffff;
479    if ax == 0 {
480        return DoubleDouble::new(0., 0.);
481    }
482    let e: i32 = (ax >> 52) as i32;
483    let m0 = (ix & 0x000fffffffffffff) | (1u64 << 52);
484    let sgn: i64 = (ix as i64) >> 63;
485    let m = ((m0 as i64) ^ sgn).wrapping_sub(sgn);
486    let mut s: i32 = 1063i32.wrapping_sub(e);
487    if s < 0 {
488        s = -s - 1;
489        if s > 10 {
490            return DoubleDouble::new(0., f64::copysign(0.0, x));
491        }
492        let iq: u64 = (m as u64).wrapping_shl(s as u32);
493        if (iq & 2047) == 0 {
494            return DoubleDouble::new(0., f64::copysign(0.0, x));
495        }
496    }
497
498    if ax <= 0x3fa2000000000000u64 {
499        // |x| <= 0.03515625
500        const PI: DoubleDouble = DoubleDouble::new(
501            f64::from_bits(0x3ca1a62633145c07),
502            f64::from_bits(0x400921fb54442d18),
503        );
504
505        if ax < 0x3c90000000000000 {
506            // for x near zero, sinpi(x) = pi*x + O(x^3), thus worst cases are those
507            // of the function pi*x, and if x is a worst case, then 2*x is another
508            // one in the next binade. For this reason, worst cases are only included
509            // for the binade [2^-1022, 2^-1021). For larger binades,
510            // up to [2^-54,2^-53), worst cases should be deduced by multiplying
511            // by some power of 2.
512            if ax < 0x0350000000000000 {
513                let t = x * f64::from_bits(0x4690000000000000);
514                let z = DoubleDouble::quick_mult_f64(PI, t);
515                let r = z.to_f64();
516                let rs = r * f64::from_bits(0x3950000000000000);
517                let rt = rs * f64::from_bits(0x4690000000000000);
518                return DoubleDouble::new(
519                    0.,
520                    dyad_fmla((z.hi - rt) + z.lo, f64::from_bits(0x3950000000000000), rs),
521                );
522            }
523            let z = DoubleDouble::quick_mult_f64(PI, x);
524            return z;
525        }
526
527        /*
528           Poly generated by Sollya:
529           d = [0, 0.03515625];
530           f_sin = sin(y*pi)/y;
531           Q = fpminimax(f_sin, [|0, 2, 4, 6, 8, 10|], [|107, 107, D...|], d, relative, floating);
532
533           See ./notes/sinpi_zero_fast_dd.sollya
534        */
535        const C: [u64; 4] = [
536            0xbfe32d2cce62bd85,
537            0x3fb50783487eb73d,
538            0xbf7e3074f120ad1f,
539            0x3f3e8d9011340e5a,
540        ];
541
542        let x2 = DoubleDouble::from_exact_mult(x, x);
543
544        const C_PI: DoubleDouble =
545            DoubleDouble::from_bit_pair((0x3ca1a626331457a4, 0x400921fb54442d18));
546
547        let p = f_polyeval4(
548            x2.hi,
549            f64::from_bits(C[0]),
550            f64::from_bits(C[1]),
551            f64::from_bits(C[2]),
552            f64::from_bits(C[3]),
553        );
554        let mut r = DoubleDouble::mul_f64_add(
555            x2,
556            p,
557            DoubleDouble::from_bit_pair((0xbc96dd7ae221e58c, 0x400466bc6775aae2)),
558        );
559        r = DoubleDouble::mul_add(
560            x2,
561            r,
562            DoubleDouble::from_bit_pair((0x3cb05511c8a6c478, 0xc014abbce625be53)),
563        );
564        r = DoubleDouble::mul_add(r, x2, C_PI);
565        r = DoubleDouble::quick_mult_f64(r, x);
566        let k = DoubleDouble::from_exact_add(r.hi, r.lo);
567        return k;
568    }
569
570    let si = e.wrapping_sub(1011);
571    if si >= 0 && (m0.wrapping_shl(si.wrapping_add(1) as u32)) == 0 {
572        // x is integer or half-integer
573        if (m0.wrapping_shl(si as u32)) == 0 {
574            return DoubleDouble::new(0., f64::copysign(0.0, x)); // x is integer
575        }
576        let t = (m0.wrapping_shl((si - 1) as u32)) >> 63;
577        // t = 0 if |x| = 1/2 mod 2, t = 1 if |x| = 3/2 mod 2
578        return DoubleDouble::new(
579            0.,
580            if t == 0 {
581                f64::copysign(1.0, x)
582            } else {
583                -f64::copysign(1.0, x)
584            },
585        );
586    }
587
588    let (y, k) = reduce_pi_64(x);
589
590    // // cos(k * pi/64) = sin(k * pi/64 + pi/2) = sin((k + 32) * pi/64).
591    let sin_k = DoubleDouble::from_bit_pair(SINPI_K_PI_OVER_64[((k as u64) & 127) as usize]);
592    let cos_k = DoubleDouble::from_bit_pair(
593        SINPI_K_PI_OVER_64[((k as u64).wrapping_add(32) & 127) as usize],
594    );
595
596    let r_sincos = sincospi_eval_extended(y);
597
598    let sin_k_cos_y = DoubleDouble::quick_mult(sin_k, r_sincos.v_cos);
599    let cos_k_sin_y = DoubleDouble::quick_mult(cos_k, r_sincos.v_sin);
600
601    // sin_k_cos_y is always >> cos_k_sin_y
602    let mut rr = DoubleDouble::from_exact_add(sin_k_cos_y.hi, cos_k_sin_y.hi);
603    rr.lo += sin_k_cos_y.lo + cos_k_sin_y.lo;
604    DoubleDouble::from_exact_add(rr.hi, rr.lo)
605}
606
607#[inline(always)]
608fn sinpi_gen_impl<B: SinCosPiBackend>(x: f64, backend: B) -> f64 {
609    let ix = x.to_bits();
610    let ax = ix & 0x7fff_ffff_ffff_ffff;
611    if ax == 0 {
612        return x;
613    }
614    let e: i32 = (ax >> 52) as i32;
615    let m0 = (ix & 0x000fffffffffffff) | (1u64 << 52);
616    let sgn: i64 = (ix as i64) >> 63;
617    let m = ((m0 as i64) ^ sgn).wrapping_sub(sgn);
618    let mut s: i32 = 1063i32.wrapping_sub(e);
619    if s < 0 {
620        if e == 0x7ff {
621            if (ix << 12) == 0 {
622                return f64::NAN;
623            }
624            return x + x; // case x=NaN
625        }
626        s = -s - 1;
627        if s > 10 {
628            return f64::copysign(0.0, x);
629        }
630        let iq: u64 = (m as u64).wrapping_shl(s as u32);
631        if (iq & 2047) == 0 {
632            return f64::copysign(0.0, x);
633        }
634    }
635
636    if ax <= 0x3fa2000000000000u64 {
637        // |x| <= 0.03515625
638        const PI: DoubleDouble = DoubleDouble::new(
639            f64::from_bits(0x3ca1a62633145c07),
640            f64::from_bits(0x400921fb54442d18),
641        );
642
643        if ax < 0x3c90000000000000 {
644            // for x near zero, sinpi(x) = pi*x + O(x^3), thus worst cases are those
645            // of the function pi*x, and if x is a worst case, then 2*x is another
646            // one in the next binade. For this reason, worst cases are only included
647            // for the binade [2^-1022, 2^-1021). For larger binades,
648            // up to [2^-54,2^-53), worst cases should be deduced by multiplying
649            // by some power of 2.
650            if ax < 0x0350000000000000 {
651                let t = x * f64::from_bits(0x4690000000000000);
652                let z = backend.quick_mult_f64(PI, t);
653                let r = z.to_f64();
654                let rs = r * f64::from_bits(0x3950000000000000);
655                let rt = rs * f64::from_bits(0x4690000000000000);
656                return backend.dyad_fma(
657                    (z.hi - rt) + z.lo,
658                    f64::from_bits(0x3950000000000000),
659                    rs,
660                );
661            }
662            let z = backend.quick_mult_f64(PI, x);
663            return z.to_f64();
664        }
665
666        /*
667           Poly generated by Sollya:
668           d = [0, 0.03515625];
669           f_sin = sin(y*pi)/y;
670           Q = fpminimax(f_sin, [|0, 2, 4, 6, 8, 10|], [|107, D...|], d, relative, floating);
671
672           See ./notes/sinpi_zero.sollya
673        */
674
675        let x2 = x * x;
676        let x3 = x2 * x;
677        let x4 = x2 * x2;
678
679        let eps = x * backend.fma(
680            x2,
681            f64::from_bits(0x3d00000000000000), // 2^-47
682            f64::from_bits(0x3bd0000000000000), // 2^-66
683        );
684
685        const C: [u64; 4] = [
686            0xc014abbce625be51,
687            0x400466bc67754b46,
688            0xbfe32d2cc12a51f4,
689            0x3fb5060540058476,
690        ];
691
692        const C_PI: DoubleDouble =
693            DoubleDouble::from_bit_pair((0x3ca1a67088eb1a46, 0x400921fb54442d18));
694
695        let mut z = backend.quick_mult_f64(C_PI, x);
696
697        let zl0 = backend.fma(x2, f64::from_bits(C[1]), f64::from_bits(C[0]));
698        let zl1 = backend.fma(x2, f64::from_bits(C[3]), f64::from_bits(C[2]));
699
700        z.lo = backend.fma(x3, backend.fma(x4, zl1, zl0), z.lo);
701        let lb = z.hi + (z.lo - eps);
702        let ub = z.hi + (z.lo + eps);
703        if lb == ub {
704            return lb;
705        }
706        return as_sinpi_zero(x, &backend);
707    }
708
709    let si = e.wrapping_sub(1011);
710    if si >= 0 && (m0.wrapping_shl(si.wrapping_add(1) as u32)) == 0 {
711        // x is integer or half-integer
712        if (m0.wrapping_shl(si as u32)) == 0 {
713            return f64::copysign(0.0, x); // x is integer
714        }
715        let t = (m0.wrapping_shl((si - 1) as u32)) >> 63;
716        // t = 0 if |x| = 1/2 mod 2, t = 1 if |x| = 3/2 mod 2
717        return if t == 0 {
718            f64::copysign(1.0, x)
719        } else {
720            -f64::copysign(1.0, x)
721        };
722    }
723
724    let (y, k) = backend.arg_reduce_pi_64(x);
725
726    // cos(k * pi/64) = sin(k * pi/64 + pi/2) = sin((k + 32) * pi/64).
727    let sin_k = DoubleDouble::from_bit_pair(SINPI_K_PI_OVER_64[((k as u64) & 127) as usize]);
728    let cos_k = DoubleDouble::from_bit_pair(
729        SINPI_K_PI_OVER_64[((k as u64).wrapping_add(32) & 127) as usize],
730    );
731
732    let r_sincos = sincospi_eval(y, &backend);
733
734    let sin_k_cos_y = backend.quick_mult(sin_k, r_sincos.v_cos);
735    let cos_k_sin_y = backend.quick_mult(cos_k, r_sincos.v_sin);
736
737    // sin_k_cos_y is always >> cos_k_sin_y
738    let mut rr = DoubleDouble::from_exact_add(sin_k_cos_y.hi, cos_k_sin_y.hi);
739    rr.lo += sin_k_cos_y.lo + cos_k_sin_y.lo;
740
741    let ub = rr.hi + (rr.lo + r_sincos.err); // (rr.lo + ERR);
742    let lb = rr.hi + (rr.lo - r_sincos.err); // (rr.lo - ERR);
743
744    if ub == lb {
745        return rr.to_f64();
746    }
747    sinpi_dd(y, sin_k, cos_k, &backend)
748}
749
750#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
751#[target_feature(enable = "avx", enable = "fma")]
752unsafe fn sinpi_fma_impl(x: f64) -> f64 {
753    sinpi_gen_impl(x, FmaSinCosPiBackend {})
754}
755
756/// Computes sin(PI*x)
757///
758/// Max ULP 0.5
759pub fn f_sinpi(x: f64) -> f64 {
760    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
761    {
762        sinpi_gen_impl(x, GenSinCosPiBackend {})
763    }
764    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
765    {
766        use std::sync::OnceLock;
767        static EXECUTOR: OnceLock<unsafe fn(f64) -> f64> = OnceLock::new();
768        let q = EXECUTOR.get_or_init(|| {
769            if std::arch::is_x86_feature_detected!("avx")
770                && std::arch::is_x86_feature_detected!("fma")
771            {
772                sinpi_fma_impl
773            } else {
774                fn def_sinpi(x: f64) -> f64 {
775                    sinpi_gen_impl(x, GenSinCosPiBackend {})
776                }
777                def_sinpi
778            }
779        });
780        unsafe { q(x) }
781    }
782}
783
784#[inline(always)]
785fn cospi_gen_impl<B: SinCosPiBackend>(x: f64, backend: B) -> f64 {
786    let ix = x.to_bits();
787    let ax = ix & 0x7fff_ffff_ffff_ffff;
788    if ax == 0 {
789        return 1.0;
790    }
791    let e: i32 = (ax >> 52) as i32;
792    // e is the unbiased exponent, we have 2^(e-1023) <= |x| < 2^(e-1022)
793    let m: i64 = ((ix & 0x000fffffffffffff) | (1u64 << 52)) as i64;
794    let mut s = 1063i32.wrapping_sub(e); // 2^(40-s) <= |x| < 2^(41-s)
795    if s < 0 {
796        // |x| >= 2^41
797        if e == 0x7ff {
798            // NaN or Inf
799            if ix.wrapping_shl(12) == 0 {
800                return f64::NAN;
801            }
802            return x + x; // NaN
803        }
804        s = -s - 1; // now 2^(41+s) <= |x| < 2^(42+s)
805        if s > 11 {
806            return 1.0;
807        } // |x| >= 2^53
808        let iq: u64 = (m as u64).wrapping_shl(s as u32).wrapping_add(1024);
809        if (iq & 2047) == 0 {
810            return 0.0;
811        }
812    }
813    if ax <= 0x3f30000000000000u64 {
814        // |x| <= 2^-12, |x| <= 0.000244140625
815        if ax <= 0x3e2ccf6429be6621u64 {
816            return 1.0 - f64::from_bits(0x3c80000000000000);
817        }
818        let x2 = x * x;
819        let x4 = x2 * x2;
820        let eps = x2 * f64::from_bits(0x3cfa000000000000);
821
822        /*
823            Generated by Sollya:
824            d = [0, 0.000244140625];
825            f_cos = cos(y*pi);
826            Q = fpminimax(f_cos, [|0, 2, 4, 6, 8|], [|107, 107, D...|], d, relative, floating);
827
828            See ./notes/cospi.sollya
829        */
830
831        const C: [u64; 4] = [
832            0xc013bd3cc9be45de,
833            0x40103c1f081b5ac4,
834            0xbff55d3c7ff79b60,
835            0x3fd24c7b6f7d0690,
836        ];
837
838        let p0 = backend.fma(x2, f64::from_bits(C[3]), f64::from_bits(C[2]));
839        let p1 = backend.fma(x2, f64::from_bits(C[1]), f64::from_bits(C[0]));
840
841        let p = x2 * backend.fma(x4, p0, p1);
842        let lb = (p - eps) + 1.;
843        let ub = (p + eps) + 1.;
844        if lb == ub {
845            return lb;
846        }
847        return as_cospi_zero(x, &backend);
848    }
849
850    let si: i32 = e.wrapping_sub(1011);
851    if si >= 0 && ((m as u64).wrapping_shl(si as u32) ^ 0x8000000000000000u64) == 0 {
852        return 0.0;
853    }
854
855    let (y, k) = backend.arg_reduce_pi_64(x);
856
857    // cos(k * pi/64) = sin(k * pi/64 + pi/2) = sin((k + 32) * pi/64).
858    let msin_k = DoubleDouble::from_bit_pair(
859        SINPI_K_PI_OVER_64[((k as u64).wrapping_add(64) & 127) as usize],
860    );
861    let cos_k = DoubleDouble::from_bit_pair(
862        SINPI_K_PI_OVER_64[((k as u64).wrapping_add(32) & 127) as usize],
863    );
864
865    let r_sincos = sincospi_eval(y, &backend);
866
867    let cos_k_cos_y = backend.quick_mult(r_sincos.v_cos, cos_k);
868    let cos_k_msin_y = backend.quick_mult(r_sincos.v_sin, msin_k);
869
870    // cos_k_cos_y is always >> cos_k_msin_y
871    let mut rr = DoubleDouble::from_exact_add(cos_k_cos_y.hi, cos_k_msin_y.hi);
872    rr.lo += cos_k_cos_y.lo + cos_k_msin_y.lo;
873
874    let ub = rr.hi + (rr.lo + r_sincos.err); // (rr.lo + ERR);
875    let lb = rr.hi + (rr.lo - r_sincos.err); // (rr.lo - ERR);
876
877    if ub == lb {
878        return rr.to_f64();
879    }
880    sinpi_dd(y, cos_k, msin_k, &backend)
881}
882
883#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
884#[target_feature(enable = "avx", enable = "fma")]
885unsafe fn cospi_fma_impl(x: f64) -> f64 {
886    cospi_gen_impl(x, FmaSinCosPiBackend {})
887}
888
889/// Computes cos(PI*x)
890///
891/// Max found ULP 0.5
892pub fn f_cospi(x: f64) -> f64 {
893    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
894    {
895        cospi_gen_impl(x, GenSinCosPiBackend {})
896    }
897    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
898    {
899        use std::sync::OnceLock;
900        static EXECUTOR: OnceLock<unsafe fn(f64) -> f64> = OnceLock::new();
901        let q = EXECUTOR.get_or_init(|| {
902            if std::arch::is_x86_feature_detected!("avx")
903                && std::arch::is_x86_feature_detected!("fma")
904            {
905                cospi_fma_impl
906            } else {
907                fn def_cospi(x: f64) -> f64 {
908                    cospi_gen_impl(x, GenSinCosPiBackend {})
909                }
910                def_cospi
911            }
912        });
913        unsafe { q(x) }
914    }
915}
916
917#[inline(always)]
918fn sincospi_gen_impl<B: SinCosPiBackend>(x: f64, backend: B) -> (f64, f64) {
919    let ix = x.to_bits();
920    let ax = ix & 0x7fff_ffff_ffff_ffff;
921    if ax == 0 {
922        return (x, 1.0);
923    }
924    let e: i32 = (ax >> 52) as i32;
925    // e is the unbiased exponent, we have 2^(e-1023) <= |x| < 2^(e-1022)
926    let m0 = (ix & 0x000fffffffffffff) | (1u64 << 52);
927    let m: i64 = ((ix & 0x000fffffffffffff) | (1u64 << 52)) as i64;
928    let mut s = 1063i32.wrapping_sub(e); // 2^(40-s) <= |x| < 2^(41-s)
929    if s < 0 {
930        // |x| >= 2^41
931        if e == 0x7ff {
932            // NaN or Inf
933            if ix.wrapping_shl(12) == 0 {
934                return (f64::NAN, f64::NAN);
935            }
936            return (x + x, x + x); // NaN
937        }
938        s = -s - 1;
939        if s > 10 {
940            static CF: [f64; 2] = [1., -1.];
941            let is_odd = backend.odd_integer(f64::from_bits(ax));
942            let cos_x = CF[is_odd as usize];
943            return (f64::copysign(0.0, x), cos_x);
944        } // |x| >= 2^53
945        let iq: u64 = (m as u64).wrapping_shl(s as u32);
946
947        // sinpi = 0 when multiple of 2048
948        let sin_zero = (iq & 2047) == 0;
949
950        // cospi = 0 when offset-by-half multiple of 2048
951        let cos_zero = ((m as u64).wrapping_shl(s as u32).wrapping_add(1024) & 2047) == 0;
952
953        if sin_zero && cos_zero {
954            // both zero (only possible if NaN or something degenerate)
955        } else if sin_zero {
956            static CF: [f64; 2] = [1., -1.];
957            let is_odd = backend.odd_integer(f64::from_bits(ax));
958            let cos_x = CF[is_odd as usize];
959            return (0.0, cos_x); // sin = 0, cos = ±1
960        } else if cos_zero {
961            // x = k / 2 * PI
962            let si = e.wrapping_sub(1011);
963            let t = (m0.wrapping_shl((si - 1) as u32)) >> 63;
964            // making sin decision based on quadrant
965            return if t == 0 {
966                (f64::copysign(1.0, x), 0.0)
967            } else {
968                (-f64::copysign(1.0, x), 0.0)
969            }; // sin = ±1, cos = 0
970        }
971    }
972
973    if ax <= 0x3f30000000000000u64 {
974        // |x| <= 2^-12, |x| <= 0.000244140625
975        if ax <= 0x3c90000000000000u64 {
976            const PI: DoubleDouble = DoubleDouble::new(
977                f64::from_bits(0x3ca1a62633145c07),
978                f64::from_bits(0x400921fb54442d18),
979            );
980            let sin_x = if ax < 0x0350000000000000 {
981                let t = x * f64::from_bits(0x4690000000000000);
982                let z = backend.quick_mult_f64(PI, t);
983                let r = z.to_f64();
984                let rs = r * f64::from_bits(0x3950000000000000);
985                let rt = rs * f64::from_bits(0x4690000000000000);
986                backend.dyad_fma((z.hi - rt) + z.lo, f64::from_bits(0x3950000000000000), rs)
987            } else {
988                let z = backend.quick_mult_f64(PI, x);
989                z.to_f64()
990            };
991            return (sin_x, 1.0 - f64::from_bits(0x3c80000000000000));
992        }
993        let x2 = x * x;
994        let x4 = x2 * x2;
995        let cos_eps = x2 * f64::from_bits(0x3cfa000000000000);
996
997        /*
998            Generated by Sollya:
999            d = [0, 0.000244140625];
1000            f_cos = cos(y*pi);
1001            Q = fpminimax(f_cos, [|0, 2, 4, 6, 8|], [|107, 107, D...|], d, relative, floating);
1002
1003            See ./notes/cospi.sollya
1004        */
1005
1006        const COS_C: [u64; 4] = [
1007            0xc013bd3cc9be45de,
1008            0x40103c1f081b5ac4,
1009            0xbff55d3c7ff79b60,
1010            0x3fd24c7b6f7d0690,
1011        ];
1012
1013        let p0 = backend.fma(x2, f64::from_bits(COS_C[3]), f64::from_bits(COS_C[2]));
1014        let p1 = backend.fma(x2, f64::from_bits(COS_C[1]), f64::from_bits(COS_C[0]));
1015
1016        let p = x2 * backend.fma(x4, p0, p1);
1017        let cos_lb = (p - cos_eps) + 1.;
1018        let cos_ub = (p + cos_eps) + 1.;
1019        let cos_x = if cos_lb == cos_ub {
1020            cos_lb
1021        } else {
1022            as_cospi_zero(x, &backend)
1023        };
1024
1025        /*
1026            Poly generated by Sollya:
1027            d = [0, 0.03515625];
1028            f_sin = sin(y*pi)/y;
1029            Q = fpminimax(f_sin, [|0, 2, 4, 6, 8, 10|], [|107, D...|], d, relative, floating);
1030
1031            See ./notes/sinpi_zero.sollya
1032        */
1033
1034        const SIN_C: [u64; 4] = [
1035            0xc014abbce625be51,
1036            0x400466bc67754b46,
1037            0xbfe32d2cc12a51f4,
1038            0x3fb5060540058476,
1039        ];
1040
1041        const C_PI: DoubleDouble =
1042            DoubleDouble::from_bit_pair((0x3ca1a67088eb1a46, 0x400921fb54442d18));
1043
1044        let mut z = backend.quick_mult_f64(C_PI, x);
1045
1046        let x3 = x2 * x;
1047
1048        let zl0 = backend.fma(x2, f64::from_bits(SIN_C[1]), f64::from_bits(SIN_C[0]));
1049        let zl1 = backend.fma(x2, f64::from_bits(SIN_C[3]), f64::from_bits(SIN_C[2]));
1050
1051        let sin_eps = x * backend.fma(
1052            x2,
1053            f64::from_bits(0x3d00000000000000), // 2^-47
1054            f64::from_bits(0x3bd0000000000000), // 2^-66
1055        );
1056
1057        z.lo = backend.fma(x3, backend.fma(x4, zl1, zl0), z.lo);
1058        let sin_lb = z.hi + (z.lo - sin_eps);
1059        let sin_ub = z.hi + (z.lo + sin_eps);
1060        let sin_x = if sin_lb == sin_ub {
1061            sin_lb
1062        } else {
1063            as_sinpi_zero(x, &backend)
1064        };
1065        return (sin_x, cos_x);
1066    }
1067
1068    let si = e.wrapping_sub(1011);
1069    if si >= 0 && (m0.wrapping_shl(si.wrapping_add(1) as u32)) == 0 {
1070        // x is integer or half-integer
1071        if (m0.wrapping_shl(si as u32)) == 0 {
1072            static CF: [f64; 2] = [1., -1.];
1073            let is_odd = backend.odd_integer(f64::from_bits(ax));
1074            let cos_x = CF[is_odd as usize];
1075            return (f64::copysign(0.0, x), cos_x); // x is integer
1076        }
1077        // x is half-integer
1078        let t = (m0.wrapping_shl((si - 1) as u32)) >> 63;
1079        // t = 0 if |x| = 1/2 mod 2, t = 1 if |x| = 3/2 mod 2
1080        return if t == 0 {
1081            (f64::copysign(1.0, x), 0.0)
1082        } else {
1083            (-f64::copysign(1.0, x), 0.0)
1084        };
1085    }
1086
1087    let (y, k) = backend.arg_reduce_pi_64(x);
1088
1089    // cos(k * pi/64) = sin(k * pi/64 + pi/2) = sin((k + 32) * pi/64).
1090    let sin_k = DoubleDouble::from_bit_pair(SINPI_K_PI_OVER_64[((k as u64) & 127) as usize]);
1091    let cos_k = DoubleDouble::from_bit_pair(
1092        SINPI_K_PI_OVER_64[((k as u64).wrapping_add(32) & 127) as usize],
1093    );
1094    let msin_k = -sin_k;
1095
1096    let r_sincos = sincospi_eval(y, &backend);
1097
1098    let sin_k_cos_y = backend.quick_mult(sin_k, r_sincos.v_cos);
1099    let cos_k_sin_y = backend.quick_mult(cos_k, r_sincos.v_sin);
1100
1101    let cos_k_cos_y = backend.quick_mult(r_sincos.v_cos, cos_k);
1102    let msin_k_sin_y = backend.quick_mult(r_sincos.v_sin, msin_k);
1103
1104    // sin_k_cos_y is always >> cos_k_sin_y
1105    let mut rr_sin = DoubleDouble::from_exact_add(sin_k_cos_y.hi, cos_k_sin_y.hi);
1106    rr_sin.lo += sin_k_cos_y.lo + cos_k_sin_y.lo;
1107
1108    let sin_ub = rr_sin.hi + (rr_sin.lo + r_sincos.err); // (rr.lo + ERR);
1109    let sin_lb = rr_sin.hi + (rr_sin.lo - r_sincos.err); // (rr.lo - ERR);
1110
1111    let mut rr_cos = DoubleDouble::from_exact_add(cos_k_cos_y.hi, msin_k_sin_y.hi);
1112    rr_cos.lo += cos_k_cos_y.lo + msin_k_sin_y.lo;
1113
1114    let cos_ub = rr_cos.hi + (rr_cos.lo + r_sincos.err); // (rr.lo + ERR);
1115    let cos_lb = rr_cos.hi + (rr_cos.lo - r_sincos.err); // (rr.lo - ERR);
1116
1117    if sin_ub == sin_lb && cos_lb == cos_ub {
1118        return (rr_sin.to_f64(), rr_cos.to_f64());
1119    }
1120
1121    sincospi_dd(y, sin_k, cos_k, cos_k, msin_k, &backend)
1122}
1123
1124#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
1125#[target_feature(enable = "avx", enable = "fma")]
1126unsafe fn sincospi_fma_impl(x: f64) -> (f64, f64) {
1127    sincospi_gen_impl(x, FmaSinCosPiBackend {})
1128}
1129
1130/// Computes sin(PI*x) and cos(PI*x)
1131///
1132/// Max found ULP 0.5
1133pub fn f_sincospi(x: f64) -> (f64, f64) {
1134    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
1135    {
1136        sincospi_gen_impl(x, GenSinCosPiBackend {})
1137    }
1138    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
1139    {
1140        use std::sync::OnceLock;
1141        static EXECUTOR: OnceLock<unsafe fn(f64) -> (f64, f64)> = OnceLock::new();
1142        let q = EXECUTOR.get_or_init(|| {
1143            if std::arch::is_x86_feature_detected!("avx")
1144                && std::arch::is_x86_feature_detected!("fma")
1145            {
1146                sincospi_fma_impl
1147            } else {
1148                fn def_sincospi(x: f64) -> (f64, f64) {
1149                    sincospi_gen_impl(x, GenSinCosPiBackend {})
1150                }
1151                def_sincospi
1152            }
1153        });
1154        unsafe { q(x) }
1155    }
1156}
1157
1158#[cfg(test)]
1159mod tests {
1160    use super::*;
1161
1162    #[test]
1163    fn test_sinpi() {
1164        assert_eq!(f_sinpi(262143.50006870925), -0.9999999767029883);
1165        assert_eq!(f_sinpi(7124076477593855.), 0.);
1166        assert_eq!(f_sinpi(-11235582092889474000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.), -0.);
1167        assert_eq!(f_sinpi(-2.7430620343968443e303), -0.0);
1168        assert_eq!(f_sinpi(0.00003195557007273919), 0.00010039138401316004);
1169        assert_eq!(f_sinpi(-0.038357843137253766), -0.12021328061499763);
1170        assert_eq!(f_sinpi(1.0156097449358867), -0.04901980680173724);
1171        assert_eq!(f_sinpi(74.8593852519989), 0.42752597787896457);
1172        assert_eq!(f_sinpi(0.500091552734375), 0.9999999586369661);
1173        assert_eq!(f_sinpi(0.5307886532952182), 0.9953257438106751);
1174        assert_eq!(f_sinpi(3.1415926535897936), -0.43030121700009316);
1175        assert_eq!(f_sinpi(-0.5305172747685276), -0.9954077178320563);
1176        assert_eq!(f_sinpi(-0.03723630312089732), -0.1167146713267927);
1177        assert_eq!(
1178            f_sinpi(0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000022946074000077123),
1179            0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000007208721750737005
1180        );
1181        assert_eq!(
1182            f_sinpi(0.000000000000000000000000000000000000007413093439574428),
1183            2.3288919890141717e-38
1184        );
1185        assert_eq!(f_sinpi(0.0031909299901270445), 0.0100244343161398578);
1186        assert_eq!(f_sinpi(0.11909245901270445), 0.36547215190661003);
1187        assert_eq!(f_sinpi(0.99909245901270445), 0.0028511202357662186);
1188        assert!(f_sinpi(f64::INFINITY).is_nan());
1189        assert!(f_sinpi(f64::NEG_INFINITY).is_nan());
1190        assert!(f_sinpi(f64::NAN).is_nan());
1191    }
1192
1193    #[test]
1194    fn test_sincospi() {
1195        let v0 = f_sincospi(1.0156097449358867);
1196        assert_eq!(v0.0, f_sinpi(1.0156097449358867));
1197        assert_eq!(v0.1, f_cospi(1.0156097449358867));
1198
1199        let v1 = f_sincospi(4503599627370496.);
1200        assert_eq!(v1.0, f_sinpi(4503599627370496.));
1201        assert_eq!(v1.1, f_cospi(4503599627370496.));
1202
1203        let v1 = f_sincospi(-108.);
1204        assert_eq!(v1.0, f_sinpi(-108.));
1205        assert_eq!(v1.1, f_cospi(-108.));
1206
1207        let v1 = f_sincospi(3.);
1208        assert_eq!(v1.0, f_sinpi(3.));
1209        assert_eq!(v1.1, f_cospi(3.));
1210
1211        let v1 = f_sincospi(13.5);
1212        assert_eq!(v1.0, f_sinpi(13.5));
1213        assert_eq!(v1.1, f_cospi(13.5));
1214
1215        let v1 = f_sincospi(7124076477593855.);
1216        assert_eq!(v1.0, f_sinpi(7124076477593855.));
1217        assert_eq!(v1.1, f_cospi(7124076477593855.));
1218
1219        let v1 = f_sincospi(2533419148247186.5);
1220        assert_eq!(v1.0, f_sinpi(2533419148247186.5));
1221        assert_eq!(v1.1, f_cospi(2533419148247186.5));
1222
1223        let v1 = f_sincospi(2.2250653705240375E-308);
1224        assert_eq!(v1.0, f_sinpi(2.2250653705240375E-308));
1225        assert_eq!(v1.1, f_cospi(2.2250653705240375E-308));
1226
1227        let v1 = f_sincospi(2533420818956351.);
1228        assert_eq!(v1.0, f_sinpi(2533420818956351.));
1229        assert_eq!(v1.1, f_cospi(2533420818956351.));
1230
1231        let v1 = f_sincospi(2533822406803233.5);
1232        assert_eq!(v1.0, f_sinpi(2533822406803233.5));
1233        assert_eq!(v1.1, f_cospi(2533822406803233.5));
1234
1235        let v1 = f_sincospi(-3040685725640478.5);
1236        assert_eq!(v1.0, f_sinpi(-3040685725640478.5));
1237        assert_eq!(v1.1, f_cospi(-3040685725640478.5));
1238
1239        let v1 = f_sincospi(2533419148247186.5);
1240        assert_eq!(v1.0, f_sinpi(2533419148247186.5));
1241        assert_eq!(v1.1, f_cospi(2533419148247186.5));
1242
1243        let v1 = f_sincospi(2533420819267583.5);
1244        assert_eq!(v1.0, f_sinpi(2533420819267583.5));
1245        assert_eq!(v1.1, f_cospi(2533420819267583.5));
1246
1247        let v1 = f_sincospi(6979704728846336.);
1248        assert_eq!(v1.0, f_sinpi(6979704728846336.));
1249        assert_eq!(v1.1, f_cospi(6979704728846336.));
1250
1251        let v1 = f_sincospi(7124076477593855.);
1252        assert_eq!(v1.0, f_sinpi(7124076477593855.));
1253        assert_eq!(v1.1, f_cospi(7124076477593855.));
1254
1255        let v1 = f_sincospi(-0.00000000002728839192371484);
1256        assert_eq!(v1.0, f_sinpi(-0.00000000002728839192371484));
1257        assert_eq!(v1.1, f_cospi(-0.00000000002728839192371484));
1258
1259        let v1 = f_sincospi(0.00002465398569495569);
1260        assert_eq!(v1.0, f_sinpi(0.00002465398569495569));
1261        assert_eq!(v1.1, f_cospi(0.00002465398569495569));
1262    }
1263
1264    #[test]
1265    fn test_cospi() {
1266        assert_eq!(0.9999497540959953, f_cospi(0.0031909299901270445));
1267        assert_eq!(0.9308216542079669, f_cospi(0.11909299901270445));
1268        assert_eq!(-0.1536194873288318, f_cospi(0.54909299901270445));
1269        assert!(f_cospi(f64::INFINITY).is_nan());
1270        assert!(f_cospi(f64::NEG_INFINITY).is_nan());
1271        assert!(f_cospi(f64::NAN).is_nan());
1272    }
1273}
pxfm/sincospi.rs

pxfm/
sincospi.rs
