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]
170pub(crate) fn sincospif_eval(x: f64) -> SinCosf {
171    let argument_reduction = ArgumentReducerPi { x };
172
173    let (y, k) = argument_reduction.reduce();
174
175    let y_sqr = y * y;
176
177    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
178    // So k is an integer and -0.5 <= y <= 0.5.
179    // Then sin(x) = sin((k + y)*pi/32)
180    //             = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
181
182    let sin_k = f64::from_bits(SIN_K_PI_OVER32[((k as u64) & 63) as usize]);
183    // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).
184    // cos_k = cos(k * pi/32)
185    let cos_k = f64::from_bits(SIN_K_PI_OVER32[((k.wrapping_add(16) as u64) & 63) as usize]);
186
187    // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
188    // with:
189    // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
190    //
191    // See ./notes/sincosf.sollya
192    let sin_y = y * f_polyeval4(
193        y_sqr,
194        f64::from_bits(0x3fb921fb54442d18),
195        f64::from_bits(0xbf24abbce625abb1),
196        f64::from_bits(0x3e7466bc624f2776),
197        f64::from_bits(0xbdb32c3a619d4a7e),
198    );
199    // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
200    // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);
201    // Note that cosm1_y = cos(y*pi/32) - 1.
202    //
203    // See ./notes/sincosf.sollya
204    let cosm1_y = y_sqr
205        * f_polyeval3(
206            y_sqr,
207            f64::from_bits(0xbf73bd3cc9be430b),
208            f64::from_bits(0x3ed03c1f070c2e27),
209            f64::from_bits(0xbe155cc84bd94200),
210        );
211
212    SinCosf {
213        sin_k,
214        cos_k,
215        sin_y,
216        cosm1_y,
217    }
218}
219
220/**
221Poly for sin(pi*y) on [-0.25; 0.25].
222Generated by Sollya:
223```text
224d = [0,0.25];
225f_sin = sin(y*pi)/y;
226Q = fpminimax(f_sin, [|0, 2, 4, 6, 8, 10|], [|D...|], d, relative, floating);
227```
228See ./notes/sinpif_0p25.sollya
229**/
230#[inline]
231pub(crate) fn sinpif_eval(y: f64) -> f64 {
232    let y2 = y * y;
233    f_polyeval6(
234        y2,
235        f64::from_bits(0x400921fb54442cf8),
236        f64::from_bits(0xc014abbce6257778),
237        f64::from_bits(0x400466bc670fa464),
238        f64::from_bits(0xbfe32d2c627f47da),
239        f64::from_bits(0x3fb5071be0a2d3da),
240        f64::from_bits(0xbf7dd4e0ae5b9fcd),
241    ) * y
242}
243
244/**
245Poly for sin(pi*y) on [-0.25; 0.25].
246Generated by Sollya:
247```text
248d = [0,0.25];
249f_sin = sin(y*pi)/y;
250Q = fpminimax(f_sin, [|0, 2, 4, 6, 8, 10, 12|], [|D...|], d, relative, floating);
251```
252See ./notes/sinpif_0p25_2.sollya
253**/
254#[inline]
255pub(crate) fn sinpif_eval2(y: f64) -> f64 {
256    let y2 = y * y;
257    f_estrin_polyeval7(
258        y2,
259        f64::from_bits(0x400921fb54442d18),
260        f64::from_bits(0xc014abbce625be09),
261        f64::from_bits(0x400466bc67754fff),
262        f64::from_bits(0xbfe32d2ccdfe9424),
263        f64::from_bits(0x3fb50782d5f14825),
264        f64::from_bits(0xbf7e2fe76fdffd2b),
265        f64::from_bits(0x3f3e357ef99eb0bb),
266    ) * y
267}
268
269/**
270Poly for cos(pi*y) on [-0.25; 0.25].
271Generated by Sollya:
272```text
273d = [0, 0.25];
274f_cos = cos(y*pi);
275Q = fpminimax(f_cos, [|0, 2, 4, 6, 8, 10|], [|1, D...|], d, relative, floating);
276```
277See ./notes/cospif_0p25.sollya
278**/
279#[inline]
280pub(crate) fn cospif_eval(y: f64) -> f64 {
281    let y2 = y * y;
282    f_estrin_polyeval7(
283        y2,
284        f64::from_bits(0x3ff0000000000000),
285        f64::from_bits(0xc013bd3cc9be459b),
286        f64::from_bits(0x40103c1f081b0c98),
287        f64::from_bits(0xbff55d3c7dc25308),
288        f64::from_bits(0x3fce1f4fb6e12563),
289        f64::from_bits(0xbf9a6c9c101c1182),
290        f64::from_bits(0x3f5f3d9faffda250),
291    )
292}
pxfm/sin_cosf/sincosf_eval.rs

pxfm/sin_cosf/
sincosf_eval.rs
