pxfm/bessel/
k0e.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::bessel::i0::bessel_rsqrt_hard;
30use crate::bessel::i0_exp;
31use crate::bessel::k0::k0_small_dd;
32use crate::double_double::DoubleDouble;
33use crate::dyadic_float::{DyadicFloat128, DyadicSign};
34
35/// Modified exponentially scaled Bessel of the first kind of order 0
36///
37/// Computes K0(x)exp(x)
38pub fn f_k0e(x: f64) -> f64 {
39    let ix = x.to_bits();
40
41    if ix >= 0x7ffu64 << 52 || ix == 0 {
42        // |x| == NaN, x == inf, |x| == 0, x < 0
43        if ix.wrapping_shl(1) == 0 {
44            // |x| == 0
45            return f64::INFINITY;
46        }
47        if x.is_infinite() {
48            return if x.is_sign_positive() { 0. } else { f64::NAN };
49        }
50        return x + f64::NAN; // x == NaN
51    }
52
53    let xb = x.to_bits();
54
55    if xb <= 0x3ff0000000000000 {
56        // x <= 1
57        let v_k0 = k0_small_dd(x);
58        let v_exp = i0_exp(x);
59        return DoubleDouble::quick_mult(v_exp, v_k0).to_f64();
60    }
61
62    k0e_asympt(x)
63}
64
65/**
66Generated in Wolfram
67
68Computes sqrt(x)*exp(x)*K0(x)=Pn(1/x)/Qm(1/x)
69hence
70K0(x)exp(x) = Pn(1/x)/Qm(1/x) / sqrt(x)
71
72```text
73<<FunctionApproximations`
74ClearAll["Global`*"]
75f[x_]:=Sqrt[x] Exp[x] BesselK[0,x]
76g[z_]:=f[1/z]
77{err, approx}=MiniMaxApproximation[g[z],{z,{0.0000000000001,1},11,11},WorkingPrecision->60]
78poly=Numerator[approx][[1]];
79coeffs=CoefficientList[poly,z];
80TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
81poly=Denominator[approx][[1]];
82coeffs=CoefficientList[poly,z];
83TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
84```
85**/
86#[inline]
87fn k0e_asympt(x: f64) -> f64 {
88    let recip = DoubleDouble::from_quick_recip(x);
89    let r_sqrt = DoubleDouble::from_sqrt(x);
90
91    const P: [(u64, u64); 12] = [
92        (0xbc9a6a11d237114e, 0x3ff40d931ff62706),
93        (0x3cdd614ddf4929e5, 0x4040645168c3e483),
94        (0xbd1ecf9ea0af6ab2, 0x40757419a703a2ab),
95        (0xbd3da3551fb27770, 0x409d4e65365522a2),
96        (0xbd564d58855d1a46, 0x40b6dd32f5a199d9),
97        (0xbd6cf055ca963a8e, 0x40c4fd2368f19618),
98        (0x3d4b6cdfbdc058df, 0x40c68faa11ebcd59),
99        (0x3d5b4ce4665bfa46, 0x40bb6fbe08e0a8ea),
100        (0xbd4316909063be15, 0x40a1953103a5be31),
101        (0x3d12f3f8edf41af0, 0x4074d2cb001e175c),
102        (0xbcd7bba36540264f, 0x40316cffcad5f8f9),
103        (0xbc6bf28dfdd5d37d, 0x3fc2f487fe78b8d7),
104    ];
105
106    let x2 = DoubleDouble::quick_mult(recip, recip);
107    let x4 = DoubleDouble::quick_mult(x2, x2);
108    let x8 = DoubleDouble::quick_mult(x4, x4);
109
110    let e0 = DoubleDouble::mul_add(
111        recip,
112        DoubleDouble::from_bit_pair(P[1]),
113        DoubleDouble::from_bit_pair(P[0]),
114    );
115    let e1 = DoubleDouble::mul_add(
116        recip,
117        DoubleDouble::from_bit_pair(P[3]),
118        DoubleDouble::from_bit_pair(P[2]),
119    );
120    let e2 = DoubleDouble::mul_add(
121        recip,
122        DoubleDouble::from_bit_pair(P[5]),
123        DoubleDouble::from_bit_pair(P[4]),
124    );
125    let e3 = DoubleDouble::mul_add(
126        recip,
127        DoubleDouble::from_bit_pair(P[7]),
128        DoubleDouble::from_bit_pair(P[6]),
129    );
130    let e4 = DoubleDouble::mul_add(
131        recip,
132        DoubleDouble::from_bit_pair(P[9]),
133        DoubleDouble::from_bit_pair(P[8]),
134    );
135    let e5 = DoubleDouble::mul_add(
136        recip,
137        DoubleDouble::from_bit_pair(P[11]),
138        DoubleDouble::from_bit_pair(P[10]),
139    );
140
141    let f0 = DoubleDouble::mul_add(x2, e1, e0);
142    let f1 = DoubleDouble::mul_add(x2, e3, e2);
143    let f2 = DoubleDouble::mul_add(x2, e5, e4);
144
145    let g0 = DoubleDouble::mul_add(x4, f1, f0);
146
147    let p_num = DoubleDouble::mul_add(x8, f2, g0);
148
149    const Q: [(u64, u64); 12] = [
150        (0x0000000000000000, 0x3ff0000000000000),
151        (0xbcb9e8a5b17e696a, 0x403a485acd64d64a),
152        (0x3cd2e2e9c87f71f7, 0x4071518092320ecb),
153        (0xbd0d05bdb9431a2f, 0x4097e57e4c22c08e),
154        (0x3d5207068ab19ba9, 0x40b2ebadb2db62f9),
155        (0xbd64e37674083471, 0x40c1c0e4e9d6493d),
156        (0x3d3efb7a9a62b020, 0x40c3b94e8d62cdc7),
157        (0x3d47d6ce80a2114b, 0x40b93c2fd39e868e),
158        (0xbd1dfda61f525861, 0x40a1877a53a7f8d8),
159        (0x3d1236ff523dfcfa, 0x4077c3a10c2827de),
160        (0xbcc889997c9b0fe7, 0x4039a1d80b11c4a1),
161        (0x3c7ded0e8d73dddc, 0x3fdafe4913722123),
162    ];
163
164    let e0 = DoubleDouble::mul_add_f64(
165        recip,
166        DoubleDouble::from_bit_pair(Q[1]),
167        f64::from_bits(0x3ff0000000000000),
168    );
169    let e1 = DoubleDouble::mul_add(
170        recip,
171        DoubleDouble::from_bit_pair(Q[3]),
172        DoubleDouble::from_bit_pair(Q[2]),
173    );
174    let e2 = DoubleDouble::mul_add(
175        recip,
176        DoubleDouble::from_bit_pair(Q[5]),
177        DoubleDouble::from_bit_pair(Q[4]),
178    );
179    let e3 = DoubleDouble::mul_add(
180        recip,
181        DoubleDouble::from_bit_pair(Q[7]),
182        DoubleDouble::from_bit_pair(Q[6]),
183    );
184    let e4 = DoubleDouble::mul_add(
185        recip,
186        DoubleDouble::from_bit_pair(Q[9]),
187        DoubleDouble::from_bit_pair(Q[8]),
188    );
189    let e5 = DoubleDouble::mul_add(
190        recip,
191        DoubleDouble::from_bit_pair(Q[11]),
192        DoubleDouble::from_bit_pair(Q[10]),
193    );
194
195    let f0 = DoubleDouble::mul_add(x2, e1, e0);
196    let f1 = DoubleDouble::mul_add(x2, e3, e2);
197    let f2 = DoubleDouble::mul_add(x2, e5, e4);
198
199    let g0 = DoubleDouble::mul_add(x4, f1, f0);
200
201    let p_den = DoubleDouble::mul_add(x8, f2, g0);
202
203    let z = DoubleDouble::div(p_num, p_den);
204
205    let r = DoubleDouble::div(z, r_sqrt);
206
207    let err = r.hi * f64::from_bits(0x3c10000000000000); // 2^-62
208    let ub = r.hi + (r.lo + err);
209    let lb = r.hi + (r.lo - err);
210    if ub != lb {
211        return k0e_asympt_hard(x);
212    }
213    r.to_f64()
214}
215
216/**
217Generated in Wolfram
218
219Computes sqrt(x)*exp(x)*K0(x)=Pn(1/x)/Qm(1/x)
220hence
221K0(x)exp(x) = Pn(1/x)/Qm(1/x) / sqrt(x)
222
223```text
224<<FunctionApproximations`
225ClearAll["Global`*"]
226f[x_]:=Sqrt[x] Exp[x] BesselK[0,x]
227g[z_]:=f[1/z]
228{err, approx}=MiniMaxApproximation[g[z],{z,{0.0000000000001,1},14,14},WorkingPrecision->90]
229poly=Numerator[approx][[1]];
230coeffs=CoefficientList[poly,z];
231TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
232poly=Denominator[approx][[1]];
233coeffs=CoefficientList[poly,z];
234TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
235```
236**/
237#[inline(never)]
238#[cold]
239fn k0e_asympt_hard(x: f64) -> f64 {
240    static P: [DyadicFloat128; 15] = [
241        DyadicFloat128 {
242            sign: DyadicSign::Pos,
243            exponent: -127,
244            mantissa: 0xa06c98ff_b1382cb2_be520f51_a7b8f970_u128,
245        },
246        DyadicFloat128 {
247            sign: DyadicSign::Pos,
248            exponent: -122,
249            mantissa: 0xc84d8d0c_7faeef84_e56abccc_3d70f8a2_u128,
250        },
251        DyadicFloat128 {
252            sign: DyadicSign::Pos,
253            exponent: -118,
254            mantissa: 0xd1a71096_3da22280_35768c9e_0d3ddf42_u128,
255        },
256        DyadicFloat128 {
257            sign: DyadicSign::Pos,
258            exponent: -115,
259            mantissa: 0xf202e474_3698aabb_05688da0_ea1a088d_u128,
260        },
261        DyadicFloat128 {
262            sign: DyadicSign::Pos,
263            exponent: -112,
264            mantissa: 0xaaa01830_8138af4d_1137b2dd_11a216f5_u128,
265        },
266        DyadicFloat128 {
267            sign: DyadicSign::Pos,
268            exponent: -110,
269            mantissa: 0x99e0649f_320bca1a_c7adadb3_f5d8498e_u128,
270        },
271        DyadicFloat128 {
272            sign: DyadicSign::Pos,
273            exponent: -109,
274            mantissa: 0xb4d81657_de1baf00_918cbc76_c6974e96_u128,
275        },
276        DyadicFloat128 {
277            sign: DyadicSign::Pos,
278            exponent: -108,
279            mantissa: 0x8a9a28c8_a61c2c7a_12416d56_51c0b3d3_u128,
280        },
281        DyadicFloat128 {
282            sign: DyadicSign::Pos,
283            exponent: -108,
284            mantissa: 0x88a079f1_d9bd4582_6353316c_3aeb9dc9_u128,
285        },
286        DyadicFloat128 {
287            sign: DyadicSign::Pos,
288            exponent: -109,
289            mantissa: 0xa82e10eb_9dc6225a_ef6223e7_54aa254d_u128,
290        },
291        DyadicFloat128 {
292            sign: DyadicSign::Pos,
293            exponent: -111,
294            mantissa: 0xf5fc07fe_6b652e8a_0b9e74ba_d0c56118_u128,
295        },
296        DyadicFloat128 {
297            sign: DyadicSign::Pos,
298            exponent: -113,
299            mantissa: 0xc5288444_c7354b24_4a4e1663_86488928_u128,
300        },
301        DyadicFloat128 {
302            sign: DyadicSign::Pos,
303            exponent: -116,
304            mantissa: 0x96d3d226_a220ae6e_d6cca1ae_40f01e27_u128,
305        },
306        DyadicFloat128 {
307            sign: DyadicSign::Pos,
308            exponent: -121,
309            mantissa: 0xa7ab931b_499c4964_499c1091_4ab9673d_u128,
310        },
311        DyadicFloat128 {
312            sign: DyadicSign::Pos,
313            exponent: -129,
314            mantissa: 0xf8373d1a_9ff3f9c6_e5cfbe0a_85ccc131_u128,
315        },
316    ];
317
318    static Q: [DyadicFloat128; 15] = [
319        DyadicFloat128 {
320            sign: DyadicSign::Pos,
321            exponent: -127,
322            mantissa: 0x80000000_00000000_00000000_00000000_u128,
323        },
324        DyadicFloat128 {
325            sign: DyadicSign::Pos,
326            exponent: -122,
327            mantissa: 0xa05190f4_dcf0d35c_277e0f21_0635c538_u128,
328        },
329        DyadicFloat128 {
330            sign: DyadicSign::Pos,
331            exponent: -118,
332            mantissa: 0xa8837381_94c38992_86c0995d_5e5fa474_u128,
333        },
334        DyadicFloat128 {
335            sign: DyadicSign::Pos,
336            exponent: -115,
337            mantissa: 0xc3a4f279_9297e905_f59cc065_75959de8_u128,
338        },
339        DyadicFloat128 {
340            sign: DyadicSign::Pos,
341            exponent: -112,
342            mantissa: 0x8b05ade4_03432e06_881ce37d_a907216d_u128,
343        },
344        DyadicFloat128 {
345            sign: DyadicSign::Pos,
346            exponent: -111,
347            mantissa: 0xfd77f85e_35626f21_355ae728_01b78cbe_u128,
348        },
349        DyadicFloat128 {
350            sign: DyadicSign::Pos,
351            exponent: -109,
352            mantissa: 0x972ed117_254970eb_661121dc_a4462d2f_u128,
353        },
354        DyadicFloat128 {
355            sign: DyadicSign::Pos,
356            exponent: -109,
357            mantissa: 0xec9d204a_9294ab57_2ef500d5_59d701b7_u128,
358        },
359        DyadicFloat128 {
360            sign: DyadicSign::Pos,
361            exponent: -109,
362            mantissa: 0xf033522d_cae45860_53a01453_c56da895_u128,
363        },
364        DyadicFloat128 {
365            sign: DyadicSign::Pos,
366            exponent: -109,
367            mantissa: 0x9a33640c_9896ead5_1ce040d7_b36544f3_u128,
368        },
369        DyadicFloat128 {
370            sign: DyadicSign::Pos,
371            exponent: -111,
372            mantissa: 0xefe714fa_49da0166_fdf8bc68_57b77fa0_u128,
373        },
374        DyadicFloat128 {
375            sign: DyadicSign::Pos,
376            exponent: -113,
377            mantissa: 0xd323b84c_214196b0_e25b8931_930fea0d_u128,
378        },
379        DyadicFloat128 {
380            sign: DyadicSign::Pos,
381            exponent: -116,
382            mantissa: 0xbbb5240b_346642d8_010383cb_1e8a607e_u128,
383        },
384        DyadicFloat128 {
385            sign: DyadicSign::Pos,
386            exponent: -120,
387            mantissa: 0x88dcfa2a_f9f7d2ab_dd017994_8fae7e87_u128,
388        },
389        DyadicFloat128 {
390            sign: DyadicSign::Pos,
391            exponent: -127,
392            mantissa: 0xc891477c_526e0f5e_74c4ae9f_9d8732b5_u128,
393        },
394    ];
395
396    let recip = DyadicFloat128::accurate_reciprocal(x);
397    let r_sqrt = bessel_rsqrt_hard(x, recip);
398
399    let mut p0 = P[14];
400    for i in (0..14).rev() {
401        p0 = recip * p0 + P[i];
402    }
403
404    let mut q = Q[14];
405    for i in (0..14).rev() {
406        q = recip * q + Q[i];
407    }
408
409    let v = p0 * q.reciprocal();
410    let r = v * r_sqrt;
411    r.fast_as_f64()
412}
413
414#[cfg(test)]
415mod tests {
416    use super::*;
417
418    #[test]
419    fn test_k0() {
420        assert_eq!(f_k0e(0.00060324324), 7.533665613459802);
421        assert_eq!(f_k0e(0.11), 2.6045757643537244);
422        assert_eq!(f_k0e(0.643), 1.3773725807788395);
423        assert_eq!(f_k0e(0.964), 1.1625987432322884);
424        assert_eq!(f_k0e(2.964), 0.7017119941259377);
425        assert_eq!(f_k0e(423.43), 0.06088931243251448);
426        assert_eq!(f_k0e(4324235240321.43), 6.027056776336986e-7);
427        assert_eq!(k0e_asympt_hard(423.43), 0.06088931243251448);
428        assert_eq!(f_k0e(0.), f64::INFINITY);
429        assert_eq!(f_k0e(-0.), f64::INFINITY);
430        assert!(f_k0e(-0.5).is_nan());
431        assert!(f_k0e(f64::NEG_INFINITY).is_nan());
432        assert_eq!(f_k0e(f64::INFINITY), 0.);
433    }
434}
pxfm/bessel/k0e.rs

pxfm/bessel/
k0e.rs
