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

Source
/*
 * // Copyright (c) Radzivon Bartoshyk 7/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::bessel::j0f_coeffs::{J0_ZEROS, J0_ZEROS_VALUE, J0F_COEFFS};
use crate::bessel::trigo_bessel::cos_small;
use crate::double_double::DoubleDouble;
use crate::polyeval::{f_polyeval9, f_polyeval10, f_polyeval12, f_polyeval14};
use crate::rounding::CpuCeil;
use crate::sincos_reduce::rem2pif_any;

/// Bessel of the first kind of order 0
///
/// Max ulp 0.5
pub fn f_j0f(x: f32) -> f32 {
    let ux = x.to_bits().wrapping_shl(1);
    if ux >= 0xffu32 << 24 || ux <= 0x6800_0000u32 {
        // |x| == 0, |x| == inf, |x| == NaN, |x| <= f32::EPSILON
        if ux == 0 {
            // |x| == 0
            return f64::from_bits(0x3ff0000000000000) as f32;
        }
        if x.is_infinite() {
            return 0.;
        }

        if ux <= 0x6800_0000u32 {
            // |x| < f32::EPSILON
            // taylor series for J0(x) ~ 1 - x^2/4 + O(x^4)
            #[cfg(any(
                all(
                    any(target_arch = "x86", target_arch = "x86_64"),
                    target_feature = "fma"
                ),
                target_arch = "aarch64"
            ))]
            {
                use crate::common::f_fmlaf;
                return f_fmlaf(x, -x * 0.25, 1.);
            }
            #[cfg(not(any(
                all(
                    any(target_arch = "x86", target_arch = "x86_64"),
                    target_feature = "fma"
                ),
                target_arch = "aarch64"
            )))]
            {
                use crate::common::f_fmla;
                let dx = x as f64;
                return f_fmla(dx, -dx * 0.25, 1.) as f32;
            }
        }

        return x + f32::NAN; // x == NaN
    }

    let x_abs = x.to_bits() & 0x7fff_ffff;

    if x_abs <= 0x4295999au32 {
        // |x| <= 74.8
        if x_abs <= 0x3e800000u32 {
            // |x| <= 0.25
            return j0f_maclaurin_series(x);
        }

        if x_abs == 0x401a42e8u32 {
            return f32::from_bits(0xbb3b2f69u32);
        }

        return small_argument_path(x);
    }

    // Exceptions
    if x_abs == 0x65ce46e4 {
        return f32::from_bits(0x1eed85c4);
    } else if x_abs == 0x7e3dcda0 {
        return f32::from_bits(0x92b81111);
    } else if x_abs == 0x76d84625 {
        return f32::from_bits(0x95d7a68b);
    } else if x_abs == 0x6bf68a7b {
        return f32::from_bits(0x1dc70a09);
    } else if x_abs == 0x7842c820 {
        return f32::from_bits(0x17ebf13e);
    } else if x_abs == 0x4ba332e9 {
        return f32::from_bits(0x27250206);
    }

    j0f_asympt(f32::from_bits(x_abs))
}

/**
Generated by SageMath:
```python
# Maclaurin series for j0
def print_expansion_at_0_f():
    print(f"pub(crate) const J0_MACLAURIN_SERIES: [u64; 9] = [")
    from mpmath import mp, j0, taylor
    mp.prec = 60
    poly = taylor(lambda val: j0(val), 0, 18)
    z = 0
    for i in range(0, 18, 2):
        print(f"{double_to_hex(poly[i])},")
    print("];")

    print(f"poly {poly}")

print_expansion_at_0_f()
```
**/
#[inline]
fn j0f_maclaurin_series(x: f32) -> f32 {
    pub(crate) const C: [u64; 9] = [
        0x3ff0000000000000,
        0xbfd0000000000000,
        0x3f90000000000000,
        0xbf3c71c71c71c71c,
        0x3edc71c71c71c71c,
        0xbe723456789abcdf,
        0x3e002e85c0898b71,
        0xbd8522a43f65486a,
        0x3d0522a43f65486a,
    ];
    let dx = x as f64;
    f_polyeval9(
        dx * dx,
        f64::from_bits(C[0]),
        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]),
    ) as f32
}

/// This method on small range searches for nearest zero or extremum.
/// Then picks stored series expansion at the point end evaluates the poly at the point.
#[inline]
fn small_argument_path(x: f32) -> f32 {
    let x_abs = f32::from_bits(x.to_bits() & 0x7fff_ffff) as f64;

    // let avg_step = 74.6145 / 47.0;
    // let inv_step = 1.0 / avg_step;

    const INV_STEP: f64 = 0.6299043751549631;

    let fx = x_abs * INV_STEP;
    const J0_ZEROS_COUNT: f64 = (J0_ZEROS.len() - 1) as f64;
    let idx0 = unsafe { fx.min(J0_ZEROS_COUNT).to_int_unchecked::<usize>() };
    let idx1 = unsafe {
        fx.cpu_ceil()
            .min(J0_ZEROS_COUNT)
            .to_int_unchecked::<usize>()
    };

    let found_zero0 = DoubleDouble::from_bit_pair(J0_ZEROS[idx0]);
    let found_zero1 = DoubleDouble::from_bit_pair(J0_ZEROS[idx1]);

    let dist0 = (found_zero0.hi - x_abs).abs();
    let dist1 = (found_zero1.hi - x_abs).abs();

    let (found_zero, idx, dist) = if dist0 < dist1 {
        (found_zero0, idx0, dist0)
    } else {
        (found_zero1, idx1, dist1)
    };

    if idx == 0 {
        return j0f_maclaurin_series(x);
    }

    // We hit exact zero, value, better to return it directly
    if dist == 0. {
        return f64::from_bits(J0_ZEROS_VALUE[idx]) as f32;
    }

    let c = &J0F_COEFFS[idx - 1];

    let r = (x_abs - found_zero.hi) - found_zero.lo;

    let p = f_polyeval14(
        r,
        f64::from_bits(c[0]),
        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]),
        f64::from_bits(c[10]),
        f64::from_bits(c[11]),
        f64::from_bits(c[12]),
        f64::from_bits(c[13]),
    );

    p as f32
}

#[inline]
pub(crate) fn j1f_rsqrt(x: f64) -> f64 {
    (1. / x) * x.sqrt()
}

/*
   Evaluates:
   J1 = sqrt(2/(PI*x)) * beta(x) * cos(x - PI/4 - alpha(x))
*/
#[inline]
fn j0f_asympt(x: f32) -> f32 {
    let dx = x as f64;

    let alpha = j0f_asympt_alpha(dx);
    let beta = j0f_asympt_beta(dx);

    let angle = rem2pif_any(x);

    const SQRT_2_OVER_PI: f64 = f64::from_bits(0x3fe9884533d43651);
    const MPI_OVER_4: f64 = f64::from_bits(0xbfe921fb54442d18);

    let x0pi34 = MPI_OVER_4 - alpha;
    let r0 = angle + x0pi34;

    let m_cos = cos_small(r0);

    let z0 = beta * m_cos;
    let scale = SQRT_2_OVER_PI * j1f_rsqrt(dx);

    (scale * z0) as f32
}

/**
Note expansion generation below: this is negative series expressed in Sage as positive,
so before any real evaluation `x=1/x` should be applied.

Generated by SageMath:
```python
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)

R_series = (-Q/P)

# alpha is atan(R_series) so we're doing Taylor series atan expansion on R_series

arctan_series_Z = sum([QQ(-1)**k * x**(QQ(2)*k+1) / RealField(700)(RealField(700)(2)*k+1) for k in range(25)])
alpha_series = arctan_series_Z(R_series)

# see the series
print(alpha_series)
```
**/
#[inline]
pub(crate) fn j0f_asympt_alpha(x: f64) -> f64 {
    const C: [u64; 12] = [
        0x3fc0000000000000,
        0xbfb0aaaaaaaaaaab,
        0x3fcad33333333333,
        0xbffa358492492492,
        0x403779a1f8e38e39,
        0xc080bd1fc8b1745d,
        0x40d16b51e66c789e,
        0xc128ecc3af33ab37,
        0x418779dae2b8512f,
        0xc1ec296336955c7f,
        0x4254f5ee683b6432,
        0xc2c2f51eced6693f,
    ];
    let recip = 1. / x;
    let x2 = recip * recip;
    let p = f_polyeval12(
        x2,
        f64::from_bits(C[0]),
        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]),
        f64::from_bits(C[10]),
        f64::from_bits(C[11]),
    );
    p * recip
}

/**
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 j0f_asympt_beta(x: f64) -> f64 {
    const C: [u64; 10] = [
        0x3ff0000000000000,
        0xbfb0000000000000,
        0x3fba800000000000,
        0xbfe15f0000000000,
        0x4017651180000000,
        0xc05ab8c13b800000,
        0x40a730492f262000,
        0xc0fc73a7acd696f0,
        0x41577458dd9fce68,
        0xc1b903ab9b27e18f,
    ];
    let recip = 1. / x;
    let x2 = recip * recip;
    f_polyeval10(
        x2,
        f64::from_bits(C[0]),
        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]),
    )
}

#[cfg(test)]
mod tests {
    use crate::f_j0f;

    #[test]
    fn test_j0f() {
        println!("0x{:8x}", f32::EPSILON.to_bits().wrapping_shl(1));
        assert_eq!(f_j0f(-3123.), 0.012329336);
        assert_eq!(f_j0f(-0.1), 0.99750155);
        assert_eq!(f_j0f(-15.1), -0.03456193);
        assert_eq!(f_j0f(3123.), 0.012329336);
        assert_eq!(f_j0f(0.1), 0.99750155);
        assert_eq!(f_j0f(15.1), -0.03456193);
        assert_eq!(f_j0f(f32::INFINITY), 0.);
        assert_eq!(f_j0f(f32::NEG_INFINITY), 0.);
        assert!(f_j0f(f32::NAN).is_nan());
    }
}

Coverage Report

Created: 2026-01-10 07:01

Line	Count	Source
1		/*
2		* // Copyright (c) Radzivon Bartoshyk 7/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::bessel::j0f_coeffs::{J0_ZEROS, J0_ZEROS_VALUE, J0F_COEFFS};
30		use crate::bessel::trigo_bessel::cos_small;
31		use crate::double_double::DoubleDouble;
32		use crate::polyeval::{f_polyeval9, f_polyeval10, f_polyeval12, f_polyeval14};
33		use crate::rounding::CpuCeil;
34		use crate::sincos_reduce::rem2pif_any;
35
36		/// Bessel of the first kind of order 0
37		///
38		/// Max ulp 0.5
39	0	pub fn f_j0f(x: f32) -> f32 {
40	0	let ux = x.to_bits().wrapping_shl(1);
41	0	if ux >= 0xffu32 << 24 \|\| ux <= 0x6800_0000u32 {
42		// \|x\| == 0, \|x\| == inf, \|x\| == NaN, \|x\| <= f32::EPSILON
43	0	if ux == 0 {
44		// \|x\| == 0
45	0	return f64::from_bits(0x3ff0000000000000) as f32;
46	0	}
47	0	if x.is_infinite() {
48	0	return 0.;
49	0	}
50
51	0	if ux <= 0x6800_0000u32 {
52		// \|x\| < f32::EPSILON
53		// taylor series for J0(x) ~ 1 - x^2/4 + O(x^4)
54		#[cfg(any(
55		all(
56		any(target_arch = "x86", target_arch = "x86_64"),
57		target_feature = "fma"
58		),
59		target_arch = "aarch64"
60		))]
61		{
62		use crate::common::f_fmlaf;
63		return f_fmlaf(x, -x * 0.25, 1.);
64		}
65		#[cfg(not(any(
66		all(
67		any(target_arch = "x86", target_arch = "x86_64"),
68		target_feature = "fma"
69		),
70		target_arch = "aarch64"
71		)))]
72		{
73		use crate::common::f_fmla;
74	0	let dx = x as f64;
75	0	return f_fmla(dx, -dx * 0.25, 1.) as f32;
76		}
77	0	}
78
79	0	return x + f32::NAN; // x == NaN
80	0	}
81
82	0	let x_abs = x.to_bits() & 0x7fff_ffff;
83
84	0	if x_abs <= 0x4295999au32 {
85		// \|x\| <= 74.8
86	0	if x_abs <= 0x3e800000u32 {
87		// \|x\| <= 0.25
88	0	return j0f_maclaurin_series(x);
89	0	}
90
91	0	if x_abs == 0x401a42e8u32 {
92	0	return f32::from_bits(0xbb3b2f69u32);
93	0	}
94
95	0	return small_argument_path(x);
96	0	}
97
98		// Exceptions
99	0	if x_abs == 0x65ce46e4 {
100	0	return f32::from_bits(0x1eed85c4);
101	0	} else if x_abs == 0x7e3dcda0 {
102	0	return f32::from_bits(0x92b81111);
103	0	} else if x_abs == 0x76d84625 {
104	0	return f32::from_bits(0x95d7a68b);
105	0	} else if x_abs == 0x6bf68a7b {
106	0	return f32::from_bits(0x1dc70a09);
107	0	} else if x_abs == 0x7842c820 {
108	0	return f32::from_bits(0x17ebf13e);
109	0	} else if x_abs == 0x4ba332e9 {
110	0	return f32::from_bits(0x27250206);
111	0	}
112
113	0	j0f_asympt(f32::from_bits(x_abs))
114	0	}
115
116		/**
117		Generated by SageMath:
118		```python
119		# Maclaurin series for j0
120		def print_expansion_at_0_f():
121		print(f"pub(crate) const J0_MACLAURIN_SERIES: [u64; 9] = [")
122		from mpmath import mp, j0, taylor
123		mp.prec = 60
124		poly = taylor(lambda val: j0(val), 0, 18)
125		z = 0
126		for i in range(0, 18, 2):
127		print(f"{double_to_hex(poly[i])},")
128		print("];")
129
130		print(f"poly {poly}")
131
132		print_expansion_at_0_f()
133		```
134		**/
135		#[inline]
136	0	fn j0f_maclaurin_series(x: f32) -> f32 {
137		pub(crate) const C: [u64; 9] = [
138		0x3ff0000000000000,
139		0xbfd0000000000000,
140		0x3f90000000000000,
141		0xbf3c71c71c71c71c,
142		0x3edc71c71c71c71c,
143		0xbe723456789abcdf,
144		0x3e002e85c0898b71,
145		0xbd8522a43f65486a,
146		0x3d0522a43f65486a,
147		];
148	0	let dx = x as f64;
149	0	f_polyeval9(
150	0	dx * dx,
151	0	f64::from_bits(C[0]),
152	0	f64::from_bits(C[1]),
153	0	f64::from_bits(C[2]),
154	0	f64::from_bits(C[3]),
155	0	f64::from_bits(C[4]),
156	0	f64::from_bits(C[5]),
157	0	f64::from_bits(C[6]),
158	0	f64::from_bits(C[7]),
159	0	f64::from_bits(C[8]),
160	0	) as f32
161	0	}
162
163		/// This method on small range searches for nearest zero or extremum.
164		/// Then picks stored series expansion at the point end evaluates the poly at the point.
165		#[inline]
166	0	fn small_argument_path(x: f32) -> f32 {
167	0	let x_abs = f32::from_bits(x.to_bits() & 0x7fff_ffff) as f64;
168
169		// let avg_step = 74.6145 / 47.0;
170		// let inv_step = 1.0 / avg_step;
171
172		const INV_STEP: f64 = 0.6299043751549631;
173
174	0	let fx = x_abs * INV_STEP;
175		const J0_ZEROS_COUNT: f64 = (J0_ZEROS.len() - 1) as f64;
176	0	let idx0 = unsafe { fx.min(J0_ZEROS_COUNT).to_int_unchecked::<usize>() };
177	0	let idx1 = unsafe {
178	0	fx.cpu_ceil()
179	0	.min(J0_ZEROS_COUNT)
180	0	.to_int_unchecked::<usize>()
181		};
182
183	0	let found_zero0 = DoubleDouble::from_bit_pair(J0_ZEROS[idx0]);
184	0	let found_zero1 = DoubleDouble::from_bit_pair(J0_ZEROS[idx1]);
185
186	0	let dist0 = (found_zero0.hi - x_abs).abs();
187	0	let dist1 = (found_zero1.hi - x_abs).abs();
188
189	0	let (found_zero, idx, dist) = if dist0 < dist1 {
190	0	(found_zero0, idx0, dist0)
191		} else {
192	0	(found_zero1, idx1, dist1)
193		};
194
195	0	if idx == 0 {
196	0	return j0f_maclaurin_series(x);
197	0	}
198
199		// We hit exact zero, value, better to return it directly
200	0	if dist == 0. {
201	0	return f64::from_bits(J0_ZEROS_VALUE[idx]) as f32;
202	0	}
203
204	0	let c = &J0F_COEFFS[idx - 1];
205
206	0	let r = (x_abs - found_zero.hi) - found_zero.lo;
207
208	0	let p = f_polyeval14(
209	0	r,
210	0	f64::from_bits(c[0]),
211	0	f64::from_bits(c[1]),
212	0	f64::from_bits(c[2]),
213	0	f64::from_bits(c[3]),
214	0	f64::from_bits(c[4]),
215	0	f64::from_bits(c[5]),
216	0	f64::from_bits(c[6]),
217	0	f64::from_bits(c[7]),
218	0	f64::from_bits(c[8]),
219	0	f64::from_bits(c[9]),
220	0	f64::from_bits(c[10]),
221	0	f64::from_bits(c[11]),
222	0	f64::from_bits(c[12]),
223	0	f64::from_bits(c[13]),
224		);
225
226	0	p as f32
227	0	}
228
229		#[inline]
230	0	pub(crate) fn j1f_rsqrt(x: f64) -> f64 {
231	0	(1. / x) * x.sqrt()
232	0	}
233
234		/*
235		Evaluates:
236		J1 = sqrt(2/(PIx)) beta(x) * cos(x - PI/4 - alpha(x))
237		*/
238		#[inline]
239	0	fn j0f_asympt(x: f32) -> f32 {
240	0	let dx = x as f64;
241
242	0	let alpha = j0f_asympt_alpha(dx);
243	0	let beta = j0f_asympt_beta(dx);
244
245	0	let angle = rem2pif_any(x);
246
247		const SQRT_2_OVER_PI: f64 = f64::from_bits(0x3fe9884533d43651);
248		const MPI_OVER_4: f64 = f64::from_bits(0xbfe921fb54442d18);
249
250	0	let x0pi34 = MPI_OVER_4 - alpha;
251	0	let r0 = angle + x0pi34;
252
253	0	let m_cos = cos_small(r0);
254
255	0	let z0 = beta * m_cos;
256	0	let scale = SQRT_2_OVER_PI * j1f_rsqrt(dx);
257
258	0	(scale * z0) as f32
259	0	}
260
261		/**
262		Note expansion generation below: this is negative series expressed in Sage as positive,
263		so before any real evaluation `x=1/x` should be applied.
264
265		Generated by SageMath:
266		```python
267		def binomial_like(n, m):
268		prod = QQ(1)
269		z = QQ(4)(n*2)
270		for k in range(1,m + 1):
271		prod = (z - (2k - 1)**2)
272		return prod / (QQ(2)*(2m) * (ZZ(m).factorial()))
273
274		R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
275		x = R.gen()
276
277		def Pn_asymptotic(n, y, terms=10):
278		# now y = 1/x
279		return sum( (-1)*m binomial_like(n, 2m) / (QQ(2)(2m)) * y*(QQ(2)m) for m in range(terms) )
280
281		def Qn_asymptotic(n, y, terms=10):
282		return sum( (-1)*m binomial_like(n, 2m + 1) / (QQ(2)(2m + 1)) * y*(QQ(2)m + 1) for m in range(terms) )
283
284		P = Pn_asymptotic(0, x, 50)
285		Q = Qn_asymptotic(0, x, 50)
286
287		R_series = (-Q/P)
288
289		# alpha is atan(R_series) so we're doing Taylor series atan expansion on R_series
290
291		arctan_series_Z = sum([QQ(-1)*k x*(QQ(2)k+1) / RealField(700)(RealField(700)(2)*k+1) for k in range(25)])
292		alpha_series = arctan_series_Z(R_series)
293
294		# see the series
295		print(alpha_series)
296		```
297		**/
298		#[inline]
299	0	pub(crate) fn j0f_asympt_alpha(x: f64) -> f64 {
300		const C: [u64; 12] = [
301		0x3fc0000000000000,
302		0xbfb0aaaaaaaaaaab,
303		0x3fcad33333333333,
304		0xbffa358492492492,
305		0x403779a1f8e38e39,
306		0xc080bd1fc8b1745d,
307		0x40d16b51e66c789e,
308		0xc128ecc3af33ab37,
309		0x418779dae2b8512f,
310		0xc1ec296336955c7f,
311		0x4254f5ee683b6432,
312		0xc2c2f51eced6693f,
313		];
314	0	let recip = 1. / x;
315	0	let x2 = recip * recip;
316	0	let p = f_polyeval12(
317	0	x2,
318	0	f64::from_bits(C[0]),
319	0	f64::from_bits(C[1]),
320	0	f64::from_bits(C[2]),
321	0	f64::from_bits(C[3]),
322	0	f64::from_bits(C[4]),
323	0	f64::from_bits(C[5]),
324	0	f64::from_bits(C[6]),
325	0	f64::from_bits(C[7]),
326	0	f64::from_bits(C[8]),
327	0	f64::from_bits(C[9]),
328	0	f64::from_bits(C[10]),
329	0	f64::from_bits(C[11]),
330		);
331	0	p * recip
332	0	}
333
334		/**
335		Beta series
336
337		Generated by SageMath:
338		```python
339		#generate b series
340		def binomial_like(n, m):
341		prod = QQ(1)
342		z = QQ(4)(n*2)
343		for k in range(1,m + 1):
344		prod = (z - (2k - 1)**2)
345		return prod / (QQ(2)*(2m) * (ZZ(m).factorial()))
346
347		R = LaurentSeriesRing(RealField(300), 'x', default_prec=300)
348		x = R.gen()
349
350		def Pn_asymptotic(n, y, terms=10):
351		# now y = 1/x
352		return sum( (-1)*m binomial_like(n, 2m) / (QQ(2)(2m)) * y*(QQ(2)m) for m in range(terms) )
353
354		def Qn_asymptotic(n, y, terms=10):
355		return sum( (-1)*m binomial_like(n, 2m + 1) / (QQ(2)(2m + 1)) * y*(QQ(2)m + 1) for m in range(terms) )
356
357		P = Pn_asymptotic(0, x, 50)
358		Q = Qn_asymptotic(0, x, 50)
359
360		def sqrt_series(s):
361		val = S.valuation()
362		lc = S[val] # Leading coefficient
363		b = lc.sqrt() * x**(val // 2)
364
365		for _ in range(5):
366		b = (b + S / b) / 2
367		b = b
368		return b
369
370		S = (P2 + Q2).truncate(50)
371
372		b_series = sqrt_series(S).truncate(30)
373		#see the series
374		print(b_series)
375		```
376		**/
377		#[inline]
378	0	pub(crate) fn j0f_asympt_beta(x: f64) -> f64 {
379		const C: [u64; 10] = [
380		0x3ff0000000000000,
381		0xbfb0000000000000,
382		0x3fba800000000000,
383		0xbfe15f0000000000,
384		0x4017651180000000,
385		0xc05ab8c13b800000,
386		0x40a730492f262000,
387		0xc0fc73a7acd696f0,
388		0x41577458dd9fce68,
389		0xc1b903ab9b27e18f,
390		];
391	0	let recip = 1. / x;
392	0	let x2 = recip * recip;
393	0	f_polyeval10(
394	0	x2,
395	0	f64::from_bits(C[0]),
396	0	f64::from_bits(C[1]),
397	0	f64::from_bits(C[2]),
398	0	f64::from_bits(C[3]),
399	0	f64::from_bits(C[4]),
400	0	f64::from_bits(C[5]),
401	0	f64::from_bits(C[6]),
402	0	f64::from_bits(C[7]),
403	0	f64::from_bits(C[8]),
404	0	f64::from_bits(C[9]),
405		)
406	0	}
407
408		#[cfg(test)]
409		mod tests {
410		use crate::f_j0f;
411
412		#[test]
413		fn test_j0f() {
414		println!("0x{:8x}", f32::EPSILON.to_bits().wrapping_shl(1));
415		assert_eq!(f_j0f(-3123.), 0.012329336);
416		assert_eq!(f_j0f(-0.1), 0.99750155);
417		assert_eq!(f_j0f(-15.1), -0.03456193);
418		assert_eq!(f_j0f(3123.), 0.012329336);
419		assert_eq!(f_j0f(0.1), 0.99750155);
420		assert_eq!(f_j0f(15.1), -0.03456193);
421		assert_eq!(f_j0f(f32::INFINITY), 0.);
422		assert_eq!(f_j0f(f32::NEG_INFINITY), 0.);
423		assert!(f_j0f(f32::NAN).is_nan());
424		}
425		}