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

pxfm/bessel/
alpha1.rs
