/rust/registry/src/index.crates.io-1949cf8c6b5b557f/pxfm-0.1.27/src/bessel/j1f.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::j1_coeffs::{J1_ZEROS, J1_ZEROS_VALUE};
use crate::bessel::j1f_coeffs::J1F_COEFFS;
use crate::bessel::trigo_bessel::sin_small;
use crate::double_double::DoubleDouble;
use crate::polyeval::{f_polyeval7, f_polyeval10, f_polyeval12, f_polyeval14};
use crate::rounding::CpuCeil;
use crate::sincos_reduce::rem2pif_any;

/// Bessel of the first kind of order 1
///
/// Max ULP 0.5
///
/// Note about accuracy:
/// - Close to zero Bessel have tiny values such that testing against MPFR must be done exactly
///   in the same precision, since any nearest representable number have ULP > 0.5.
///   For example `J1(0.000000000000000000000000000000000000023509886)` in single precision
///   have an error at least 0.72 ULP for any number with extended precision,
///   that would be represented in f32.
pub fn f_j1f(x: f32) -> f32 {
    let ux = x.to_bits().wrapping_shl(1);
    if ux >= 0xffu32 << 24 || ux == 0 {
        // |x| == 0, |x| == inf, |x| == NaN
        if ux == 0 {
            // |x| == 0
            return x;
        }
        if x.is_infinite() {
            return 0.;
        }
        return x + f32::NAN; // x == NaN
    }

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

    if ax < 0x429533c2u32 {
        // |x| <= 74.60109
        if ax < 0x3e800000u32 {
            // |x| <= 0.25
            if ax <= 0x34000000u32 {
                // |x| <= f32::EPSILON
                // taylor series for J1(x) ~ x/2 + O(x^3)
                return x * 0.5;
            }
            return poly_near_zero(x);
        }
        return small_argument_path(x);
    }

    // Exceptional cases:
    if ax == 0x6ef9be45 {
        return if x.is_sign_negative() {
            f32::from_bits(0x187d8a8f)
        } else {
            -f32::from_bits(0x187d8a8f)
        };
    } else if ax == 0x7f0e5a38 {
        return if x.is_sign_negative() {
            -f32::from_bits(0x131f680b)
        } else {
            f32::from_bits(0x131f680b)
        };
    }

    j1f_asympt(x) as f32
}

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

/*
   Evaluates:
   J1 = sqrt(2/(PI*x)) * beta(x) * cos(x - 3*PI/4 - alpha(x))
   discarding 1*PI/2 using identities gives:
   J1 = sqrt(2/(PI*x)) * beta(x) * sin(x - PI/4 - alpha(x))

   to avoid squashing small (-PI/4 - alpha(x)) into a large x actual expansion is:

   J1 = sqrt(2/(PI*x)) * beta(x) * sin((x mod 2*PI) - PI/4 - alpha(x))
*/
#[inline]
fn j1f_asympt(x: f32) -> f64 {
    static SGN: [f64; 2] = [1., -1.];
    let sign_scale = SGN[x.is_sign_negative() as usize];
    let x = f32::from_bits(x.to_bits() & 0x7fff_ffff);

    let dx = x as f64;

    let alpha = j1f_asympt_alpha(dx);
    let beta = j1f_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_sin = sin_small(r0);

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

    scale * z0 * sign_scale
}

/**
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(1, x, 50)
Q = Qn_asymptotic(1, 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)
```

See notes/bessel_asympt.ipynb for generation
**/
#[inline]
pub(crate) fn j1f_asympt_alpha(x: f64) -> f64 {
    const C: [u64; 12] = [
        0xbfd8000000000000,
        0x3fc5000000000000,
        0xbfd7bccccccccccd,
        0x4002f486db6db6db,
        0xc03e9fbf40000000,
        0x4084997b55945d17,
        0xc0d4a914195269d9,
        0x412cd1b53816aec1,
        0xc18aa4095d419351,
        0x41ef809305f11b9d,
        0xc2572e6809ed618b,
        0x42c4c5b6057839f9,
    ];
    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
}

/**
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(1, x, 50)
Q = Qn_asymptotic(1, 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 beta series
print(b_series)
```

See notes/bessel_asympt.ipynb for generation
**/
#[inline]
pub(crate) fn j1f_asympt_beta(x: f64) -> f64 {
    const C: [u64; 10] = [
        0x3ff0000000000000,
        0x3fc8000000000000,
        0xbfc8c00000000000,
        0x3fe9c50000000000,
        0xc01ef5b680000000,
        0x40609860dd400000,
        0xc0abae9b7a06e000,
        0x41008711d41c1428,
        0xc15ab70164c8be6e,
        0x41bc1055e24f297f,
    ];
    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]),
    )
}

/**
Generated in Sollya:
```python
pretty = proc(u) {
  return ~(floor(u*1000)/1000);
};

bessel_j1 = library("./cmake-build-release/libbessel_sollya.dylib");

f = bessel_j1(x)/x;
d = [0, 0.921];
w = 1;
pf = fpminimax(f, [|0,2,4,6,8,10,12|], [|D...|], d, absolute, floating);

w = 1;
or_f = bessel_j1(x);
pf1 = pf * x;
err_p = -log2(dirtyinfnorm(pf1*w-or_f, d));
print ("relative error:", pretty(err_p));

for i from 0 to degree(pf) by 2 do {
    print("'", coeff(pf, i), "',");
};
```
See ./notes/bessel_sollya/bessel_j1f_at_zero.sollya
**/
#[inline]
fn poly_near_zero(x: f32) -> f32 {
    let dx = x as f64;
    let x2 = dx * dx;
    let p = f_polyeval7(
        x2,
        f64::from_bits(0x3fe0000000000000),
        f64::from_bits(0xbfaffffffffffffc),
        f64::from_bits(0x3f65555555554089),
        f64::from_bits(0xbf0c71c71c2a74ae),
        f64::from_bits(0x3ea6c16bbd1dc5c1),
        f64::from_bits(0xbe384562afb69e7d),
        f64::from_bits(0x3dc248d0d0221cd0),
    );
    (p * dx) 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 {
    static SIGN: [f64; 2] = [1., -1.];
    let sign_scale = SIGN[x.is_sign_negative() as usize];
    let x_abs = f32::from_bits(x.to_bits() & 0x7fff_ffff) as f64;

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

    const INV_STEP: f64 = 0.6300176043004198;

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

    let found_zero0 = DoubleDouble::from_bit_pair(J1_ZEROS[idx0]);
    let found_zero1 = DoubleDouble::from_bit_pair(J1_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 poly_near_zero(x);
    }

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

    let c = &J1F_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 * sign_scale) as f32
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_f_j1f() {
        assert_eq!(
            f_j1f(77.743162408196766932633181568235159),
            0.09049267898021947
        );
        assert_eq!(
            f_j1f(-0.000000000000000000000000000000000000008827127),
            -0.0000000000000000000000000000000000000044135635
        );
        assert_eq!(
            f_j1f(0.000000000000000000000000000000000000008827127),
            0.0000000000000000000000000000000000000044135635
        );
        assert_eq!(f_j1f(5.4), -0.3453447907795863);
        assert_eq!(
            f_j1f(84.027189586293545175976760219782591),
            0.0870430264022591
        );
        assert_eq!(f_j1f(f32::INFINITY), 0.);
        assert_eq!(f_j1f(f32::NEG_INFINITY), 0.);
        assert!(f_j1f(f32::NAN).is_nan());
        assert_eq!(f_j1f(-1.7014118e38), 0.000000000000000000006856925);
    }
}

Coverage Report

Created: 2025-12-11 07:11

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::j1_coeffs::{J1_ZEROS, J1_ZEROS_VALUE};
30		use crate::bessel::j1f_coeffs::J1F_COEFFS;
31		use crate::bessel::trigo_bessel::sin_small;
32		use crate::double_double::DoubleDouble;
33		use crate::polyeval::{f_polyeval7, f_polyeval10, f_polyeval12, f_polyeval14};
34		use crate::rounding::CpuCeil;
35		use crate::sincos_reduce::rem2pif_any;
36
37		/// Bessel of the first kind of order 1
38		///
39		/// Max ULP 0.5
40		///
41		/// Note about accuracy:
42		/// - Close to zero Bessel have tiny values such that testing against MPFR must be done exactly
43		/// in the same precision, since any nearest representable number have ULP > 0.5.
44		/// For example `J1(0.000000000000000000000000000000000000023509886)` in single precision
45		/// have an error at least 0.72 ULP for any number with extended precision,
46		/// that would be represented in f32.
47	0	pub fn f_j1f(x: f32) -> f32 {
48	0	let ux = x.to_bits().wrapping_shl(1);
49	0	if ux >= 0xffu32 << 24 \|\| ux == 0 {
50		// \|x\| == 0, \|x\| == inf, \|x\| == NaN
51	0	if ux == 0 {
52		// \|x\| == 0
53	0	return x;
54	0	}
55	0	if x.is_infinite() {
56	0	return 0.;
57	0	}
58	0	return x + f32::NAN; // x == NaN
59	0	}
60
61	0	let ax = x.to_bits() & 0x7fff_ffff;
62
63	0	if ax < 0x429533c2u32 {
64		// \|x\| <= 74.60109
65	0	if ax < 0x3e800000u32 {
66		// \|x\| <= 0.25
67	0	if ax <= 0x34000000u32 {
68		// \|x\| <= f32::EPSILON
69		// taylor series for J1(x) ~ x/2 + O(x^3)
70	0	return x * 0.5;
71	0	}
72	0	return poly_near_zero(x);
73	0	}
74	0	return small_argument_path(x);
75	0	}
76
77		// Exceptional cases:
78	0	if ax == 0x6ef9be45 {
79	0	return if x.is_sign_negative() {
80	0	f32::from_bits(0x187d8a8f)
81		} else {
82	0	-f32::from_bits(0x187d8a8f)
83		};
84	0	} else if ax == 0x7f0e5a38 {
85	0	return if x.is_sign_negative() {
86	0	-f32::from_bits(0x131f680b)
87		} else {
88	0	f32::from_bits(0x131f680b)
89		};
90	0	}
91
92	0	j1f_asympt(x) as f32
93	0	}
94
95		#[inline]
96	0	fn j1f_rsqrt(x: f64) -> f64 {
97	0	(1. / x) * x.sqrt()
98	0	}
99
100		/*
101		Evaluates:
102		J1 = sqrt(2/(PIx)) beta(x) * cos(x - 3*PI/4 - alpha(x))
103		discarding 1*PI/2 using identities gives:
104		J1 = sqrt(2/(PIx)) beta(x) * sin(x - PI/4 - alpha(x))
105
106		to avoid squashing small (-PI/4 - alpha(x)) into a large x actual expansion is:
107
108		J1 = sqrt(2/(PIx)) beta(x) * sin((x mod 2*PI) - PI/4 - alpha(x))
109		*/
110		#[inline]
111	0	fn j1f_asympt(x: f32) -> f64 {
112		static SGN: [f64; 2] = [1., -1.];
113	0	let sign_scale = SGN[x.is_sign_negative() as usize];
114	0	let x = f32::from_bits(x.to_bits() & 0x7fff_ffff);
115
116	0	let dx = x as f64;
117
118	0	let alpha = j1f_asympt_alpha(dx);
119	0	let beta = j1f_asympt_beta(dx);
120
121	0	let angle = rem2pif_any(x);
122
123		const SQRT_2_OVER_PI: f64 = f64::from_bits(0x3fe9884533d43651);
124		const MPI_OVER_4: f64 = f64::from_bits(0xbfe921fb54442d18);
125
126	0	let x0pi34 = MPI_OVER_4 - alpha;
127	0	let r0 = angle + x0pi34;
128
129	0	let m_sin = sin_small(r0);
130
131	0	let z0 = beta * m_sin;
132	0	let scale = SQRT_2_OVER_PI * j1f_rsqrt(dx);
133
134	0	scale * z0 * sign_scale
135	0	}
136
137		/**
138		Note expansion generation below: this is negative series expressed in Sage as positive,
139		so before any real evaluation `x=1/x` should be applied.
140
141		Generated by SageMath:
142		```python
143		def binomial_like(n, m):
144		prod = QQ(1)
145		z = QQ(4)(n*2)
146		for k in range(1,m + 1):
147		prod = (z - (2k - 1)**2)
148		return prod / (QQ(2)*(2m) * (ZZ(m).factorial()))
149
150		R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
151		x = R.gen()
152
153		def Pn_asymptotic(n, y, terms=10):
154		# now y = 1/x
155		return sum( (-1)*m binomial_like(n, 2m) / (QQ(2)(2m)) * y*(QQ(2)m) for m in range(terms) )
156
157		def Qn_asymptotic(n, y, terms=10):
158		return sum( (-1)*m binomial_like(n, 2m + 1) / (QQ(2)(2m + 1)) * y*(QQ(2)m + 1) for m in range(terms) )
159
160		P = Pn_asymptotic(1, x, 50)
161		Q = Qn_asymptotic(1, x, 50)
162
163		R_series = (-Q/P)
164
165		# alpha is atan(R_series) so we're doing Taylor series atan expansion on R_series
166
167		arctan_series_Z = sum([QQ(-1)*k x*(QQ(2)k+1) / RealField(700)(RealField(700)(2)*k+1) for k in range(25)])
168		alpha_series = arctan_series_Z(R_series)
169
170		# see the series
171		print(alpha_series)
172		```
173
174		See notes/bessel_asympt.ipynb for generation
175		**/
176		#[inline]
177	0	pub(crate) fn j1f_asympt_alpha(x: f64) -> f64 {
178		const C: [u64; 12] = [
179		0xbfd8000000000000,
180		0x3fc5000000000000,
181		0xbfd7bccccccccccd,
182		0x4002f486db6db6db,
183		0xc03e9fbf40000000,
184		0x4084997b55945d17,
185		0xc0d4a914195269d9,
186		0x412cd1b53816aec1,
187		0xc18aa4095d419351,
188		0x41ef809305f11b9d,
189		0xc2572e6809ed618b,
190		0x42c4c5b6057839f9,
191		];
192	0	let recip = 1. / x;
193	0	let x2 = recip * recip;
194	0	let p = f_polyeval12(
195	0	x2,
196	0	f64::from_bits(C[0]),
197	0	f64::from_bits(C[1]),
198	0	f64::from_bits(C[2]),
199	0	f64::from_bits(C[3]),
200	0	f64::from_bits(C[4]),
201	0	f64::from_bits(C[5]),
202	0	f64::from_bits(C[6]),
203	0	f64::from_bits(C[7]),
204	0	f64::from_bits(C[8]),
205	0	f64::from_bits(C[9]),
206	0	f64::from_bits(C[10]),
207	0	f64::from_bits(C[11]),
208		);
209	0	p * recip
210	0	}
211
212		/**
213		Note expansion generation below: this is negative series expressed in Sage as positive,
214		so before any real evaluation `x=1/x` should be applied
215
216		Generated by SageMath:
217		```python
218		def binomial_like(n, m):
219		prod = QQ(1)
220		z = QQ(4)(n*2)
221		for k in range(1,m + 1):
222		prod = (z - (2k - 1)**2)
223		return prod / (QQ(2)*(2m) * (ZZ(m).factorial()))
224
225		R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
226		x = R.gen()
227
228		def Pn_asymptotic(n, y, terms=10):
229		# now y = 1/x
230		return sum( (-1)*m binomial_like(n, 2m) / (QQ(2)(2m)) * y*(QQ(2)m) for m in range(terms) )
231
232		def Qn_asymptotic(n, y, terms=10):
233		return sum( (-1)*m binomial_like(n, 2m + 1) / (QQ(2)(2m + 1)) * y*(QQ(2)m + 1) for m in range(terms) )
234
235		P = Pn_asymptotic(1, x, 50)
236		Q = Qn_asymptotic(1, x, 50)
237
238		def sqrt_series(s):
239		val = S.valuation()
240		lc = S[val] # Leading coefficient
241		b = lc.sqrt() * x**(val // 2)
242
243		for _ in range(5):
244		b = (b + S / b) / 2
245		b = b
246		return b
247
248		S = (P2 + Q2).truncate(50)
249
250		b_series = sqrt_series(S).truncate(30)
251		# see the beta series
252		print(b_series)
253		```
254
255		See notes/bessel_asympt.ipynb for generation
256		**/
257		#[inline]
258	0	pub(crate) fn j1f_asympt_beta(x: f64) -> f64 {
259		const C: [u64; 10] = [
260		0x3ff0000000000000,
261		0x3fc8000000000000,
262		0xbfc8c00000000000,
263		0x3fe9c50000000000,
264		0xc01ef5b680000000,
265		0x40609860dd400000,
266		0xc0abae9b7a06e000,
267		0x41008711d41c1428,
268		0xc15ab70164c8be6e,
269		0x41bc1055e24f297f,
270		];
271	0	let recip = 1. / x;
272	0	let x2 = recip * recip;
273	0	f_polyeval10(
274	0	x2,
275	0	f64::from_bits(C[0]),
276	0	f64::from_bits(C[1]),
277	0	f64::from_bits(C[2]),
278	0	f64::from_bits(C[3]),
279	0	f64::from_bits(C[4]),
280	0	f64::from_bits(C[5]),
281	0	f64::from_bits(C[6]),
282	0	f64::from_bits(C[7]),
283	0	f64::from_bits(C[8]),
284	0	f64::from_bits(C[9]),
285		)
286	0	}
287
288		/**
289		Generated in Sollya:
290		```python
291		pretty = proc(u) {
292		return ~(floor(u*1000)/1000);
293		};
294
295		bessel_j1 = library("./cmake-build-release/libbessel_sollya.dylib");
296
297		f = bessel_j1(x)/x;
298		d = [0, 0.921];
299		w = 1;
300		pf = fpminimax(f, [\|0,2,4,6,8,10,12\|], [\|D...\|], d, absolute, floating);
301
302		w = 1;
303		or_f = bessel_j1(x);
304		pf1 = pf * x;
305		err_p = -log2(dirtyinfnorm(pf1*w-or_f, d));
306		print ("relative error:", pretty(err_p));
307
308		for i from 0 to degree(pf) by 2 do {
309		print("'", coeff(pf, i), "',");
310		};
311		```
312		See ./notes/bessel_sollya/bessel_j1f_at_zero.sollya
313		**/
314		#[inline]
315	0	fn poly_near_zero(x: f32) -> f32 {
316	0	let dx = x as f64;
317	0	let x2 = dx * dx;
318	0	let p = f_polyeval7(
319	0	x2,
320	0	f64::from_bits(0x3fe0000000000000),
321	0	f64::from_bits(0xbfaffffffffffffc),
322	0	f64::from_bits(0x3f65555555554089),
323	0	f64::from_bits(0xbf0c71c71c2a74ae),
324	0	f64::from_bits(0x3ea6c16bbd1dc5c1),
325	0	f64::from_bits(0xbe384562afb69e7d),
326	0	f64::from_bits(0x3dc248d0d0221cd0),
327		);
328	0	(p * dx) as f32
329	0	}
330
331		/// This method on small range searches for nearest zero or extremum.
332		/// Then picks stored series expansion at the point end evaluates the poly at the point.
333		#[inline]
334	0	fn small_argument_path(x: f32) -> f32 {
335		static SIGN: [f64; 2] = [1., -1.];
336	0	let sign_scale = SIGN[x.is_sign_negative() as usize];
337	0	let x_abs = f32::from_bits(x.to_bits() & 0x7fff_ffff) as f64;
338
339		// let avg_step = 74.60109 / 47.0;
340		// let inv_step = 1.0 / avg_step;
341
342		const INV_STEP: f64 = 0.6300176043004198;
343
344	0	let fx = x_abs * INV_STEP;
345		const J1_ZEROS_COUNT: f64 = (J1_ZEROS.len() - 1) as f64;
346	0	let idx0 = unsafe { fx.min(J1_ZEROS_COUNT).to_int_unchecked::<usize>() };
347	0	let idx1 = unsafe {
348	0	fx.cpu_ceil()
349	0	.min(J1_ZEROS_COUNT)
350	0	.to_int_unchecked::<usize>()
351		};
352
353	0	let found_zero0 = DoubleDouble::from_bit_pair(J1_ZEROS[idx0]);
354	0	let found_zero1 = DoubleDouble::from_bit_pair(J1_ZEROS[idx1]);
355
356	0	let dist0 = (found_zero0.hi - x_abs).abs();
357	0	let dist1 = (found_zero1.hi - x_abs).abs();
358
359	0	let (found_zero, idx, dist) = if dist0 < dist1 {
360	0	(found_zero0, idx0, dist0)
361		} else {
362	0	(found_zero1, idx1, dist1)
363		};
364
365	0	if idx == 0 {
366	0	return poly_near_zero(x);
367	0	}
368
369		// We hit exact zero, value, better to return it directly
370	0	if dist == 0. {
371	0	return (f64::from_bits(J1_ZEROS_VALUE[idx]) * sign_scale) as f32;
372	0	}
373
374	0	let c = &J1F_COEFFS[idx - 1];
375
376	0	let r = (x_abs - found_zero.hi) - found_zero.lo;
377
378	0	let p = f_polyeval14(
379	0	r,
380	0	f64::from_bits(c[0]),
381	0	f64::from_bits(c[1]),
382	0	f64::from_bits(c[2]),
383	0	f64::from_bits(c[3]),
384	0	f64::from_bits(c[4]),
385	0	f64::from_bits(c[5]),
386	0	f64::from_bits(c[6]),
387	0	f64::from_bits(c[7]),
388	0	f64::from_bits(c[8]),
389	0	f64::from_bits(c[9]),
390	0	f64::from_bits(c[10]),
391	0	f64::from_bits(c[11]),
392	0	f64::from_bits(c[12]),
393	0	f64::from_bits(c[13]),
394		);
395
396	0	(p * sign_scale) as f32
397	0	}
398
399		#[cfg(test)]
400		mod tests {
401		use super::*;
402
403		#[test]
404		fn test_f_j1f() {
405		assert_eq!(
406		f_j1f(77.743162408196766932633181568235159),
407		0.09049267898021947
408		);
409		assert_eq!(
410		f_j1f(-0.000000000000000000000000000000000000008827127),
411		-0.0000000000000000000000000000000000000044135635
412		);
413		assert_eq!(
414		f_j1f(0.000000000000000000000000000000000000008827127),
415		0.0000000000000000000000000000000000000044135635
416		);
417		assert_eq!(f_j1f(5.4), -0.3453447907795863);
418		assert_eq!(
419		f_j1f(84.027189586293545175976760219782591),
420		0.0870430264022591
421		);
422		assert_eq!(f_j1f(f32::INFINITY), 0.);
423		assert_eq!(f_j1f(f32::NEG_INFINITY), 0.);
424		assert!(f_j1f(f32::NAN).is_nan());
425		assert_eq!(f_j1f(-1.7014118e38), 0.000000000000000000006856925);
426		}
427		}