pxfm/bessel/
k1ef.rs

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 */
29use crate::bessel::j0f::j1f_rsqrt;
30use crate::common::f_fmla;
31use crate::exponents::core_expf;
32use crate::logs::fast_logf;
33use crate::polyeval::{f_estrin_polyeval8, f_polyeval3, f_polyeval4};
34
35/// Modified exponentially scaled Bessel of the second kind of order 1
36///
37/// Computes K1(x)exp(x)
38///
39/// Max ULP 0.5
40pub fn f_k1ef(x: f32) -> f32 {
41    let ux = x.to_bits();
42    if ux >= 0xffu32 << 23 || ux == 0 {
43        // |x| == 0, |x| == inf, |x| == NaN, x < 0
44        if ux.wrapping_shl(1) == 0 {
45            return f32::INFINITY;
46        }
47        if x.is_infinite() {
48            return if x.is_sign_positive() { 0. } else { f32::NAN };
49        }
50        return x + f32::NAN; // x == NaN
51    }
52
53    let xb = x.to_bits();
54
55    if xb <= 0x3f800000u32 {
56        // x <= 1.0
57        if xb <= 0x34000000u32 {
58            // |x| <= f32::EPSILON
59            let dx = x as f64;
60            let leading_term = 1. / dx + 1.;
61            if xb <= 0x3109705fu32 {
62                // |x| <= 2e-9
63                // taylor series for tiny K1(x)exp(x) ~ 1/x + 1 + O(x)
64                return leading_term as f32;
65            }
66            // taylor series for small K1(x)exp(x) ~ 1/x+1+1/4 (1+2 EulerGamma-2 Log[2]+2 Log[x]) x + O(x^3)
67            const C: f64 = f64::from_bits(0xbffd8773039049e8); // 1 + 2 EulerGamma-2 Log[2]
68            let log_x = fast_logf(x);
69            let r = f_fmla(log_x, 2., C);
70            let w0 = f_fmla(dx * 0.25, r, leading_term);
71            return w0 as f32;
72        }
73        return k1ef_small(x);
74    }
75
76    k1ef_asympt(x)
77}
78
79/**
80Computes
81I1(x) = x/2 * (1 + 1 * (x/2)^2 + (x/2)^4 * P((x/2)^2))
82
83Generated by Woflram Mathematica:
84
85```text
86<<FunctionApproximations`
87ClearAll["Global`*"]
88f[x_]:=(BesselI[1,x]*2/x-1-1/2(x/2)^2)/(x/2)^4
89g[z_]:=f[2 Sqrt[z]]
90{err, approx}=MiniMaxApproximation[g[z],{z,{0.000000001,1},3,2},WorkingPrecision->60]
91poly=Numerator[approx][[1]];
92coeffs=CoefficientList[poly,z];
93TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
94poly=Denominator[approx][[1]];
95coeffs=CoefficientList[poly,z];
96TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
97```
98**/
99#[inline]
100fn i1f_small(x: f32) -> f64 {
101    let dx = x as f64;
102    let x_over_two = dx * 0.5;
103    let x_over_two_sqr = x_over_two * x_over_two;
104    let x_over_two_p4 = x_over_two_sqr * x_over_two_sqr;
105
106    let p_num = f_polyeval4(
107        x_over_two_sqr,
108        f64::from_bits(0x3fb5555555555355),
109        f64::from_bits(0x3f6ebf07f0dbc49b),
110        f64::from_bits(0x3f1fdc02bf28a8d9),
111        f64::from_bits(0x3ebb5e7574c700a6),
112    );
113    let p_den = f_polyeval3(
114        x_over_two_sqr,
115        f64::from_bits(0x3ff0000000000000),
116        f64::from_bits(0xbfa39b64b6135b5a),
117        f64::from_bits(0x3f3fa729bbe951f9),
118    );
119    let p = p_num / p_den;
120
121    let p1 = f_fmla(0.5, x_over_two_sqr, 1.);
122    let p2 = f_fmla(x_over_two_p4, p, p1);
123    p2 * x_over_two
124}
125
126/**
127Series for
128f(x) := BesselK(1, x) - Log(x)*BesselI(1, x) - 1/x
129
130Generated by Wolfram Mathematica:
131```text
132<<FunctionApproximations`
133ClearAll["Global`*"]
134f[x_]:=(BesselK[1, x]-Log[x]BesselI[1,x]-1/x)/x
135g[z_]:=f[Sqrt[z]]
136{err, approx}=MiniMaxApproximation[g[z],{z,{0.000000001,1},3,3},WorkingPrecision->60]
137poly=Numerator[approx][[1]];
138coeffs=CoefficientList[poly,z];
139TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
140poly=Denominator[approx][[1]];
141coeffs=CoefficientList[poly,z];
142TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
143```
144**/
145#[inline]
146fn k1ef_small(x: f32) -> f32 {
147    let dx = x as f64;
148    let rcp = 1. / dx;
149    let x2 = dx * dx;
150    let p_num = f_polyeval4(
151        x2,
152        f64::from_bits(0xbfd3b5b6028a83d6),
153        f64::from_bits(0xbfb3fde2c83f7cca),
154        f64::from_bits(0xbf662b2e5defbe8c),
155        f64::from_bits(0xbefa2a63cc5c4feb),
156    );
157    let p_den = f_polyeval4(
158        x2,
159        f64::from_bits(0x3ff0000000000000),
160        f64::from_bits(0xbf9833197207a7c6),
161        f64::from_bits(0x3f315663bc7330ef),
162        f64::from_bits(0xbeb9211958f6b8c3),
163    );
164    let p = p_num / p_den;
165
166    let v_exp = core_expf(x);
167    let lg = fast_logf(x);
168    let v_i = i1f_small(x);
169    let z = f_fmla(lg, v_i, rcp);
170    let z0 = f_fmla(p, dx, z);
171    (z0 * v_exp) as f32
172}
173
174/**
175Generated by Wolfram Mathematica:
176```text
177<<FunctionApproximations`
178ClearAll["Global`*"]
179f[x_]:=Sqrt[x] Exp[x] BesselK[1,x]
180g[z_]:=f[1/z]
181{err, approx}=MiniMaxApproximation[g[z],{z,{0.000000001,1},7,7},WorkingPrecision->60]
182poly=Numerator[approx][[1]];
183coeffs=CoefficientList[poly,z];
184TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
185poly=Denominator[approx][[1]];
186coeffs=CoefficientList[poly,z];
187TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
188```
189**/
190#[inline]
191fn k1ef_asympt(x: f32) -> f32 {
192    let dx = x as f64;
193    let recip = 1. / dx;
194    let r_sqrt = j1f_rsqrt(dx);
195    let p_num = f_estrin_polyeval8(
196        recip,
197        f64::from_bits(0x3ff40d931ff6270d),
198        f64::from_bits(0x402d250670ed7a6c),
199        f64::from_bits(0x404e517b9b494d38),
200        f64::from_bits(0x405cb02b7433a838),
201        f64::from_bits(0x405a03e606a1b871),
202        f64::from_bits(0x4045c98d4308dbcd),
203        f64::from_bits(0x401d115c4ce0540c),
204        f64::from_bits(0x3fd4213e72b24b3a),
205    );
206    let p_den = f_estrin_polyeval8(
207        recip,
208        f64::from_bits(0x3ff0000000000000),
209        f64::from_bits(0x402681096aa3a87d),
210        f64::from_bits(0x404623ab8d72ceea),
211        f64::from_bits(0x40530af06ea802b2),
212        f64::from_bits(0x404d526906fb9cec),
213        f64::from_bits(0x403281caca389f1b),
214        f64::from_bits(0x3ffdb93996948bb4),
215        f64::from_bits(0x3f9a009da07eb989),
216    );
217    let v = p_num / p_den;
218    let pp = v * r_sqrt;
219    pp as f32
220}
221
222#[cfg(test)]
223mod tests {
224    use super::*;
225
226    #[test]
227    fn test_k1f() {
228        assert_eq!(f_k1ef(0.00000000005423), 18439980000.0);
229        assert_eq!(f_k1ef(0.0000000043123), 231894820.0);
230        assert_eq!(f_k1ef(0.3), 4.125158);
231        assert_eq!(f_k1ef(1.89), 1.0710458);
232        assert_eq!(f_k1ef(5.89), 0.5477655);
233        assert_eq!(f_k1ef(101.89), 0.12461915);
234        assert_eq!(f_k1ef(0.), f32::INFINITY);
235        assert_eq!(f_k1ef(-0.), f32::INFINITY);
236        assert!(f_k1ef(-0.5).is_nan());
237        assert!(f_k1ef(f32::NEG_INFINITY).is_nan());
238        assert_eq!(f_k1ef(f32::INFINITY), 0.);
239    }
240}
pxfm/bessel/k1ef.rs

pxfm/bessel/
k1ef.rs
