pxfm/bessel/
beta1.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
59def sqrt_series(s):
60    val = S.valuation()
61    lc = S[val]  # Leading coefficient
62    b = lc.sqrt() * x**(val // 2)
63
64    for _ in range(5):
65        b = (b + S / b) / 2
66        b = b
67    return b
68
69S = (P**2 + Q**2).truncate(50)
70
71b_series = sqrt_series(S).truncate(30)
72# see the beta series
73print(b_series)
74```
75
76See notes/bessel_asympt.ipynb for generation
77**/
78#[inline]
79pub(crate) fn bessel_1_asympt_beta_fast(recip: DoubleDouble) -> DoubleDouble {
80    const C: [u64; 10] = [
81        0x3ff0000000000000,
82        0x3fc8000000000000,
83        0xbfc8c00000000000,
84        0x3fe9c50000000000,
85        0xc01ef5b680000000,
86        0x40609860dd400000,
87        0xc0abae9b7a06e000,
88        0x41008711d41c1428,
89        0xc15ab70164c8be6e,
90        0x41bc1055e24f297f,
91    ];
92
93    // Doing (1/x)*(1/x) instead (1/(x*x)) to avoid spurious overflow/underflow
94    let x2 = DoubleDouble::quick_mult(recip, recip);
95
96    let p = f_polyeval9(
97        x2.hi,
98        f64::from_bits(C[1]),
99        f64::from_bits(C[2]),
100        f64::from_bits(C[3]),
101        f64::from_bits(C[4]),
102        f64::from_bits(C[5]),
103        f64::from_bits(C[6]),
104        f64::from_bits(C[7]),
105        f64::from_bits(C[8]),
106        f64::from_bits(C[9]),
107    );
108
109    DoubleDouble::mul_f64_add_f64(x2, p, f64::from_bits(C[0]))
110}
111
112/**
113Note expansion generation below: this is negative series expressed in Sage as positive,
114so before any real evaluation `x=1/x` should be applied
115
116Generated by SageMath:
117```python
118def binomial_like(n, m):
119    prod = QQ(1)
120    z = QQ(4)*(n**2)
121    for k in range(1,m + 1):
122        prod *= (z - (2*k - 1)**2)
123    return prod / (QQ(2)**(2*m) * (ZZ(m).factorial()))
124
125R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
126x = R.gen()
127
128def Pn_asymptotic(n, y, terms=10):
129    # now y = 1/x
130    return sum( (-1)**m * binomial_like(n, 2*m) / (QQ(2)**(2*m)) * y**(QQ(2)*m) for m in range(terms) )
131
132def Qn_asymptotic(n, y, terms=10):
133    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) )
134
135P = Pn_asymptotic(1, x, 50)
136Q = Qn_asymptotic(1, x, 50)
137
138def sqrt_series(s):
139    val = S.valuation()
140    lc = S[val]  # Leading coefficient
141    b = lc.sqrt() * x**(val // 2)
142
143    for _ in range(5):
144        b = (b + S / b) / 2
145        b = b
146    return b
147
148S = (P**2 + Q**2).truncate(50)
149
150b_series = sqrt_series(S).truncate(30)
151# see the beta series
152print(b_series)
153```
154
155See notes/bessel_asympt.ipynb for generation
156**/
157#[inline]
158pub(crate) fn bessel_1_asympt_beta(recip: DoubleDouble) -> DoubleDouble {
159    const C: [(u64, u64); 10] = [
160        (0x0000000000000000, 0x3ff0000000000000), // 1
161        (0x0000000000000000, 0x3fc8000000000000), // 2
162        (0x0000000000000000, 0xbfc8c00000000000), // 3
163        (0x0000000000000000, 0x3fe9c50000000000), // 4
164        (0x0000000000000000, 0xc01ef5b680000000), // 5
165        (0x0000000000000000, 0x40609860dd400000), // 6
166        (0x0000000000000000, 0xc0abae9b7a06e000), // 7
167        (0x0000000000000000, 0x41008711d41c1428), // 8
168        (0xbdf7a00000000000, 0xc15ab70164c8be6e),
169        (0xbe40e1f000000000, 0x41bc1055e24f297f),
170    ];
171
172    // Doing (1/x)*(1/x) instead (1/(x*x)) to avoid spurious overflow/underflow
173    let x2 = DoubleDouble::quick_mult(recip, recip);
174
175    let mut p = DoubleDouble::mul_add(
176        x2,
177        DoubleDouble::from_bit_pair(C[9]),
178        DoubleDouble::from_bit_pair(C[8]),
179    );
180
181    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[7].1)); // 8
182    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[6].1)); // 7
183    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[5].1)); // 6
184    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[4].1)); // 5
185    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[3].1)); // 4
186    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[2].1)); // 3
187    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[1].1)); // 2
188    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[0].1)); // 1
189    p
190}
191
192/// see [bessel_1_asympt_beta] for more info
193pub(crate) fn bessel_1_asympt_beta_hard(recip: DyadicFloat128) -> DyadicFloat128 {
194    static C: [DyadicFloat128; 12] = [
195        DyadicFloat128 {
196            sign: DyadicSign::Pos,
197            exponent: -127,
198            mantissa: 0x80000000_00000000_00000000_00000000_u128,
199        },
200        DyadicFloat128 {
201            sign: DyadicSign::Pos,
202            exponent: -130,
203            mantissa: 0xc0000000_00000000_00000000_00000000_u128,
204        },
205        DyadicFloat128 {
206            sign: DyadicSign::Neg,
207            exponent: -130,
208            mantissa: 0xc6000000_00000000_00000000_00000000_u128,
209        },
210        DyadicFloat128 {
211            sign: DyadicSign::Pos,
212            exponent: -128,
213            mantissa: 0xce280000_00000000_00000000_00000000_u128,
214        },
215        DyadicFloat128 {
216            sign: DyadicSign::Neg,
217            exponent: -125,
218            mantissa: 0xf7adb400_00000000_00000000_00000000_u128,
219        },
220        DyadicFloat128 {
221            sign: DyadicSign::Pos,
222            exponent: -120,
223            mantissa: 0x84c306ea_00000000_00000000_00000000_u128,
224        },
225        DyadicFloat128 {
226            sign: DyadicSign::Neg,
227            exponent: -116,
228            mantissa: 0xdd74dbd0_37000000_00000000_00000000_u128,
229        },
230        DyadicFloat128 {
231            sign: DyadicSign::Pos,
232            exponent: -110,
233            mantissa: 0x84388ea0_e0a14000_00000000_00000000_u128,
234        },
235        DyadicFloat128 {
236            sign: DyadicSign::Neg,
237            exponent: -105,
238            mantissa: 0xd5b80b26_45f372f4_00000000_00000000_u128,
239        },
240        DyadicFloat128 {
241            sign: DyadicSign::Pos,
242            exponent: -99,
243            mantissa: 0xe082af12_794bf6f1_e1000000_00000000_u128,
244        },
245        DyadicFloat128 {
246            sign: DyadicSign::Neg,
247            exponent: -92,
248            mantissa: 0x94a06149_f30146bc_fe8ed000_00000000_u128,
249        },
250        DyadicFloat128 {
251            sign: DyadicSign::Pos,
252            exponent: -86,
253            mantissa: 0xf212edfc_42a62526_4fac2b0c_00000000_u128,
254        },
255    ];
256
257    let x2 = recip * recip;
258
259    let mut p = C[11];
260    for i in (0..11).rev() {
261        p = x2 * p + C[i];
262    }
263    p
264}
pxfm/bessel/beta1.rs

pxfm/bessel/
beta1.rs
