pxfm/bessel/
alpha0.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::double_double::DoubleDouble;
30use crate::dyadic_float::{DyadicFloat128, DyadicSign};
31use crate::polyeval::f_polyeval9;
32
33//
34/// See [bessel_0_asympt_alpha] for the info
35pub(crate) fn bessel_0_asympt_alpha_hard(reciprocal: DyadicFloat128) -> DyadicFloat128 {
36    static C: [DyadicFloat128; 18] = [
37        DyadicFloat128 {
38            sign: DyadicSign::Pos,
39            exponent: -130,
40            mantissa: 0x80000000_00000000_00000000_00000000_u128,
41        },
42        DyadicFloat128 {
43            sign: DyadicSign::Neg,
44            exponent: -131,
45            mantissa: 0x85555555_55555555_55555555_55555555_u128,
46        },
47        DyadicFloat128 {
48            sign: DyadicSign::Pos,
49            exponent: -130,
50            mantissa: 0xd6999999_99999999_99999999_9999999a_u128,
51        },
52        DyadicFloat128 {
53            sign: DyadicSign::Neg,
54            exponent: -127,
55            mantissa: 0xd1ac2492_49249249_24924924_92492492_u128,
56        },
57        DyadicFloat128 {
58            sign: DyadicSign::Pos,
59            exponent: -123,
60            mantissa: 0xbbcd0fc7_1c71c71c_71c71c71_c71c71c7_u128,
61        },
62        DyadicFloat128 {
63            sign: DyadicSign::Neg,
64            exponent: -118,
65            mantissa: 0x85e8fe45_8ba2e8ba_2e8ba2e8_ba2e8ba3_u128,
66        },
67        DyadicFloat128 {
68            sign: DyadicSign::Pos,
69            exponent: -113,
70            mantissa: 0x8b5a8f33_63c4ec4e_c4ec4ec4_ec4ec4ec_u128,
71        },
72        DyadicFloat128 {
73            sign: DyadicSign::Neg,
74            exponent: -108,
75            mantissa: 0xc7661d79_9d59b555_55555555_55555555_u128,
76        },
77        DyadicFloat128 {
78            sign: DyadicSign::Pos,
79            exponent: -102,
80            mantissa: 0xbbced715_c2897a28_78787878_78787878_u128,
81        },
82        DyadicFloat128 {
83            sign: DyadicSign::Neg,
84            exponent: -96,
85            mantissa: 0xe14b19b4_aae3f7fe_be1af286_bca1af28_u128,
86        },
87        DyadicFloat128 {
88            sign: DyadicSign::Pos,
89            exponent: -89,
90            mantissa: 0xa7af7341_db2192db_975e0c30_c30c30c3_u128,
91        },
92        DyadicFloat128 {
93            sign: DyadicSign::Neg,
94            exponent: -82,
95            mantissa: 0x97a8f676_b349f6fc_5cefd338_590b2164_u128,
96        },
97        DyadicFloat128 {
98            sign: DyadicSign::Pos,
99            exponent: -75,
100            mantissa: 0xa3d299fb_6f304d73_86e15f12_0fd70a3d_u128,
101        },
102        DyadicFloat128 {
103            sign: DyadicSign::Neg,
104            exponent: -68,
105            mantissa: 0xd050b737_cbc044ef_e8807e3c_87f43da1_u128,
106        },
107        DyadicFloat128 {
108            sign: DyadicSign::Pos,
109            exponent: -60,
110            mantissa: 0x9a02379b_daa7e492_854f42de_6d3dffe6_u128,
111        },
112        DyadicFloat128 {
113            sign: DyadicSign::Neg,
114            exponent: -52,
115            mantissa: 0x83011a39_380e467d_de6b70ec_b92ce0cc_u128,
116        },
117        DyadicFloat128 {
118            sign: DyadicSign::Pos,
119            exponent: -45,
120            mantissa: 0xfe16521f_c79e5d9a_a5bed653_e3844e9a_u128,
121        },
122        DyadicFloat128 {
123            sign: DyadicSign::Neg,
124            exponent: -36,
125            mantissa: 0x8b54b13d_3fb3e1c4_15dbb880_0bb32218_u128,
126        },
127    ];
128
129    let x2 = reciprocal * reciprocal;
130
131    let mut p = C[17];
132    for i in (0..17).rev() {
133        p = x2 * p + C[i];
134    }
135
136    p * reciprocal
137}
138
139/**
140Note expansion generation below: this is negative series expressed in Sage as positive,
141so before any real evaluation `x=1/x` should be applied.
142
143Generated by SageMath:
144```python
145def binomial_like(n, m):
146    prod = QQ(1)
147    z = QQ(4)*(n**2)
148    for k in range(1,m + 1):
149        prod *= (z - (2*k - 1)**2)
150    return prod / (QQ(2)**(2*m) * (ZZ(m).factorial()))
151
152R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
153x = R.gen()
154
155def Pn_asymptotic(n, y, terms=10):
156    # now y = 1/x
157    return sum( (-1)**m * binomial_like(n, 2*m) / (QQ(2)**(2*m)) * y**(QQ(2)*m) for m in range(terms) )
158
159def Qn_asymptotic(n, y, terms=10):
160    return sum( (-1)**m * binomial_like(n, 2*m + 1) / (QQ(2)**(2*m + 1)) * y**(QQ(2)*m + 1) for m in range(terms) )
161
162P = Pn_asymptotic(0, x, 50)
163Q = Qn_asymptotic(0, x, 50)
164
165R_series = (-Q/P)
166
167# alpha is atan(R_series) so we're doing Taylor series atan expansion on R_series
168
169arctan_series_Z = sum([QQ(-1)**k * x**(QQ(2)*k+1) / RealField(700)(RealField(700)(2)*k+1) for k in range(25)])
170alpha_series = arctan_series_Z(R_series)
171
172# see the series
173print(alpha_series)
174```
175**/
176#[inline]
177pub(crate) fn bessel_0_asympt_alpha(recip: DoubleDouble) -> DoubleDouble {
178    const C: [(u64, u64); 12] = [
179        (0x0000000000000000, 0x3fc0000000000000),
180        (0x3c55555555555555, 0xbfb0aaaaaaaaaaab),
181        (0x3c5999999999999a, 0x3fcad33333333333),
182        (0xbc92492492492492, 0xbffa358492492492),
183        (0xbcbc71c71c71c71c, 0x403779a1f8e38e39),
184        (0xbd0745d1745d1746, 0xc080bd1fc8b1745d),
185        (0xbd7d89d89d89d89e, 0x40d16b51e66c789e),
186        (0x3dc5555555555555, 0xc128ecc3af33ab37),
187        (0x3e2143c3c3c3c3c4, 0x418779dae2b8512f),
188        (0x3df41e50d79435e5, 0xc1ec296336955c7f),
189        (0x3ef6dcbaf0618618, 0x4254f5ee683b6432),
190        (0x3f503a3102cc7a6f, 0xc2c2f51eced6693f),
191    ];
192
193    // Doing (1/x)*(1/x) instead (1/(x*x)) to avoid spurious overflow/underflow
194    let x2 = DoubleDouble::quick_mult(recip, recip);
195
196    let mut p = DoubleDouble::mul_add(
197        x2,
198        DoubleDouble::from_bit_pair(C[11]),
199        DoubleDouble::from_bit_pair(C[10]),
200    );
201
202    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[9]));
203    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[8]));
204    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[7]));
205    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[6]));
206    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[5]));
207    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[4]));
208    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[3]));
209    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[2]));
210    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[1]));
211    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[0].1));
212
213    let z = DoubleDouble::quick_mult(p, recip);
214
215    DoubleDouble::from_exact_add(z.hi, z.lo)
216}
217
218/**
219Note expansion generation below: this is negative series expressed in Sage as positive,
220so before any real evaluation `x=1/x` should be applied.
221
222Generated by SageMath:
223```python
224def binomial_like(n, m):
225    prod = QQ(1)
226    z = QQ(4)*(n**2)
227    for k in range(1,m + 1):
228        prod *= (z - (2*k - 1)**2)
229    return prod / (QQ(2)**(2*m) * (ZZ(m).factorial()))
230
231R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
232x = R.gen()
233
234def Pn_asymptotic(n, y, terms=10):
235    # now y = 1/x
236    return sum( (-1)**m * binomial_like(n, 2*m) / (QQ(2)**(2*m)) * y**(QQ(2)*m) for m in range(terms) )
237
238def Qn_asymptotic(n, y, terms=10):
239    return sum( (-1)**m * binomial_like(n, 2*m + 1) / (QQ(2)**(2*m + 1)) * y**(QQ(2)*m + 1) for m in range(terms) )
240
241P = Pn_asymptotic(0, x, 50)
242Q = Qn_asymptotic(0, x, 50)
243
244R_series = (-Q/P)
245
246# alpha is atan(R_series) so we're doing Taylor series atan expansion on R_series
247
248arctan_series_Z = sum([QQ(-1)**k * x**(QQ(2)*k+1) / RealField(700)(RealField(700)(2)*k+1) for k in range(25)])
249alpha_series = arctan_series_Z(R_series)
250
251# see the series
252print(alpha_series)
253```
254**/
255#[inline]
256pub(crate) fn bessel_0_asympt_alpha_fast(recip: DoubleDouble) -> DoubleDouble {
257    const C: [u64; 12] = [
258        0x3fc0000000000000,
259        0xbfb0aaaaaaaaaaab,
260        0x3fcad33333333333,
261        0xbffa358492492492,
262        0x403779a1f8e38e39,
263        0xc080bd1fc8b1745d,
264        0x40d16b51e66c789e,
265        0xc128ecc3af33ab37,
266        0x418779dae2b8512f,
267        0xc1ec296336955c7f,
268        0x4254f5ee683b6432,
269        0xc2c2f51eced6693f,
270    ];
271
272    // Doing (1/x)*(1/x) instead (1/(x*x)) to avoid spurious overflow/underflow
273    let x2 = DoubleDouble::quick_mult(recip, recip);
274
275    let p = f_polyeval9(
276        x2.hi,
277        f64::from_bits(C[3]),
278        f64::from_bits(C[4]),
279        f64::from_bits(C[5]),
280        f64::from_bits(C[6]),
281        f64::from_bits(C[7]),
282        f64::from_bits(C[8]),
283        f64::from_bits(C[9]),
284        f64::from_bits(C[10]),
285        f64::from_bits(C[11]),
286    );
287
288    let mut z = DoubleDouble::mul_f64_add_f64(x2, p, f64::from_bits(C[2]));
289    z = DoubleDouble::mul_add_f64(x2, z, f64::from_bits(C[1]));
290    z = DoubleDouble::mul_add_f64(x2, z, f64::from_bits(C[0]));
291    DoubleDouble::quick_mult(z, recip)
292}
pxfm/bessel/alpha0.rs

pxfm/bessel/
alpha0.rs
