/rust/registry/src/index.crates.io-1949cf8c6b5b557f/pxfm-0.1.27/src/bessel/beta0.rs

Source
/*
 * // Copyright (c) Radzivon Bartoshyk 8/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
use crate::double_double::DoubleDouble;
use crate::dyadic_float::{DyadicFloat128, DyadicSign};
use crate::polyeval::f_polyeval9;

/**
Beta series

Generated by SageMath:
```python
#generate b series
def binomial_like(n, m):
    prod = QQ(1)
    z = QQ(4)*(n**2)
    for k in range(1,m + 1):
        prod *= (z - (2*k - 1)**2)
    return prod / (QQ(2)**(2*m) * (ZZ(m).factorial()))

R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
x = R.gen()

def Pn_asymptotic(n, y, terms=10):
    # now y = 1/x
    return sum( (-1)**m * binomial_like(n, 2*m) / (QQ(2)**(2*m)) * y**(QQ(2)*m) for m in range(terms) )

def Qn_asymptotic(n, y, terms=10):
    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) )

P = Pn_asymptotic(0, x, 50)
Q = Qn_asymptotic(0, x, 50)

def sqrt_series(s):
    val = S.valuation()
    lc = S[val]  # Leading coefficient
    b = lc.sqrt() * x**(val // 2)

    for _ in range(5):
        b = (b + S / b) / 2
        b = b
    return b

S = (P**2 + Q**2).truncate(50)

b_series = sqrt_series(S).truncate(30)
#see the series
print(b_series)
```
**/
#[inline]
pub(crate) fn bessel_0_asympt_beta(recip: DoubleDouble) -> DoubleDouble {
    const C: [(u64, u64); 10] = [
        (0x0000000000000000, 0x3ff0000000000000),
        (0x0000000000000000, 0xbfb0000000000000),
        (0x0000000000000000, 0x3fba800000000000),
        (0x0000000000000000, 0xbfe15f0000000000),
        (0x0000000000000000, 0x4017651180000000),
        (0x0000000000000000, 0xc05ab8c13b800000),
        (0x0000000000000000, 0x40a730492f262000),
        (0x0000000000000000, 0xc0fc73a7acd696f0),
        (0xbdf3a00000000000, 0x41577458dd9fce68),
        (0xbe4ba6b000000000, 0xc1b903ab9b27e18f),
    ];

    // Doing (1/x)*(1/x) instead (1/(x*x)) to avoid spurious overflow/underflow
    let x2 = DoubleDouble::quick_mult(recip, recip);

    let mut p = DoubleDouble::mul_add(
        x2,
        DoubleDouble::from_bit_pair(C[9]),
        DoubleDouble::from_bit_pair(C[8]),
    );

    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[7].1));
    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[6].1));
    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[5].1));
    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[4].1));
    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[3].1));
    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[2].1));
    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[1].1));
    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[0].1));
    p
}

/**
Beta series

Generated by SageMath:
```python
#generate b series
def binomial_like(n, m):
    prod = QQ(1)
    z = QQ(4)*(n**2)
    for k in range(1,m + 1):
        prod *= (z - (2*k - 1)**2)
    return prod / (QQ(2)**(2*m) * (ZZ(m).factorial()))

R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
x = R.gen()

def Pn_asymptotic(n, y, terms=10):
    # now y = 1/x
    return sum( (-1)**m * binomial_like(n, 2*m) / (QQ(2)**(2*m)) * y**(QQ(2)*m) for m in range(terms) )

def Qn_asymptotic(n, y, terms=10):
    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) )

P = Pn_asymptotic(0, x, 50)
Q = Qn_asymptotic(0, x, 50)

def sqrt_series(s):
    val = S.valuation()
    lc = S[val]  # Leading coefficient
    b = lc.sqrt() * x**(val // 2)

    for _ in range(5):
        b = (b + S / b) / 2
        b = b
    return b

S = (P**2 + Q**2).truncate(50)

b_series = sqrt_series(S).truncate(30)
#see the series
print(b_series)
```
**/
#[inline]
pub(crate) fn bessel_0_asympt_beta_fast(recip: DoubleDouble) -> DoubleDouble {
    const C: [u64; 10] = [
        0x3ff0000000000000,
        0xbfb0000000000000,
        0x3fba800000000000,
        0xbfe15f0000000000,
        0x4017651180000000,
        0xc05ab8c13b800000,
        0x40a730492f262000,
        0xc0fc73a7acd696f0,
        0x41577458dd9fce68,
        0xc1b903ab9b27e18f,
    ];

    // Doing (1/x)*(1/x) instead (1/(x*x)) to avoid spurious overflow/underflow
    let x2 = DoubleDouble::quick_mult(recip, recip);

    let p = f_polyeval9(
        x2.hi,
        f64::from_bits(C[1]),
        f64::from_bits(C[2]),
        f64::from_bits(C[3]),
        f64::from_bits(C[4]),
        f64::from_bits(C[5]),
        f64::from_bits(C[6]),
        f64::from_bits(C[7]),
        f64::from_bits(C[8]),
        f64::from_bits(C[9]),
    );

    DoubleDouble::mul_f64_add_f64(x2, p, f64::from_bits(C[0]))
}

/// see [bessel_0_asympt_beta] for more info
pub(crate) fn bessel_0_asympt_beta_hard(recip: DyadicFloat128) -> DyadicFloat128 {
    static C: [DyadicFloat128; 12] = [
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -127,
            mantissa: 0x80000000_00000000_00000000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -131,
            mantissa: 0x80000000_00000000_00000000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -131,
            mantissa: 0xd4000000_00000000_00000000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -128,
            mantissa: 0x8af80000_00000000_00000000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -125,
            mantissa: 0xbb288c00_00000000_00000000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -121,
            mantissa: 0xd5c609dc_00000000_00000000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -116,
            mantissa: 0xb9824979_31000000_00000000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -111,
            mantissa: 0xe39d3d66_b4b78000_00000000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -105,
            mantissa: 0xbba2c6ec_fe733d8c_00000000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -99,
            mantissa: 0xc81d5cd9_3f0c79ba_6b000000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -92,
            mantissa: 0x86118ddf_c1ffc100_0ee1b000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -86,
            mantissa: 0xdc7ccfa9_930b874d_52df3464_00000000_u128,
        },
    ];

    let x2 = recip * recip;

    let mut p = C[11];
    for i in (0..11).rev() {
        p = x2 * p + C[i];
    }
    p
}

Coverage Report

Created: 2026-01-22 07:28

Line	Count	Source
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		*/
29		use crate::double_double::DoubleDouble;
30		use crate::dyadic_float::{DyadicFloat128, DyadicSign};
31		use crate::polyeval::f_polyeval9;
32
33		/**
34		Beta series
35
36		Generated by SageMath:
37		```python
38		#generate b series
39		def 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 - (2k - 1)**2)
44		return prod / (QQ(2)*(2m) * (ZZ(m).factorial()))
45
46		R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
47		x = R.gen()
48
49		def Pn_asymptotic(n, y, terms=10):
50		# now y = 1/x
51		return sum( (-1)*m binomial_like(n, 2m) / (QQ(2)(2m)) * y*(QQ(2)m) for m in range(terms) )
52
53		def Qn_asymptotic(n, y, terms=10):
54		return sum( (-1)*m binomial_like(n, 2m + 1) / (QQ(2)(2m + 1)) * y*(QQ(2)m + 1) for m in range(terms) )
55
56		P = Pn_asymptotic(0, x, 50)
57		Q = Qn_asymptotic(0, x, 50)
58
59		def 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
69		S = (P2 + Q2).truncate(50)
70
71		b_series = sqrt_series(S).truncate(30)
72		#see the series
73		print(b_series)
74		```
75		**/
76		#[inline]
77	0	pub(crate) fn bessel_0_asympt_beta(recip: DoubleDouble) -> DoubleDouble {
78		const C: [(u64, u64); 10] = [
79		(0x0000000000000000, 0x3ff0000000000000),
80		(0x0000000000000000, 0xbfb0000000000000),
81		(0x0000000000000000, 0x3fba800000000000),
82		(0x0000000000000000, 0xbfe15f0000000000),
83		(0x0000000000000000, 0x4017651180000000),
84		(0x0000000000000000, 0xc05ab8c13b800000),
85		(0x0000000000000000, 0x40a730492f262000),
86		(0x0000000000000000, 0xc0fc73a7acd696f0),
87		(0xbdf3a00000000000, 0x41577458dd9fce68),
88		(0xbe4ba6b000000000, 0xc1b903ab9b27e18f),
89		];
90
91		// Doing (1/x)(1/x) instead (1/(xx)) to avoid spurious overflow/underflow
92	0	let x2 = DoubleDouble::quick_mult(recip, recip);
93
94	0	let mut p = DoubleDouble::mul_add(
95	0	x2,
96	0	DoubleDouble::from_bit_pair(C[9]),
97	0	DoubleDouble::from_bit_pair(C[8]),
98		);
99
100	0	p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[7].1));
101	0	p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[6].1));
102	0	p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[5].1));
103	0	p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[4].1));
104	0	p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[3].1));
105	0	p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[2].1));
106	0	p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[1].1));
107	0	p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(C[0].1));
108	0	p
109	0	}
110
111		/**
112		Beta series
113
114		Generated by SageMath:
115		```python
116		#generate b series
117		def 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 - (2k - 1)**2)
122		return prod / (QQ(2)*(2m) * (ZZ(m).factorial()))
123
124		R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
125		x = R.gen()
126
127		def Pn_asymptotic(n, y, terms=10):
128		# now y = 1/x
129		return sum( (-1)*m binomial_like(n, 2m) / (QQ(2)(2m)) * y*(QQ(2)m) for m in range(terms) )
130
131		def Qn_asymptotic(n, y, terms=10):
132		return sum( (-1)*m binomial_like(n, 2m + 1) / (QQ(2)(2m + 1)) * y*(QQ(2)m + 1) for m in range(terms) )
133
134		P = Pn_asymptotic(0, x, 50)
135		Q = Qn_asymptotic(0, x, 50)
136
137		def sqrt_series(s):
138		val = S.valuation()
139		lc = S[val] # Leading coefficient
140		b = lc.sqrt() * x**(val // 2)
141
142		for _ in range(5):
143		b = (b + S / b) / 2
144		b = b
145		return b
146
147		S = (P2 + Q2).truncate(50)
148
149		b_series = sqrt_series(S).truncate(30)
150		#see the series
151		print(b_series)
152		```
153		**/
154		#[inline]
155	0	pub(crate) fn bessel_0_asympt_beta_fast(recip: DoubleDouble) -> DoubleDouble {
156		const C: [u64; 10] = [
157		0x3ff0000000000000,
158		0xbfb0000000000000,
159		0x3fba800000000000,
160		0xbfe15f0000000000,
161		0x4017651180000000,
162		0xc05ab8c13b800000,
163		0x40a730492f262000,
164		0xc0fc73a7acd696f0,
165		0x41577458dd9fce68,
166		0xc1b903ab9b27e18f,
167		];
168
169		// Doing (1/x)(1/x) instead (1/(xx)) to avoid spurious overflow/underflow
170	0	let x2 = DoubleDouble::quick_mult(recip, recip);
171
172	0	let p = f_polyeval9(
173	0	x2.hi,
174	0	f64::from_bits(C[1]),
175	0	f64::from_bits(C[2]),
176	0	f64::from_bits(C[3]),
177	0	f64::from_bits(C[4]),
178	0	f64::from_bits(C[5]),
179	0	f64::from_bits(C[6]),
180	0	f64::from_bits(C[7]),
181	0	f64::from_bits(C[8]),
182	0	f64::from_bits(C[9]),
183		);
184
185	0	DoubleDouble::mul_f64_add_f64(x2, p, f64::from_bits(C[0]))
186	0	}
187
188		/// see [bessel_0_asympt_beta] for more info
189	0	pub(crate) fn bessel_0_asympt_beta_hard(recip: DyadicFloat128) -> DyadicFloat128 {
190		static C: [DyadicFloat128; 12] = [
191		DyadicFloat128 {
192		sign: DyadicSign::Pos,
193		exponent: -127,
194		mantissa: 0x80000000_00000000_00000000_00000000_u128,
195		},
196		DyadicFloat128 {
197		sign: DyadicSign::Neg,
198		exponent: -131,
199		mantissa: 0x80000000_00000000_00000000_00000000_u128,
200		},
201		DyadicFloat128 {
202		sign: DyadicSign::Pos,
203		exponent: -131,
204		mantissa: 0xd4000000_00000000_00000000_00000000_u128,
205		},
206		DyadicFloat128 {
207		sign: DyadicSign::Neg,
208		exponent: -128,
209		mantissa: 0x8af80000_00000000_00000000_00000000_u128,
210		},
211		DyadicFloat128 {
212		sign: DyadicSign::Pos,
213		exponent: -125,
214		mantissa: 0xbb288c00_00000000_00000000_00000000_u128,
215		},
216		DyadicFloat128 {
217		sign: DyadicSign::Neg,
218		exponent: -121,
219		mantissa: 0xd5c609dc_00000000_00000000_00000000_u128,
220		},
221		DyadicFloat128 {
222		sign: DyadicSign::Pos,
223		exponent: -116,
224		mantissa: 0xb9824979_31000000_00000000_00000000_u128,
225		},
226		DyadicFloat128 {
227		sign: DyadicSign::Neg,
228		exponent: -111,
229		mantissa: 0xe39d3d66_b4b78000_00000000_00000000_u128,
230		},
231		DyadicFloat128 {
232		sign: DyadicSign::Pos,
233		exponent: -105,
234		mantissa: 0xbba2c6ec_fe733d8c_00000000_00000000_u128,
235		},
236		DyadicFloat128 {
237		sign: DyadicSign::Neg,
238		exponent: -99,
239		mantissa: 0xc81d5cd9_3f0c79ba_6b000000_00000000_u128,
240		},
241		DyadicFloat128 {
242		sign: DyadicSign::Pos,
243		exponent: -92,
244		mantissa: 0x86118ddf_c1ffc100_0ee1b000_00000000_u128,
245		},
246		DyadicFloat128 {
247		sign: DyadicSign::Neg,
248		exponent: -86,
249		mantissa: 0xdc7ccfa9_930b874d_52df3464_00000000_u128,
250		},
251		];
252
253	0	let x2 = recip * recip;
254
255	0	let mut p = C[11];
256	0	for i in (0..11).rev() {
257	0	p = x2 * p + C[i];
258	0	}
259	0	p
260	0	}