pxfm/bessel/
trigo_bessel.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 9/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::{f_fmla, min_normal_f64};
30use crate::round::RoundFinite;
31
32// Generated by SageMath:
33// print("[")
34// for k in range(128):
35//     k = RealField(150)(k) * RealField(150).pi() / RealField(150)(64)
36//     print(double_to_hex(k.sin()) + ",")
37// print("];")
38pub(crate) static SIN_K_PI_OVER_64: [u64; 128] = [
39    0x0000000000000000,
40    0x3fa91f65f10dd814,
41    0x3fb917a6bc29b42c,
42    0x3fc2c8106e8e613a,
43    0x3fc8f8b83c69a60b,
44    0x3fcf19f97b215f1b,
45    0x3fd294062ed59f06,
46    0x3fd58f9a75ab1fdd,
47    0x3fd87de2a6aea963,
48    0x3fdb5d1009e15cc0,
49    0x3fde2b5d3806f63b,
50    0x3fe073879922ffee,
51    0x3fe1c73b39ae68c8,
52    0x3fe30ff7fce17035,
53    0x3fe44cf325091dd6,
54    0x3fe57d69348ceca0,
55    0x3fe6a09e667f3bcd,
56    0x3fe7b5df226aafaf,
57    0x3fe8bc806b151741,
58    0x3fe9b3e047f38741,
59    0x3fea9b66290ea1a3,
60    0x3feb728345196e3e,
61    0x3fec38b2f180bdb1,
62    0x3feced7af43cc773,
63    0x3fed906bcf328d46,
64    0x3fee212104f686e5,
65    0x3fee9f4156c62dda,
66    0x3fef0a7efb9230d7,
67    0x3fef6297cff75cb0,
68    0x3fefa7557f08a517,
69    0x3fefd88da3d12526,
70    0x3feff621e3796d7e,
71    0x3ff0000000000000,
72    0x3feff621e3796d7e,
73    0x3fefd88da3d12526,
74    0x3fefa7557f08a517,
75    0x3fef6297cff75cb0,
76    0x3fef0a7efb9230d7,
77    0x3fee9f4156c62dda,
78    0x3fee212104f686e5,
79    0x3fed906bcf328d46,
80    0x3feced7af43cc773,
81    0x3fec38b2f180bdb1,
82    0x3feb728345196e3e,
83    0x3fea9b66290ea1a3,
84    0x3fe9b3e047f38741,
85    0x3fe8bc806b151741,
86    0x3fe7b5df226aafaf,
87    0x3fe6a09e667f3bcd,
88    0x3fe57d69348ceca0,
89    0x3fe44cf325091dd6,
90    0x3fe30ff7fce17035,
91    0x3fe1c73b39ae68c8,
92    0x3fe073879922ffee,
93    0x3fde2b5d3806f63b,
94    0x3fdb5d1009e15cc0,
95    0x3fd87de2a6aea963,
96    0x3fd58f9a75ab1fdd,
97    0x3fd294062ed59f06,
98    0x3fcf19f97b215f1b,
99    0x3fc8f8b83c69a60b,
100    0x3fc2c8106e8e613a,
101    0x3fb917a6bc29b42c,
102    0x3fa91f65f10dd814,
103    0xb69f77598338bfdf,
104    0xbfa91f65f10dd814,
105    0xbfb917a6bc29b42c,
106    0xbfc2c8106e8e613a,
107    0xbfc8f8b83c69a60b,
108    0xbfcf19f97b215f1b,
109    0xbfd294062ed59f06,
110    0xbfd58f9a75ab1fdd,
111    0xbfd87de2a6aea963,
112    0xbfdb5d1009e15cc0,
113    0xbfde2b5d3806f63b,
114    0xbfe073879922ffee,
115    0xbfe1c73b39ae68c8,
116    0xbfe30ff7fce17035,
117    0xbfe44cf325091dd6,
118    0xbfe57d69348ceca0,
119    0xbfe6a09e667f3bcd,
120    0xbfe7b5df226aafaf,
121    0xbfe8bc806b151741,
122    0xbfe9b3e047f38741,
123    0xbfea9b66290ea1a3,
124    0xbfeb728345196e3e,
125    0xbfec38b2f180bdb1,
126    0xbfeced7af43cc773,
127    0xbfed906bcf328d46,
128    0xbfee212104f686e5,
129    0xbfee9f4156c62dda,
130    0xbfef0a7efb9230d7,
131    0xbfef6297cff75cb0,
132    0xbfefa7557f08a517,
133    0xbfefd88da3d12526,
134    0xbfeff621e3796d7e,
135    0xbff0000000000000,
136    0xbfeff621e3796d7e,
137    0xbfefd88da3d12526,
138    0xbfefa7557f08a517,
139    0xbfef6297cff75cb0,
140    0xbfef0a7efb9230d7,
141    0xbfee9f4156c62dda,
142    0xbfee212104f686e5,
143    0xbfed906bcf328d46,
144    0xbfeced7af43cc773,
145    0xbfec38b2f180bdb1,
146    0xbfeb728345196e3e,
147    0xbfea9b66290ea1a3,
148    0xbfe9b3e047f38741,
149    0xbfe8bc806b151741,
150    0xbfe7b5df226aafaf,
151    0xbfe6a09e667f3bcd,
152    0xbfe57d69348ceca0,
153    0xbfe44cf325091dd6,
154    0xbfe30ff7fce17035,
155    0xbfe1c73b39ae68c8,
156    0xbfe073879922ffee,
157    0xbfde2b5d3806f63b,
158    0xbfdb5d1009e15cc0,
159    0xbfd87de2a6aea963,
160    0xbfd58f9a75ab1fdd,
161    0xbfd294062ed59f06,
162    0xbfcf19f97b215f1b,
163    0xbfc8f8b83c69a60b,
164    0xbfc2c8106e8e613a,
165    0xbfb917a6bc29b42c,
166    0xbfa91f65f10dd814,
167];
168
169#[inline]
170pub(crate) fn reduce_small_pi64(x: f64) -> (f64, i64) {
171    // Generated in SageMath:
172    // z = RealField(300)(64) / RealField(300).pi()
173    // n = 32
174    // x_hi = RealField(n)(z)  # convert to f64
175    // x_mid = RealField(n)(z - RealField(300)(x_hi))
176    // x_lo = RealField(n)(z - RealField(300)(x_hi) - RealField(300)(x_mid))
177    // print(double_to_hex(x_hi), ",")
178    // print(double_to_hex(x_mid), ",")
179    // print(double_to_hex(x_lo), ",")
180    const MPI_OVER_SIXTY_FOUR: [u64; 3] =
181        [0xbfa921fb54400000, 0xbd80b4611a600000, 0xbb53198a2e037073];
182    const SIXTY_EIGHT_OVER_PI: f64 = f64::from_bits(0x40345f306dc9c883);
183    let prod_hi = x * SIXTY_EIGHT_OVER_PI;
184    let kd = prod_hi.round_finite();
185
186    // Let y = x - k * (pi/64)
187    // Then |y| < pi / 64
188    // With extra rounding errors, we can bound |y| < 1.6 * 2^-7.
189    let y_hi = f_fmla(kd, f64::from_bits(MPI_OVER_SIXTY_FOUR[0]), x); // Exact
190    // |u.hi| < 1.6*2^-7
191    let u_hi = f_fmla(kd, f64::from_bits(MPI_OVER_SIXTY_FOUR[1]), y_hi);
192    (u_hi, unsafe {
193        kd.to_int_unchecked::<i64>() // indeterminate values is always filtered out before this call, as well only lowest bits are used
194    })
195}
196
197struct SinCosPi64 {
198    v_sin: f64,
199    v_cos: f64,
200}
201
202#[inline]
203fn sincos_eval_pi64(x: f64) -> SinCosPi64 {
204    let x2 = x * x;
205    let x4 = x2 * x2;
206
207    // Sin poly generated by Sollya:
208    // d = [0, pi/64];
209    // f_sin = sin(x)/x;
210    // Q = fpminimax(f_sin, [|0, 2, 4, 6|], [|D...|], d);
211    const S: [u64; 4] = [
212        0x3ff0000000000000,
213        0xbfc5555555555451,
214        0x3f8111111072c563,
215        0xbf2a01321c030841,
216    ];
217
218    let s0 = f_fmla(x2, f64::from_bits(S[1]), f64::from_bits(S[0]));
219    let s1 = f_fmla(x2, f64::from_bits(S[3]), f64::from_bits(S[2]));
220    let v_sin = f_fmla(x4, s1, s0) * x;
221
222    // Cos poly generated by Sollya:
223    // d = [0, pi/64];
224    // f_cos = cos(x);
225    // Q = fpminimax(f_cos, [|0, 2, 4, 6|], [|1, D...|], d);
226    const C: [u64; 4] = [
227        0x3ff0000000000000,
228        0xbfdffffffffffb6c,
229        0x3fa5555553f117c1,
230        0xbf56c0f056672a03,
231    ];
232    let c0 = f_fmla(x2, f64::from_bits(C[1]), f64::from_bits(C[0]));
233    let c1 = f_fmla(x2, f64::from_bits(C[3]), f64::from_bits(C[2]));
234    let v_cos = f_fmla(x4, c1, c0);
235
236    SinCosPi64 { v_sin, v_cos }
237}
238
239#[inline]
240pub(crate) fn sin_small(z: f64) -> f64 {
241    let x_e = (z.to_bits() >> 52) & 0x7ff;
242    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
243
244    if x_e < E_BIAS - 26 {
245        return f_fmla(z, f64::from_bits(0xbc90000000000000), z);
246    }
247
248    let (angle_dd, k) = reduce_small_pi64(z);
249    let sin_cos = sincos_eval_pi64(angle_dd);
250
251    // cos(k * pi/64) = sin(k * pi/64 + pi/2) = sin((k + 64) * pi/64).
252    let sk = SIN_K_PI_OVER_64[((k as u64) & 127) as usize];
253    let ck = SIN_K_PI_OVER_64[(((k as u64).wrapping_add(32)) & 127) as usize];
254
255    let sin_k = f64::from_bits(sk);
256    let cos_k = f64::from_bits(ck);
257    f_fmla(sin_cos.v_cos, sin_k, sin_cos.v_sin * cos_k)
258}
259
260#[inline]
261pub(crate) fn cos_small(z: f64) -> f64 {
262    let x_e = (z.to_bits() >> 52) & 0x7ff;
263    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
264
265    if x_e < E_BIAS - 27 {
266        // Signed zeros.
267        if z == 0.0 {
268            return 1.0;
269        }
270        // For |x| < 2^-26, |sin(x) - x| < ulp(x)/2.
271        return 1.0 - min_normal_f64();
272    }
273
274    let (angle_dd, k) = reduce_small_pi64(z);
275
276    let sin_cos = sincos_eval_pi64(angle_dd);
277
278    // cos(k * pi/64) = sin(k * pi/64 + pi/2) = sin((k + 64) * pi/64).
279    let sk = SIN_K_PI_OVER_64[((k as u64).wrapping_add(64) & 127) as usize];
280    let ck = SIN_K_PI_OVER_64[(((k as u64).wrapping_add(32)) & 127) as usize];
281
282    let sin_k = f64::from_bits(sk);
283    let cos_k = f64::from_bits(ck);
284    f_fmla(sin_cos.v_cos, cos_k, sin_cos.v_sin * sin_k)
285}
pxfm/bessel/trigo_bessel.rs

pxfm/bessel/
trigo_bessel.rs
