pxfm/bessel/
k1e.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::k1::k1_small;
32use crate::double_double::DoubleDouble;
33use crate::dyadic_float::{DyadicFloat128, DyadicSign};
34
35/// Modified exponentially scaled Bessel of the second kind of order 1
36///
37/// Computes K1(x)exp(x)
38pub fn f_k1e(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_exp = i0_exp(x);
58        let v_k = k1_small(x);
59        return DoubleDouble::quick_mult(v_exp, v_k).to_f64();
60    }
61
62    k1e_asympt(x)
63}
64
65/**
66Generated by Wolfram Mathematica:
67```text
68<<FunctionApproximations`
69ClearAll["Global`*"]
70f[x_]:=Sqrt[x] Exp[x] BesselK[1,x]
71g[z_]:=f[1/z]
72{err,approx}=MiniMaxApproximation[g[z],{z,{0.0000000000001,1},11,11},WorkingPrecision->60]
73poly=Numerator[approx][[1]];
74coeffs=CoefficientList[poly,z];
75TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
76poly=Denominator[approx][[1]];
77coeffs=CoefficientList[poly,z];
78TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
79```
80**/
81#[inline]
82fn k1e_asympt(x: f64) -> f64 {
83    let recip = DoubleDouble::from_quick_recip(x);
84    let r_sqrt = DoubleDouble::from_sqrt(x);
85
86    const P: [(u64, u64); 12] = [
87        (0xbc9a6a0690becb3b, 0x3ff40d931ff62706),
88        (0xbce573e1bbf2f0b7, 0x40402cebfab5721d),
89        (0x3d11a739b7c11e7b, 0x4074f58abc0cfbf1),
90        (0xbd2682a09ded0116, 0x409c8315f8facef2),
91        (0xbd3a19e91a120168, 0x40b65f7a4caed8b9),
92        (0x3d449c3d2b834543, 0x40c4fe41fdb4e7b8),
93        (0xbd6bdd415ac7f7e1, 0x40c7aa402d035d03),
94        (0x3d528412ff0d6b24, 0x40bf68faddd7d850),
95        (0xbd48f4bb3f61dac6, 0x40a75f5650249952),
96        (0xbd1dc534b275e309, 0x4081bddd259c0582),
97        (0xbcce5103350bd226, 0x4046c7a049014484),
98        (0x3c8935f8acd6c1d0, 0x3fef7524082b1859),
99    ];
100
101    let x2 = DoubleDouble::quick_mult(recip, recip);
102    let x4 = DoubleDouble::quick_mult(x2, x2);
103    let x8 = DoubleDouble::quick_mult(x4, x4);
104
105    let e0 = DoubleDouble::mul_add(
106        recip,
107        DoubleDouble::from_bit_pair(P[1]),
108        DoubleDouble::from_bit_pair(P[0]),
109    );
110    let e1 = DoubleDouble::mul_add(
111        recip,
112        DoubleDouble::from_bit_pair(P[3]),
113        DoubleDouble::from_bit_pair(P[2]),
114    );
115    let e2 = DoubleDouble::mul_add(
116        recip,
117        DoubleDouble::from_bit_pair(P[5]),
118        DoubleDouble::from_bit_pair(P[4]),
119    );
120    let e3 = DoubleDouble::mul_add(
121        recip,
122        DoubleDouble::from_bit_pair(P[7]),
123        DoubleDouble::from_bit_pair(P[6]),
124    );
125    let e4 = DoubleDouble::mul_add(
126        recip,
127        DoubleDouble::from_bit_pair(P[9]),
128        DoubleDouble::from_bit_pair(P[8]),
129    );
130    let e5 = DoubleDouble::mul_add(
131        recip,
132        DoubleDouble::from_bit_pair(P[11]),
133        DoubleDouble::from_bit_pair(P[10]),
134    );
135
136    let f0 = DoubleDouble::mul_add(x2, e1, e0);
137    let f1 = DoubleDouble::mul_add(x2, e3, e2);
138    let f2 = DoubleDouble::mul_add(x2, e5, e4);
139
140    let g0 = DoubleDouble::mul_add(x4, f1, f0);
141
142    let p_num = DoubleDouble::mul_add(x8, f2, g0);
143
144    const Q: [(u64, u64); 12] = [
145        (0x0000000000000000, 0x3ff0000000000000),
146        (0x3cc0d2508437b3f4, 0x40396ff483adec14),
147        (0xbd130a9c9f8a5338, 0x4070225588d8c15d),
148        (0xbceceba8fa0e65a2, 0x4095481f6684e3bb),
149        (0x3d4099f3c178fd2a, 0x40afedc8a778bf42),
150        (0xbd3a7e6a6276a3e7, 0x40bc0c060112692e),
151        (0x3d11538c155b16d8, 0x40bcb12bd1101782),
152        (0xbd5f7b04cdea2c61, 0x40b07fa363202e10),
153        (0xbce444ed035b66c6, 0x4093d6fe8f44f838),
154        (0xbcf6f88fb942b610, 0x4065c99fa44030c3),
155        (0xbcbd1d2aedee5bc9, 0x40207ffabeb00eea),
156        (0xbc39a0c8091102c9, 0x3facff3d892cd57a),
157    ];
158
159    let e0 = DoubleDouble::mul_add_f64(
160        recip,
161        DoubleDouble::from_bit_pair(Q[1]),
162        f64::from_bits(0x3ff0000000000000),
163    );
164    let e1 = DoubleDouble::mul_add(
165        recip,
166        DoubleDouble::from_bit_pair(Q[3]),
167        DoubleDouble::from_bit_pair(Q[2]),
168    );
169    let e2 = DoubleDouble::mul_add(
170        recip,
171        DoubleDouble::from_bit_pair(Q[5]),
172        DoubleDouble::from_bit_pair(Q[4]),
173    );
174    let e3 = DoubleDouble::mul_add(
175        recip,
176        DoubleDouble::from_bit_pair(Q[7]),
177        DoubleDouble::from_bit_pair(Q[6]),
178    );
179    let e4 = DoubleDouble::mul_add(
180        recip,
181        DoubleDouble::from_bit_pair(Q[9]),
182        DoubleDouble::from_bit_pair(Q[8]),
183    );
184    let e5 = DoubleDouble::mul_add(
185        recip,
186        DoubleDouble::from_bit_pair(Q[11]),
187        DoubleDouble::from_bit_pair(Q[10]),
188    );
189
190    let f0 = DoubleDouble::mul_add(x2, e1, e0);
191    let f1 = DoubleDouble::mul_add(x2, e3, e2);
192    let f2 = DoubleDouble::mul_add(x2, e5, e4);
193
194    let g0 = DoubleDouble::mul_add(x4, f1, f0);
195
196    let p_den = DoubleDouble::mul_add(x8, f2, g0);
197
198    let z = DoubleDouble::div(p_num, p_den);
199
200    let r = DoubleDouble::div(z, r_sqrt);
201
202    let err = r.hi * f64::from_bits(0x3c10000000000000); // 2^-61
203    let ub = r.hi + (r.lo + err);
204    let lb = r.hi + (r.lo - err);
205    if ub != lb {
206        return k1e_asympt_hard(x);
207    }
208    r.to_f64()
209}
210
211/**
212Generated by Wolfram Mathematica:
213```text
214<<FunctionApproximations`
215ClearAll["Global`*"]
216f[x_]:=Sqrt[x] Exp[x] BesselK[1,x]
217g[z_]:=f[1/z]
218{err,approx}=MiniMaxApproximation[g[z],{z,{0.0000000000001,1},14,14},WorkingPrecision->70]
219poly=Numerator[approx][[1]];
220coeffs=CoefficientList[poly,z];
221TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
222poly=Denominator[approx][[1]];
223coeffs=CoefficientList[poly,z];
224TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
225```
226**/
227#[cold]
228#[inline(never)]
229fn k1e_asympt_hard(x: f64) -> f64 {
230    static P: [DyadicFloat128; 15] = [
231        DyadicFloat128 {
232            sign: DyadicSign::Pos,
233            exponent: -127,
234            mantissa: 0xa06c98ff_b1382cb2_be5210ac_f26f25d1_u128,
235        },
236        DyadicFloat128 {
237            sign: DyadicSign::Pos,
238            exponent: -122,
239            mantissa: 0xc5f546cb_659a39d0_fafbd188_36ca05b9_u128,
240        },
241        DyadicFloat128 {
242            sign: DyadicSign::Pos,
243            exponent: -118,
244            mantissa: 0xcd0b7cfa_de158d26_7084bbe9_f1bdb66d_u128,
245        },
246        DyadicFloat128 {
247            sign: DyadicSign::Pos,
248            exponent: -115,
249            mantissa: 0xeac7be2f_957d1260_8849508a_2a5a8972_u128,
250        },
251        DyadicFloat128 {
252            sign: DyadicSign::Pos,
253            exponent: -112,
254            mantissa: 0xa4d14fec_fecc6444_4c7b0287_dad71a86_u128,
255        },
256        DyadicFloat128 {
257            sign: DyadicSign::Pos,
258            exponent: -110,
259            mantissa: 0x94e3180c_01df9932_ad2acd8b_bab59c05_u128,
260        },
261        DyadicFloat128 {
262            sign: DyadicSign::Pos,
263            exponent: -109,
264            mantissa: 0xb0de10f8_74918442_94a96368_8eaa4d0d_u128,
265        },
266        DyadicFloat128 {
267            sign: DyadicSign::Pos,
268            exponent: -108,
269            mantissa: 0x8adfea76_d6dbe5d9_46bfaf83_9341f4b5_u128,
270        },
271        DyadicFloat128 {
272            sign: DyadicSign::Pos,
273            exponent: -108,
274            mantissa: 0x8f0a4337_b69b602c_cf187222_f3a3379f_u128,
275        },
276        DyadicFloat128 {
277            sign: DyadicSign::Pos,
278            exponent: -109,
279            mantissa: 0xbd4c3ebf_c2db0fad_1b425641_cc470043_u128,
280        },
281        DyadicFloat128 {
282            sign: DyadicSign::Pos,
283            exponent: -110,
284            mantissa: 0x9b14d29f_9b97e3c8_c1a7b9d0_787f0ddb_u128,
285        },
286        DyadicFloat128 {
287            sign: DyadicSign::Pos,
288            exponent: -112,
289            mantissa: 0x93e670d2_07a553ef_a90d4895_cf1b5011_u128,
290        },
291        DyadicFloat128 {
292            sign: DyadicSign::Pos,
293            exponent: -115,
294            mantissa: 0x93e0ee0a_cb4d8910_6b4d3e37_f4f9df49_u128,
295        },
296        DyadicFloat128 {
297            sign: DyadicSign::Pos,
298            exponent: -120,
299            mantissa: 0xff0ce10d_5585abd1_e8a53a12_65131ad4_u128,
300        },
301        DyadicFloat128 {
302            sign: DyadicSign::Pos,
303            exponent: -126,
304            mantissa: 0xf020536d_822cbe51_c8de095a_03367c83_u128,
305        },
306    ];
307
308    static Q: [DyadicFloat128; 15] = [
309        DyadicFloat128 {
310            sign: DyadicSign::Pos,
311            exponent: -127,
312            mantissa: 0x80000000_00000000_00000000_00000000_u128,
313        },
314        DyadicFloat128 {
315            sign: DyadicSign::Pos,
316            exponent: -122,
317            mantissa: 0x9c729dd5_4828a918_42807f58_d485a511_u128,
318        },
319        DyadicFloat128 {
320            sign: DyadicSign::Pos,
321            exponent: -118,
322            mantissa: 0x9ff6f631_0794001d_433ab0c5_d4c682a9_u128,
323        },
324        DyadicFloat128 {
325            sign: DyadicSign::Pos,
326            exponent: -115,
327            mantissa: 0xb3f81e8b_1e0e85a6_3928342e_c83088a1_u128,
328        },
329        DyadicFloat128 {
330            sign: DyadicSign::Pos,
331            exponent: -113,
332            mantissa: 0xf6b1c203_a60d4294_239ad045_2c67c224_u128,
333        },
334        DyadicFloat128 {
335            sign: DyadicSign::Pos,
336            exponent: -111,
337            mantissa: 0xd7a98b14_7a499762_abde5c38_3a5b40e4_u128,
338        },
339        DyadicFloat128 {
340            sign: DyadicSign::Pos,
341            exponent: -110,
342            mantissa: 0xf4eb8b77_a2cdc686_afd1273f_d464c8b7_u128,
343        },
344        DyadicFloat128 {
345            sign: DyadicSign::Pos,
346            exponent: -109,
347            mantissa: 0xb4c1e12a_93ee86fc_930c6f94_cfa6ac3a_u128,
348        },
349        DyadicFloat128 {
350            sign: DyadicSign::Pos,
351            exponent: -109,
352            mantissa: 0xaaeaab88_32b776b7_fdd76b0f_24349f41_u128,
353        },
354        DyadicFloat128 {
355            sign: DyadicSign::Pos,
356            exponent: -110,
357            mantissa: 0xc8ec9d61_5bf2ee9b_878b4962_4a5cee85_u128,
358        },
359        DyadicFloat128 {
360            sign: DyadicSign::Pos,
361            exponent: -111,
362            mantissa: 0x8b97bab0_3351673f_22f10d40_fd1c9ff3_u128,
363        },
364        DyadicFloat128 {
365            sign: DyadicSign::Pos,
366            exponent: -114,
367            mantissa: 0xd31cb80a_bf8cbedc_b0dcf7e7_c599f79e_u128,
368        },
369        DyadicFloat128 {
370            sign: DyadicSign::Pos,
371            exponent: -117,
372            mantissa: 0x96b354c8_69197193_ea4f608f_81943988_u128,
373        },
374        DyadicFloat128 {
375            sign: DyadicSign::Pos,
376            exponent: -122,
377            mantissa: 0x989af1bb_e48b5c44_7cd09746_f15e935a_u128,
378        },
379        DyadicFloat128 {
380            sign: DyadicSign::Pos,
381            exponent: -130,
382            mantissa: 0xb7b51326_23c29ed5_8d3dcf5a_79bd9a4f_u128,
383        },
384    ];
385
386    let recip = DyadicFloat128::accurate_reciprocal(x);
387    let r_sqrt = bessel_rsqrt_hard(x, recip);
388
389    let mut p0 = P[14];
390    for i in (0..14).rev() {
391        p0 = recip * p0 + P[i];
392    }
393
394    let mut q0 = Q[14];
395    for i in (0..14).rev() {
396        q0 = recip * q0 + Q[i];
397    }
398
399    let v = p0 * q0.reciprocal();
400    let r = v * r_sqrt;
401    r.fast_as_f64()
402}
403
404#[cfg(test)]
405mod tests {
406    use super::*;
407    #[test]
408    fn test_k1() {
409        assert_eq!(f_k1e(0.643), 2.253195748291852);
410        assert_eq!(f_k1e(0.964), 1.6787831013451477);
411        assert_eq!(f_k1e(2.964), 0.8123854795542738);
412        assert_eq!(f_k1e(8.43), 0.4502184086111872);
413        assert_eq!(f_k1e(16.43), 0.3161307996938612);
414        assert_eq!(f_k1e(423.43), 0.06096117017402597);
415        assert_eq!(f_k1e(9044.431), 0.01317914752085687);
416        assert_eq!(k1e_asympt_hard(16.43), 0.3161307996938612);
417        assert_eq!(k1e_asympt_hard(423.43), 0.06096117017402597);
418        assert_eq!(k1e_asympt_hard(9044.431), 0.01317914752085687);
419        assert_eq!(f_k1e(0.), f64::INFINITY);
420        assert_eq!(f_k1e(-0.), f64::INFINITY);
421        assert!(f_k1e(-0.5).is_nan());
422        assert!(f_k1e(f64::NEG_INFINITY).is_nan());
423        assert_eq!(f_k1e(f64::INFINITY), 0.);
424    }
425}
pxfm/bessel/k1e.rs

pxfm/bessel/
k1e.rs
