pxfm/bessel/
k0f.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::i0f::i0f_small;
30use crate::common::f_fmla;
31use crate::exponents::core_expf;
32use crate::logs::fast_logf;
33use crate::polyeval::{f_estrin_polyeval7, f_estrin_polyeval8};
34
35/// Modified Bessel of the second kind of order 0
36///
37/// Max ULP 0.5
38///
39/// This method have exactly one exception which is not correctly rounded with FMA.
40pub fn f_k0f(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 >= 0x42cbc4fbu32 {
56        // x > 101.88473
57        return 0.;
58    }
59
60    if xb <= 0x3f800000u32 {
61        if xb <= 0x34000000u32 {
62            // |x| < f32::EPSILON
63            // taylor series for K0(x) ~ -euler_gamma + log(2) - log(x)
64            let log_x = fast_logf(x);
65            const EULER_GAMMA_PLUS_LOG2: f64 = f64::from_bits(0x3fbdadb014541eb2);
66            return (-log_x + EULER_GAMMA_PLUS_LOG2) as f32;
67        }
68        // x <= 1.0
69        return k0f_small(x);
70    }
71
72    k0f_asympt(x)
73}
74
75/**
76K0(x) + log(x) * I0(x) = P(x^2)
77hence
78K0(x) = P(x^2) - log(x)*I0(x)
79
80Polynomial generated by Wolfram Mathematica:
81```text
82<<FunctionApproximations`
83ClearAll["Global`*"]
84f[x_]:=BesselK[0,x]+Log[x]BesselI[0,x]
85g[z_]:=f[Sqrt[z]]
86{err, approx}=MiniMaxApproximation[g[z],{z,{0.000000001,1},6,0},WorkingPrecision->60]
87poly=Numerator[approx][[1]];
88coeffs=CoefficientList[poly,z];
89TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
90```
91**/
92#[inline]
93fn k0f_small(x: f32) -> f32 {
94    let v_log = fast_logf(x);
95    let i0 = i0f_small(x);
96
97    let dx = x as f64;
98
99    let p = f_estrin_polyeval7(
100        dx * dx,
101        f64::from_bits(0x3fbdadb014541ece),
102        f64::from_bits(0x3fd1dadb01453e9c),
103        f64::from_bits(0x3f99dadb01491ac7),
104        f64::from_bits(0x3f4bb90e82a4f609),
105        f64::from_bits(0x3eef4749ebd25b10),
106        f64::from_bits(0x3e85d5b5668593af),
107        f64::from_bits(0x3e15233b0788618b),
108    );
109    let c = f_fmla(-i0, v_log, p);
110    c as f32
111}
112
113/**
114Generated in Wolfram
115
116Computes sqrt(x)*exp(x)*K0(x)=Pn(1/x)/Qm(1/x)
117hence
118K0(x) = Pn(1/x)/Qm(1/x) / (sqrt(x) * exp(x))
119
120```text
121<<FunctionApproximations`
122ClearAll["Global`*"]
123f[x_]:=Sqrt[x] Exp[x] BesselK[0,x]
124g[z_]:=f[1/z]
125{err, approx}=MiniMaxApproximation[g[z],{z,{0.0000000000001,1},7,7},WorkingPrecision->60]
126poly=Numerator[approx][[1]];
127coeffs=CoefficientList[poly,z];
128TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
129poly=Denominator[approx][[1]];
130coeffs=CoefficientList[poly,z];
131TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
132```
133**/
134#[inline]
135fn k0f_asympt(x: f32) -> f32 {
136    let dx = x as f64;
137    let recip = 1. / dx;
138    let e = core_expf(x);
139    let r_sqrt = dx.sqrt();
140
141    let p_num = f_estrin_polyeval8(
142        recip,
143        f64::from_bits(0x3ff40d931ff62701),
144        f64::from_bits(0x402d8410a62d9c17),
145        f64::from_bits(0x404e9f1804dd7e54),
146        f64::from_bits(0x405c076822dcd255),
147        f64::from_bits(0x4057379c6932949f),
148        f64::from_bits(0x403ffd64a0bd54b7),
149        f64::from_bits(0x400cc53ed733fd97),
150        f64::from_bits(0x3faf8cc8756944eb),
151    );
152    let p_den = f_estrin_polyeval8(
153        recip,
154        f64::from_bits(0x3ff0000000000000),
155        f64::from_bits(0x4027ccde1d27ffc9),
156        f64::from_bits(0x40492418136fb90f),
157        f64::from_bits(0x4057be8a00983906),
158        f64::from_bits(0x4054cc77d2379b76),
159        f64::from_bits(0x403fd218713ec08d),
160        f64::from_bits(0x4011c77649d3f65f),
161        f64::from_bits(0x3fc2080a59e87324),
162    );
163    let v = p_num / p_den;
164    let pp = v / (e * r_sqrt);
165    pp as f32
166}
167
168#[cfg(test)]
169mod tests {
170    use super::*;
171
172    #[test]
173    fn test_k0f() {
174        assert_eq!(f_k0f(2.034804e-5), 10.918458);
175        assert_eq!(f_k0f(0.010260499), 4.695535);
176        assert_eq!(f_k0f(0.3260499), 1.2965646);
177        assert_eq!(f_k0f(0.72341), 0.636511734);
178        assert_eq!(f_k0f(0.), f32::INFINITY);
179        assert_eq!(f_k0f(-0.), f32::INFINITY);
180        assert!(f_k0f(-0.5).is_nan());
181        assert!(f_k0f(f32::NEG_INFINITY).is_nan());
182        assert_eq!(f_k0f(f32::INFINITY), 0.);
183    }
184}
pxfm/bessel/k0f.rs

pxfm/bessel/
k0f.rs
