pxfm/
sincos.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::*;
30use crate::double_double::DoubleDouble;
31use crate::sin::{
32    cos_accurate_near_zero, range_reduction_small, sin_accurate_near_zero, sincos_eval,
33};
34use crate::sin_helper::sincos_eval_dd;
35use crate::sin_table::SIN_K_PI_OVER_128;
36use crate::sincos_reduce::LargeArgumentReduction;
37use std::hint::black_box;
38
39/// Sine and cosine for double precision
40///
41/// ULP 0.5
42pub fn f_sincos(x: f64) -> (f64, f64) {
43    let x_e = (x.to_bits() >> 52) & 0x7ff;
44    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
45
46    let y: DoubleDouble;
47    let k;
48
49    let mut argument_reduction = LargeArgumentReduction::default();
50
51    // |x| < 2^32 (with FMA) or |x| < 2^23 (w/o FMA)
52    if x_e < E_BIAS + 16 {
53        // |x| < 2^-26
54        let ax = x.to_bits() & 0x7fff_ffff_ffff_ffff;
55        if ax <= 0x3fa921fbd34a9550 {
56            // |x| <= 0.0490874
57            if x_e < E_BIAS - 27 {
58                // Signed zeros.
59                if x == 0.0 {
60                    return (x, 1.0);
61                }
62                // For |x| < 2^-26, |sin(x) - x| < ulp(x)/2.
63                let s_sin = dyad_fmla(x, f64::from_bits(0xbc90000000000000), x);
64                let s_cos = black_box(1.0) - min_normal_f64();
65                return (s_sin, s_cos);
66            }
67
68            // Polynomial for sin(x)/x
69            // Generated by Sollya:
70            // d = [0, 0.0490874];
71            // f_sin = sin(y)/y;
72            // Q = fpminimax(f_sin, [|0, 2, 4, 6, 8|], [|1, D...|], d, relative, floating);
73            const SIN_C: [u64; 4] = [
74                0xbfc5555555555555,
75                0x3f8111111110e45a,
76                0xbf2a019ffd7fdaaf,
77                0x3ec71819b9bf01ef,
78            ];
79
80            let x2 = x * x;
81            let x4 = x2 * x2;
82
83            let p01 = f_fmla(x2, f64::from_bits(SIN_C[1]), f64::from_bits(SIN_C[0]));
84            let p23 = f_fmla(x2, f64::from_bits(SIN_C[3]), f64::from_bits(SIN_C[2]));
85            let w0 = f_fmla(x4, p23, p01);
86            let w1 = x2 * w0 * x;
87            let r_sin = DoubleDouble::from_exact_add(x, w1);
88
89            // Polynomial for cos(x)
90            // Generated by Sollya:
91            // d = [0, 0.0490874];
92            // f_cos = cos(y);
93            // Q = fpminimax(f_cos, [|0, 2, 4, 6, 8|], [|1, D...|], d, relative, floating);
94            const COS_C: [u64; 4] = [
95                0xbfe0000000000000,
96                0x3fa55555555554a4,
97                0xbf56c16c1619b84a,
98                0x3efa013d3d01cf7f,
99            ];
100            let p01 = f_fmla(x2, f64::from_bits(COS_C[1]), f64::from_bits(COS_C[0]));
101            let p23 = f_fmla(x2, f64::from_bits(COS_C[3]), f64::from_bits(COS_C[2]));
102            let w0 = f_fmla(x4, p23, p01);
103            let w1 = x2 * w0;
104            let r_cos = DoubleDouble::from_exact_add(1., w1);
105
106            let err = f_fmla(
107                x2,
108                f64::from_bits(0x3cb0000000000000), // 2^-52
109                f64::from_bits(0x3be0000000000000), // 2^-65
110            );
111
112            let sin_ub = r_sin.hi + (r_sin.lo + err);
113            let sin_lb = r_sin.hi + (r_sin.lo - err);
114            let sin_x = if sin_ub == sin_lb {
115                sin_ub
116            } else {
117                sin_accurate_near_zero(x)
118            };
119
120            let cos_ub = r_cos.hi + (r_cos.lo + err);
121            let cos_lb = r_cos.hi + (r_cos.lo - err);
122            let cos_x = if cos_ub == cos_lb {
123                cos_ub
124            } else {
125                cos_accurate_near_zero(x)
126            };
127            return (sin_x, cos_x);
128        } else {
129            // // Small range reduction.
130            (y, k) = range_reduction_small(x);
131        }
132    } else {
133        // Inf or NaN
134        if x_e > 2 * E_BIAS {
135            // sin(+-Inf) = NaN
136            return (x + f64::NAN, x + f64::NAN);
137        }
138
139        // Large range reduction.
140        (k, y) = argument_reduction.reduce(x);
141    }
142
143    let r_sincos = sincos_eval(y);
144    let (sin_y, cos_y) = (r_sincos.v_sin, r_sincos.v_cos);
145
146    // Fast look up version, but needs 256-entry table.
147    // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
148    let sk = SIN_K_PI_OVER_128[(k & 255) as usize];
149    let ck = SIN_K_PI_OVER_128[((k.wrapping_add(64)) & 255) as usize];
150    let sin_k = DoubleDouble::from_bit_pair(sk);
151    let cos_k = DoubleDouble::from_bit_pair(ck);
152
153    let msin_k = -sin_k;
154
155    // After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128).
156    // So k is an integer and -pi / 256 <= y <= pi / 256.
157    // Then sin(x) = sin((k * pi/128 + y)
158    //             = sin(y) * cos(k*pi/128) + cos(y) * sin(k*pi/128)
159    let sin_k_cos_y = DoubleDouble::quick_mult(cos_y, sin_k);
160    let cos_k_sin_y = DoubleDouble::quick_mult(sin_y, cos_k);
161    //      cos(x) = cos((k * pi/128 + y)
162    //             = cos(y) * cos(k*pi/128) - sin(y) * sin(k*pi/128)
163    let cos_k_cos_y = DoubleDouble::quick_mult(cos_y, cos_k);
164    let msin_k_sin_y = DoubleDouble::quick_mult(sin_y, msin_k);
165
166    // cos_k_sin_y is always >> sin_k_cos_y
167    let mut sin_dd = DoubleDouble::from_exact_add(sin_k_cos_y.hi, cos_k_sin_y.hi);
168    // cos_k_cos_y is always >> msin_k_sin_y
169    let mut cos_dd = DoubleDouble::from_exact_add(cos_k_cos_y.hi, msin_k_sin_y.hi);
170    sin_dd.lo += sin_k_cos_y.lo + cos_k_sin_y.lo;
171    cos_dd.lo += msin_k_sin_y.lo + cos_k_cos_y.lo;
172
173    let sin_lp = sin_dd.lo + r_sincos.err;
174    let sin_lm = sin_dd.lo - r_sincos.err;
175    let cos_lp = cos_dd.lo + r_sincos.err;
176    let cos_lm = cos_dd.lo - r_sincos.err;
177
178    let sin_upper = sin_dd.hi + sin_lp;
179    let sin_lower = sin_dd.hi + sin_lm;
180    let cos_upper = cos_dd.hi + cos_lp;
181    let cos_lower = cos_dd.hi + cos_lm;
182
183    // Ziv's rounding test.
184    if sin_upper == sin_lower && cos_upper == cos_lower {
185        return (sin_upper, cos_upper);
186    }
187
188    sincos_hard(y, sin_k, cos_k, sin_upper, sin_lower, cos_upper, cos_lower)
189}
190
191#[cold]
192#[inline(never)]
193#[allow(clippy::too_many_arguments)]
194fn sincos_hard(
195    y: DoubleDouble,
196    sin_k: DoubleDouble,
197    cos_k: DoubleDouble,
198    sin_upper: f64,
199    sin_lower: f64,
200    cos_upper: f64,
201    cos_lower: f64,
202) -> (f64, f64) {
203    let r_sincos = sincos_eval_dd(y);
204
205    let msin_k = -sin_k;
206
207    let sin_x = if sin_upper == sin_lower {
208        sin_upper
209    } else {
210        // sin(x) = sin((k * pi/128 + u)
211        //        = sin(u) * cos(k*pi/128) + cos(u) * sin(k*pi/128)
212
213        DoubleDouble::mul_add(sin_k, r_sincos.v_cos, cos_k * r_sincos.v_sin).to_f64()
214    };
215
216    let cos_x = if cos_upper == cos_lower {
217        cos_upper
218    } else {
219        // cos(x) = cos((k * pi/128 + u)
220        //        = cos(u) * cos(k*pi/128) - sin(u) * sin(k*pi/128)
221        DoubleDouble::mul_add(cos_k, r_sincos.v_cos, msin_k * r_sincos.v_sin).to_f64()
222    };
223    (sin_x, cos_x)
224}
225
226#[cfg(test)]
227mod tests {
228    use super::*;
229
230    #[test]
231    fn f_sincos_test() {
232        let subnormal = f_sincos(0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000015708065354637772);
233        assert_eq!(subnormal.0, 1.5708065354637772e-307);
234        assert_eq!(subnormal.1, 1.0);
235        let zx_0 = f_sincos(0.0);
236        assert_eq!(zx_0.0, 0.0);
237        assert_eq!(zx_0.1, 1.0);
238        let zx_1 = f_sincos(1.0);
239        assert_eq!(zx_1.0, 0.8414709848078965);
240        assert_eq!(zx_1.1, 0.5403023058681398);
241        let zx_0_p5 = f_sincos(-0.5);
242        assert_eq!(zx_0_p5.0, -0.479425538604203);
243        assert_eq!(zx_0_p5.1, 0.8775825618903728);
244
245        let zx_1 = f_sincos(0.002341235432);
246        assert_eq!(zx_1.0, 0.0023412332931324344);
247        assert_eq!(zx_1.1, 0.9999972593095778);
248
249        let zx_1 = f_sincos(0.0198676543432);
250        assert_eq!(zx_1.0, 0.019866347330026367);
251        assert_eq!(zx_1.1, 0.9998026446473137);
252    }
253}
pxfm/sincos.rs

pxfm/
sincos.rs
