pxfm/bessel/
k2f.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::common::f_fmla;
30use crate::exponents::core_expf;
31use crate::logs::fast_logf;
32use crate::polyeval::{f_estrin_polyeval5, f_estrin_polyeval8, f_polyeval4, f_polyeval11};
33
34/// Modified Bessel of the second kind of order 2
35///
36/// ulp 0.5
37pub fn f_k2f(x: f32) -> f32 {
38    let ux = x.to_bits();
39    if ux >= 0xffu32 << 23 || ux == 0 {
40        // |x| == 0, |x| == inf, |x| == NaN, x < 0
41        if ux.wrapping_shl(1) == 0 {
42            // |x| == 0
43            return f32::INFINITY;
44        }
45        if x.is_infinite() {
46            return if x.is_sign_positive() { 0. } else { f32::NAN };
47        }
48        return x + f32::NAN; // x == NaN
49    }
50
51    let xb = x.to_bits();
52
53    if xb >= 0x42cbceefu32 {
54        // |x| >= 101.90417
55        return 0.;
56    }
57
58    if xb <= 0x3f800000u32 {
59        // x <= 1.0
60        if xb <= 0x3e9eb852u32 {
61            // x <= 0.31
62            if xb <= 0x34000000u32 {
63                // x <= f32::EPSILON
64                let dx = x as f64;
65                let r = 2. / (dx * dx);
66                return r as f32;
67            }
68            return k2f_tiny(x);
69        }
70        return k2f_small(x);
71    }
72
73    k2f_asympt(x)
74}
75
76#[inline]
77fn k2f_tiny(x: f32) -> f32 {
78    // Power series at zero for K2
79    // 2.0000000000000000/x^2-0.50000000000000000-0.12500000000000000 (-0.86593151565841245+1.0000000000000000 Log[x]) x^2-0.010416666666666667 (-1.5325981823250791+1.0000000000000000 Log[x]) x^4-0.00032552083333333333 (-1.9075981823250791+1.0000000000000000 Log[x]) x^6-0.0000054253472222222222 (-2.1742648489917458+1.0000000000000000 Log[x]) x^8+O[x]^9
80    //-0.50000000000000000+2.0000000000000000/x^2 + a3 * x^8 + x^6 * a2 + x^4 * a1 + x^2 * a0
81    let dx = x as f64;
82    let log_x = fast_logf(x);
83    let a0 = f_fmla(-4.0000000000000000, log_x, 3.4637260626336498) * 0.031250000000000000;
84    let a1 = f_fmla(-12.000000000000000, log_x, 18.391178187900949) * 0.00086805555555555556;
85    let a2 = f_fmla(-24.000000000000000, log_x, 45.782356375801899) * 0.000013563368055555556;
86    let a3 = (log_x - 2.1742648489917458) * (-0.0000054253472222222222);
87    let dx_sqr = dx * dx;
88    let two_over_dx = 2. / dx_sqr;
89    let p = f_polyeval4(dx_sqr, a0, a1, a2, a3);
90    let r = f_fmla(p, dx_sqr, two_over_dx) - 0.5;
91    r as f32
92}
93
94/**
95Computes
96I2(x) = x^2 * R(x^2)
97
98Generated by Wolfram Mathematica:
99
100```text
101<<FunctionApproximations`
102ClearAll["Global`*"]
103f[x_]:=BesselI[2,x]/x^2
104g[z_]:=f[Sqrt[z]]
105{err,approx}=MiniMaxApproximation[g[z],{z,{0.3,1},4,4},WorkingPrecision->75]
106poly=Numerator[approx][[1]];
107coeffs=CoefficientList[poly,z];
108TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
109poly=Denominator[approx][[1]];
110coeffs=CoefficientList[poly,z];
111TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
112```
113**/
114#[inline]
115fn i2f_small(x: f32) -> f64 {
116    let dx = x as f64;
117    let x_sqr = dx * dx;
118
119    let p_num = f_estrin_polyeval5(
120        x_sqr,
121        f64::from_bits(0x3fc0000000000000),
122        f64::from_bits(0x3f81520c0669099e),
123        f64::from_bits(0x3f27310bf5c5e9b0),
124        f64::from_bits(0x3eb8e2947e0a6098),
125        f64::from_bits(0x3e336dfad46e2f35),
126    );
127    let p_den = f_estrin_polyeval5(
128        x_sqr,
129        f64::from_bits(0x3ff0000000000000),
130        f64::from_bits(0xbf900d253bb12edc),
131        f64::from_bits(0x3f1ed3d9ab228297),
132        f64::from_bits(0xbea14e6660c00303),
133        f64::from_bits(0x3e13eb951a6cf38f),
134    );
135    let p = p_num / p_den;
136    p * x_sqr
137}
138
139/**
140Series for
141R(x^2) := (BesselK(2, x) - Log(x)*BesselI(2, x) - 2/x^2)/(1+x^2)
142
143Generated by Wolfram Mathematica:
144```text
145<<FunctionApproximations`
146ClearAll["Global`*"]
147f[x_]:=(BesselK[2,x]-Log[x]BesselI[2,x]-2/(x^2))/(1+x^2)
148g[z_]:=f[Sqrt[z]]
149{err,approx}=MiniMaxApproximation[g[z],{z,{0.3,1.0},10,10},WorkingPrecision->60]
150poly=Numerator[approx][[1]];
151coeffs=CoefficientList[poly,z];
152TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
153poly=Denominator[approx][[1]];
154coeffs=CoefficientList[poly,z];
155TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
156```
157**/
158#[inline]
159fn k2f_small(x: f32) -> f32 {
160    let dx = x as f64;
161    let dx_sqr = dx * dx;
162    let p_num = f_polyeval11(
163        dx_sqr,
164        f64::from_bits(0xbfdff794c9ee3b5c),
165        f64::from_bits(0xc047d3276f18e5d2),
166        f64::from_bits(0xc09200ed3702875a),
167        f64::from_bits(0xc0c39f395c47be27),
168        f64::from_bits(0xc0e0ec95bd1a3192),
169        f64::from_bits(0xc0e5973cb871c8d0),
170        f64::from_bits(0xc0cdaf528de00d53),
171        f64::from_bits(0xc0afe6d3009de17c),
172        f64::from_bits(0xc098417b22844112),
173        f64::from_bits(0x4025c45260bb1b6a),
174        f64::from_bits(0x402f2bf6b95ffe0c),
175    );
176    let p_den = f_polyeval11(
177        dx_sqr,
178        f64::from_bits(0x3ff0000000000000),
179        f64::from_bits(0x405879a43b253224),
180        f64::from_bits(0x40a3a501408a0198),
181        f64::from_bits(0x40d8172abc4a8ccc),
182        f64::from_bits(0x40f9fcb05e98bdbd),
183        f64::from_bits(0x4109c45b54be586b),
184        f64::from_bits(0x4106ad7023dd0b90),
185        f64::from_bits(0x40ed7e988d2ba5a9),
186        f64::from_bits(0x40966305e1c1123a),
187        f64::from_bits(0xc090832b6a87317c),
188        f64::from_bits(0x403b48eb703f4644),
189    );
190    let p = p_num / p_den;
191
192    let two_over_dx_sqr = 2. / dx_sqr;
193
194    let lg = fast_logf(x);
195    let v_i = i2f_small(x);
196    let z = f_fmla(lg, v_i, two_over_dx_sqr);
197    let z0 = f_fmla(p, f_fmla(dx, dx, 1.), z);
198    z0 as f32
199}
200
201/**
202Generated by Wolfram Mathematica:
203```text
204<<FunctionApproximations`
205ClearAll["Global`*"]
206f[x_]:=Sqrt[x] Exp[x] BesselK[1,x]
207g[z_]:=f[1/z]
208{err, approx}=MiniMaxApproximation[g[z],{z,{0.000000001,1},7,7},WorkingPrecision->60]
209poly=Numerator[approx][[1]];
210coeffs=CoefficientList[poly,z];
211TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
212poly=Denominator[approx][[1]];
213coeffs=CoefficientList[poly,z];
214TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
215```
216**/
217#[inline]
218fn k2f_asympt(x: f32) -> f32 {
219    let dx = x as f64;
220    let recip = 1. / dx;
221    let e = core_expf(x);
222    let r_sqrt = dx.sqrt();
223    let p_num = f_estrin_polyeval8(
224        recip,
225        f64::from_bits(0x3ff40d931ff626f2),
226        f64::from_bits(0x402d954dceb445df),
227        f64::from_bits(0x405084ea6680d028),
228        f64::from_bits(0x406242344a8ea488),
229        f64::from_bits(0x406594aa56f50fea),
230        f64::from_bits(0x405aa04eb4f0af1c),
231        f64::from_bits(0x403dd3e8e63849ef),
232        f64::from_bits(0x4004e85453648d43),
233    );
234    let p_den = f_estrin_polyeval8(
235        recip,
236        f64::from_bits(0x3ff0000000000000),
237        f64::from_bits(0x4023da9f4e05358e),
238        f64::from_bits(0x4040a4e4ceb523c9),
239        f64::from_bits(0x404725c423c9f990),
240        f64::from_bits(0x403a60c00deededc),
241        f64::from_bits(0x40149975b84c3946),
242        f64::from_bits(0x3fc69439846db871),
243        f64::from_bits(0xbf6400819bac6f45),
244    );
245    let v = p_num / p_den;
246    let pp = v / (e * r_sqrt);
247    pp as f32
248}
249
250#[cfg(test)]
251mod tests {
252    use super::*;
253
254    #[test]
255    fn test_k2f() {
256        assert!(f_k2f(-1.).is_nan());
257        assert!(f_k2f(f32::NAN).is_nan());
258        assert_eq!(f_k2f(0.), f32::INFINITY);
259        assert_eq!(f_k2f(0.65), 4.3059196);
260        assert_eq!(f_k2f(1.65), 0.44830766);
261    }
262}
pxfm/bessel/k2f.rs

pxfm/bessel/
k2f.rs
