pxfm/bessel/
k0ef.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::bessel::j0f::j1f_rsqrt;
31use crate::common::f_fmla;
32use crate::exponents::core_expf;
33use crate::logs::fast_logf;
34use crate::polyeval::{f_estrin_polyeval7, f_estrin_polyeval8};
35
36/// Modified exponentially scaled Bessel of the first kind of order 0
37///
38/// Computes K0(x)exp(x)
39///
40/// Max ULP 0.5
41pub fn f_k0ef(x: f32) -> f32 {
42    let ux = x.to_bits();
43    if ux >= 0xffu32 << 23 || ux == 0 {
44        // |x| == 0, |x| == inf, |x| == NaN, x < 0
45        if ux.wrapping_shl(1) == 0 {
46            // |x| == 0
47            return f32::INFINITY;
48        }
49        if x.is_infinite() {
50            return if x.is_sign_positive() { 0. } else { f32::NAN };
51        }
52        return x + f32::NAN; // x == NaN
53    }
54
55    let xb = x.to_bits();
56
57    if xb <= 0x3f800000u32 {
58        // x <= 1.0
59        if xb <= 0x34000000u32 {
60            // |x| <= f32::EPSILON
61            // taylor series for K0(x)exp(x) ~ (-euler_gamma + log(2) - log(x)) + (-euler_gamma + log(2) - log(x)) * x
62            let dx = x as f64;
63            let log_x = fast_logf(x);
64            const M_EULER_GAMMA_P_LOG2: f64 = f64::from_bits(0x3fbdadb014541eb2);
65            let c1 = -log_x + M_EULER_GAMMA_P_LOG2;
66            return f_fmla(c1, dx, c1) as f32;
67        }
68        return k0ef_small(x);
69    }
70
71    k0ef_asympt(x)
72}
73
74/**
75K0(x) + log(x) * I0(x) = P(x^2)
76hence
77K0(x) = P(x^2) - log(x)*I0(x)
78
79Polynomial generated by Wolfram Mathematica:
80```text
81<<FunctionApproximations`
82ClearAll["Global`*"]
83f[x_]:=BesselK[0,x]+Log[x]BesselI[0,x]
84g[z_]:=f[Sqrt[z]]
85{err, approx}=MiniMaxApproximation[g[z],{z,{0.000000001,1},6,0},WorkingPrecision->60]
86poly=Numerator[approx][[1]];
87coeffs=CoefficientList[poly,z];
88TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
89```
90**/
91#[inline]
92fn k0ef_small(x: f32) -> f32 {
93    let v_log = fast_logf(x);
94    let i0 = i0f_small(x);
95
96    let v_exp = core_expf(x);
97
98    let dx = x as f64;
99
100    let p = f_estrin_polyeval7(
101        dx * dx,
102        f64::from_bits(0x3fbdadb014541ece),
103        f64::from_bits(0x3fd1dadb01453e9c),
104        f64::from_bits(0x3f99dadb01491ac7),
105        f64::from_bits(0x3f4bb90e82a4f609),
106        f64::from_bits(0x3eef4749ebd25b10),
107        f64::from_bits(0x3e85d5b5668593af),
108        f64::from_bits(0x3e15233b0788618b),
109    );
110    let c = f_fmla(-i0, v_log, p);
111    (c * v_exp) as f32
112}
113
114/**
115Generated in Wolfram
116
117Computes sqrt(x)*exp(x)*K0(x)=Pn(1/x)/Qm(1/x)
118hence
119K0(x)exp(x) = Pn(1/x)/Qm(1/x) / sqrt(x)
120
121```text
122<<FunctionApproximations`
123ClearAll["Global`*"]
124f[x_]:=Sqrt[x] Exp[x] BesselK[0,x]
125g[z_]:=f[1/z]
126{err,approx}=MiniMaxApproximation[g[z],{z,{2^-33,1},7,7},WorkingPrecision->60]
127poly=Numerator[approx][[1]];
128coeffs=CoefficientList[poly,z];
129TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
130poly=Denominator[approx][[1]];
131coeffs=CoefficientList[poly,z];
132TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
133```
134**/
135#[inline]
136fn k0ef_asympt(x: f32) -> f32 {
137    let dx = x as f64;
138    let recip = 1. / dx;
139    let r_sqrt = j1f_rsqrt(dx);
140
141    let p_num = f_estrin_polyeval8(
142        recip,
143        f64::from_bits(0x3ff40d931ff62701),
144        f64::from_bits(0x402d8410a60e2ced),
145        f64::from_bits(0x404e9f18049bf704),
146        f64::from_bits(0x405c07682282783c),
147        f64::from_bits(0x4057379c68ce6d5e),
148        f64::from_bits(0x403ffd64a0105c4e),
149        f64::from_bits(0x400cc53ed67913b4),
150        f64::from_bits(0x3faf8cc8747a5d72),
151    );
152    let p_den = f_estrin_polyeval8(
153        recip,
154        f64::from_bits(0x3ff0000000000000),
155        f64::from_bits(0x4027ccde1d0eeb14),
156        f64::from_bits(0x40492418133aa7a7),
157        f64::from_bits(0x4057be8a004d0938),
158        f64::from_bits(0x4054cc77d1dfef26),
159        f64::from_bits(0x403fd2187097af1d),
160        f64::from_bits(0x4011c77649649e55),
161        f64::from_bits(0x3fc2080a5965ef9b),
162    );
163    let v = p_num / p_den;
164    let pp = v * 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_k0ef(2.034804e-5), 10.918679);
175        assert_eq!(f_k0ef(0.010260499), 4.743962);
176        assert_eq!(f_k0ef(0.3260499), 1.7963701);
177        assert_eq!(f_k0ef(0.72341), 1.3121376);
178        assert_eq!(f_k0ef(0.), f32::INFINITY);
179        assert_eq!(f_k0ef(-0.), f32::INFINITY);
180        assert!(f_k0ef(-0.5).is_nan());
181        assert!(f_k0ef(f32::NEG_INFINITY).is_nan());
182        assert_eq!(f_k0ef(f32::INFINITY), 0.);
183    }
184}
pxfm/bessel/k0ef.rs

pxfm/bessel/
k0ef.rs
