pxfm/
sin_helper.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 7/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::double_double::DoubleDouble;
30use crate::dyadic_float::DyadicFloat128;
31use crate::round::RoundFinite;
32use crate::sin::{SinCos, get_sin_k_rational, sincos_eval};
33use crate::sin_table::SIN_K_PI_OVER_128;
34use crate::sincos_dyadic::{range_reduction_small_f128_f128, sincos_eval_dyadic};
35
36#[inline]
37fn sin_eval_dd(z: DoubleDouble) -> DoubleDouble {
38    const SIN_C: [(u64, u64); 5] = [
39        (0x0000000000000000, 0x3ff0000000000000),
40        (0xbc65555555555555, 0xbfc5555555555555),
41        (0x3c01111111111111, 0x3f81111111111111),
42        (0xbb6a01a01a01a01a, 0xbf2a01a01a01a01a),
43        (0xbb6c154f8ddc6c00, 0x3ec71de3a556c734),
44    ];
45    let x2 = DoubleDouble::quick_mult(z, z);
46    let mut p = DoubleDouble::mul_add(
47        x2,
48        DoubleDouble::from_bit_pair(SIN_C[4]),
49        DoubleDouble::from_bit_pair(SIN_C[3]),
50    );
51
52    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(SIN_C[2]));
53    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(SIN_C[1]));
54    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(SIN_C[0]));
55    DoubleDouble::quick_mult(p, z)
56}
57
58pub(crate) fn sin_dd_small(z: DoubleDouble) -> DoubleDouble {
59    let x_e = (z.hi.to_bits() >> 52) & 0x7ff;
60    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
61
62    if x_e < E_BIAS - 8 {
63        return sin_eval_dd(z);
64    }
65
66    let (u_f128, k) = range_reduction_small_dd(z);
67
68    let sin_cos = sincos_eval_dd(u_f128);
69
70    // Fast look up version, but needs 256-entry table.
71    // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
72    let sk = SIN_K_PI_OVER_128[(k & 255) as usize];
73    let ck = SIN_K_PI_OVER_128[((k.wrapping_add(64)) & 255) as usize];
74
75    let sin_k = DoubleDouble::from_bit_pair(sk);
76    let cos_k = DoubleDouble::from_bit_pair(ck);
77
78    let sin_k_cos_y = DoubleDouble::quick_mult(sin_cos.v_cos, sin_k);
79    let cos_k_sin_y = DoubleDouble::quick_mult(sin_cos.v_sin, cos_k);
80
81    // sin_k_cos_y is always >> cos_k_sin_y
82    let mut rr = DoubleDouble::from_exact_add(sin_k_cos_y.hi, cos_k_sin_y.hi);
83    rr.lo += sin_k_cos_y.lo + cos_k_sin_y.lo;
84    rr
85}
86
87pub(crate) fn sin_dd_small_fast(z: DoubleDouble) -> DoubleDouble {
88    let x_e = (z.hi.to_bits() >> 52) & 0x7ff;
89    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
90
91    if x_e < E_BIAS - 8 {
92        return sin_eval_dd(z);
93    }
94
95    let (u_f128, k) = range_reduction_small_dd(z);
96
97    let sin_cos = sincos_eval(u_f128);
98
99    // Fast look up version, but needs 256-entry table.
100    // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
101    let sk = SIN_K_PI_OVER_128[(k & 255) as usize];
102    let ck = SIN_K_PI_OVER_128[((k.wrapping_add(64)) & 255) as usize];
103
104    let sin_k = DoubleDouble::from_bit_pair(sk);
105    let cos_k = DoubleDouble::from_bit_pair(ck);
106
107    let sin_k_cos_y = DoubleDouble::quick_mult(sin_cos.v_cos, sin_k);
108    let cos_k_sin_y = DoubleDouble::quick_mult(sin_cos.v_sin, cos_k);
109
110    // sin_k_cos_y is always >> cos_k_sin_y
111    let mut rr = DoubleDouble::from_exact_add(sin_k_cos_y.hi, cos_k_sin_y.hi);
112    rr.lo += sin_k_cos_y.lo + cos_k_sin_y.lo;
113    rr
114}
115
116#[inline]
117fn cos_eval_dd(z: DoubleDouble) -> DoubleDouble {
118    let x2 = DoubleDouble::quick_mult(z, z);
119    const COS_C: [(u64, u64); 5] = [
120        (0x0000000000000000, 0x3ff0000000000000),
121        (0x0000000000000000, 0xbfe0000000000000),
122        (0x3c45555555555555, 0x3fa5555555555555),
123        (0x3bef49f49f49f49f, 0xbf56c16c16c16c17),
124        (0x3b3a01a01a01a01a, 0x3efa01a01a01a01a),
125    ];
126
127    let mut p = DoubleDouble::mul_add(
128        x2,
129        DoubleDouble::from_bit_pair(COS_C[4]),
130        DoubleDouble::from_bit_pair(COS_C[3]),
131    );
132
133    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(COS_C[2]));
134    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(COS_C[1]));
135    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(COS_C[0]));
136
137    p
138}
139
140pub(crate) fn cos_dd_small(z: DoubleDouble) -> DoubleDouble {
141    let x_e = (z.hi.to_bits() >> 52) & 0x7ff;
142    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
143
144    if x_e < E_BIAS - 8 {
145        return cos_eval_dd(z);
146    }
147
148    let (u_f128, k) = range_reduction_small_dd(z);
149
150    let sin_cos = sincos_eval_dd(u_f128);
151
152    // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
153    let sk = SIN_K_PI_OVER_128[(k.wrapping_add(128) & 255) as usize];
154    let ck = SIN_K_PI_OVER_128[((k.wrapping_add(64)) & 255) as usize];
155    let msin_k = DoubleDouble::from_bit_pair(sk);
156    let cos_k = DoubleDouble::from_bit_pair(ck);
157
158    let cos_k_cos_y = DoubleDouble::quick_mult(sin_cos.v_cos, cos_k);
159    let cos_k_msin_y = DoubleDouble::quick_mult(sin_cos.v_sin, msin_k);
160
161    // cos_k_cos_y is always >> cos_k_msin_y
162    let mut rr = DoubleDouble::from_exact_add(cos_k_cos_y.hi, cos_k_msin_y.hi);
163    rr.lo += cos_k_cos_y.lo + cos_k_msin_y.lo;
164
165    rr
166}
167
168pub(crate) fn cos_dd_small_fast(z: DoubleDouble) -> DoubleDouble {
169    let x_e = (z.hi.to_bits() >> 52) & 0x7ff;
170    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
171
172    if x_e < E_BIAS - 8 {
173        return cos_eval_dd(z);
174    }
175
176    let (u_f128, k) = range_reduction_small_dd(z);
177
178    let sin_cos = sincos_eval(u_f128);
179
180    // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
181    let sk = SIN_K_PI_OVER_128[(k.wrapping_add(128) & 255) as usize];
182    let ck = SIN_K_PI_OVER_128[((k.wrapping_add(64)) & 255) as usize];
183    let msin_k = DoubleDouble::from_bit_pair(sk);
184    let cos_k = DoubleDouble::from_bit_pair(ck);
185
186    let cos_k_cos_y = DoubleDouble::quick_mult(sin_cos.v_cos, cos_k);
187    let cos_k_msin_y = DoubleDouble::quick_mult(sin_cos.v_sin, msin_k);
188
189    // cos_k_cos_y is always >> cos_k_msin_y
190    let mut rr = DoubleDouble::from_exact_add(cos_k_cos_y.hi, cos_k_msin_y.hi);
191    rr.lo += cos_k_cos_y.lo + cos_k_msin_y.lo;
192
193    rr
194}
195
196pub(crate) fn sin_f128_small(z: DyadicFloat128) -> DyadicFloat128 {
197    let (u_f128, k) = range_reduction_small_f128_f128(z);
198
199    let sin_cos = sincos_eval_dyadic(&u_f128);
200    // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
201    let sin_k_f128 = get_sin_k_rational(k as u64);
202    let cos_k_f128 = get_sin_k_rational((k as u64).wrapping_add(64));
203
204    // sin(x) = sin(k * pi/128 + u)
205    //        = sin(u) * cos(k*pi/128) + cos(u) * sin(k*pi/128)
206    (sin_k_f128 * sin_cos.v_cos) + (cos_k_f128 * sin_cos.v_sin)
207}
208
209pub(crate) fn cos_f128_small(z: DyadicFloat128) -> DyadicFloat128 {
210    let (u_f128, k) = range_reduction_small_f128_f128(z);
211
212    let sin_cos = sincos_eval_dyadic(&u_f128);
213    // -sin(k * pi/128) = sin((k + 128) * pi/128)
214    // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
215    let msin_k_f128 = get_sin_k_rational((k as u64).wrapping_add(128));
216    let cos_k_f128 = get_sin_k_rational((k as u64).wrapping_add(64));
217
218    // cos(x) = cos((k * pi/128 + u)
219    //        = cos(u) * cos(k*pi/128) - sin(u) * sin(k*pi/128)
220    (cos_k_f128 * sin_cos.v_cos) + (msin_k_f128 * sin_cos.v_sin)
221}
222
223#[inline]
224pub(crate) fn sincos_eval_dd(u: DoubleDouble) -> SinCos {
225    const SIN_C: [(u64, u64); 5] = [
226        (0x0000000000000000, 0x3ff0000000000000),
227        (0xbc65555555555555, 0xbfc5555555555555),
228        (0x3c01111111111111, 0x3f81111111111111),
229        (0xbb6a01a01a01a01a, 0xbf2a01a01a01a01a),
230        (0xbb6c154f8ddc6c00, 0x3ec71de3a556c734),
231    ];
232    let x2 = DoubleDouble::quick_mult(u, u);
233    let mut p = DoubleDouble::quick_mul_add(
234        x2,
235        DoubleDouble::from_bit_pair(SIN_C[4]),
236        DoubleDouble::from_bit_pair(SIN_C[3]),
237    );
238
239    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(SIN_C[2]));
240    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(SIN_C[1]));
241    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(SIN_C[0]));
242    let sin_u = DoubleDouble::quick_mult(p, u);
243
244    const COS_C: [(u64, u64); 5] = [
245        (0x0000000000000000, 0x3ff0000000000000),
246        (0x0000000000000000, 0xbfe0000000000000),
247        (0x3c45555555555555, 0x3fa5555555555555),
248        (0x3bef49f49f49f49f, 0xbf56c16c16c16c17),
249        (0x3b3a01a01a01a01a, 0x3efa01a01a01a01a),
250    ];
251
252    let mut p = DoubleDouble::quick_mul_add(
253        x2,
254        DoubleDouble::from_bit_pair(COS_C[4]),
255        DoubleDouble::from_bit_pair(COS_C[3]),
256    );
257
258    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(COS_C[2]));
259    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(COS_C[1]));
260    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(COS_C[0]));
261
262    let cos_u = p;
263    SinCos {
264        v_sin: sin_u,
265        v_cos: cos_u,
266        err: 0.,
267    }
268}
269
270#[inline]
271pub(crate) fn range_reduction_small_dd(x: DoubleDouble) -> (DoubleDouble, i64) {
272    const MPI_OVER_128: [u64; 5] = [
273        0xbf9921fb54400000,
274        0xbd70b4611a600000,
275        0xbb43198a2e000000,
276        0xb91b839a25200000,
277        0xb6b2704453400000,
278    ];
279    const ONE_TWENTY_EIGHT_OVER_PI_D: f64 = f64::from_bits(0x40445f306dc9c883);
280    let prod_hi = DoubleDouble::quick_mult_f64(x, ONE_TWENTY_EIGHT_OVER_PI_D);
281    let kd = prod_hi.to_f64().round_finite();
282
283    let p_hi = f64::from_bits(MPI_OVER_128[0]); // hi
284    let p_mid = f64::from_bits(MPI_OVER_128[1]); // mid
285    let p_lo = f64::from_bits(MPI_OVER_128[2]); // lo
286    let p_lo_lo = f64::from_bits(MPI_OVER_128[3]); // lo_lo
287
288    let mut q = DoubleDouble::f64_mul_f64_add(kd, p_hi, x);
289    q = DoubleDouble::f64_mul_f64_add(kd, p_mid, q);
290    q = DoubleDouble::f64_mul_f64_add(kd, p_lo, q);
291    q = DoubleDouble::f64_mul_f64_add(kd, p_lo_lo, q);
292
293    (q, unsafe { kd.to_int_unchecked::<i64>() })
294}
pxfm/sin_helper.rs

pxfm/
sin_helper.rs
