pxfm/bessel/
k0.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::double_double::DoubleDouble;
32use crate::dyadic_float::{DyadicFloat128, DyadicSign};
33use crate::exponents::rational128_exp;
34use crate::logs::{fast_log_d_to_dd, log_dd};
35use crate::polyeval::f_polyeval3;
36
37/// Modified Bessel of the second kind of order 0
38///
39/// Max ULP 0.5
40pub fn f_k0(x: f64) -> f64 {
41    let ix = x.to_bits();
42
43    if ix >= 0x7ffu64 << 52 || ix == 0 {
44        // |x| == NaN, x == inf, |x| == 0, x < 0
45        if ix.wrapping_shl(1) == 0 {
46            // |x| == 0
47            return f64::INFINITY;
48        }
49        if x.is_infinite() {
50            return if x.is_sign_positive() { 0. } else { f64::NAN };
51        }
52        return x + f64::NAN; // x == NaN
53    }
54
55    let xb = x.to_bits();
56
57    if xb >= 0x40862e42fefa39f0u64 {
58        // x >= 709.7827128933841
59        return 0.;
60    }
61
62    if xb <= 0x3ff0000000000000 {
63        // x <= 1
64        return k0_small_dd(x).to_f64();
65    }
66
67    k0_asympt(x)
68}
69
70/**
71Computes I0 on interval [0; 1]
72as rational approximation I0 = 1 + (x/2)^2 * Pn((x/2)^2)/Qm((x/2)^2))
73
74Generated by Wolfram Mathematica:
75```python
76<<FunctionApproximations`
77ClearAll["Global`*"]
78f[x_]:=(BesselI[0,x]-1)/(x/2)^2
79g[z_]:=f[2 Sqrt[z]]
80{err,approx}=MiniMaxApproximation[g[z],{z,{0.0000000000001,1},5,5},WorkingPrecision->60]
81poly=Numerator[approx][[1]];
82coeffs=CoefficientList[poly,z];
83TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
84poly=Denominator[approx][[1]];
85coeffs=CoefficientList[poly,z];
86TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
87```
88**/
89#[inline]
90fn i0_0_to_1_fast(x: f64) -> DoubleDouble {
91    let half_x = x * 0.5; // this is exact
92    let eval_x = DoubleDouble::from_exact_mult(half_x, half_x);
93
94    const P: [(u64, u64); 3] = [
95        (0xbae20452afd5045b, 0x3ff0000000000000),
96        (0xbc5b6ff3f140da20, 0x3fc93c83592c03de),
97        (0x3c25b350e9128d49, 0x3f904f33ef2de455),
98    ];
99
100    let ps_num = f_polyeval3(
101        eval_x.hi,
102        f64::from_bits(0x3f433805a2fabaaa),
103        f64::from_bits(0x3ee5897e7f554966),
104        f64::from_bits(0x3e731401f0bb5de4),
105    );
106    let mut p_num = DoubleDouble::mul_f64_add(eval_x, ps_num, DoubleDouble::from_bit_pair(P[2]));
107    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[1]));
108    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[0]));
109
110    const Q: [(u64, u64); 3] = [
111        (0x0000000000000000, 0x3ff0000000000000),
112        (0x3c323fa63bef2b4e, 0xbfab0df29b4ff089),
113        (0x3bfedbdf64ed3110, 0x3f564662064157d2),
114    ];
115
116    let ps_den = f_polyeval3(
117        eval_x.hi,
118        f64::from_bits(0xbef6bdbb484fd0a4),
119        f64::from_bits(0x3e8d6ced53309351),
120        f64::from_bits(0xbe13cff13854e945),
121    );
122    let mut p_den = DoubleDouble::mul_f64_add(eval_x, ps_den, DoubleDouble::from_bit_pair(Q[2]));
123    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[1]));
124    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[0]));
125    let p = DoubleDouble::div(p_num, p_den);
126
127    let z = DoubleDouble::quick_mult(p, eval_x);
128
129    DoubleDouble::full_add_f64(z, 1.)
130}
131
132/**
133K0(x) + log(x) * I0(x) = P(x^2)
134hence
135K0(x) = P(x^2) - log(x)*I0(x)
136
137Generated in Wolfram Mathematica:
138```text
139<<FunctionApproximations`
140ClearAll["Global`*"]
141f[x_]:=BesselK[0,x]+Log[x]BesselI[0,x]
142g[z_]:=f[Sqrt[z]]
143{err,approx}=MiniMaxApproximation[g[z],{z,{0.000000001,1},5,5},WorkingPrecision->60]
144poly=Numerator[approx][[1]];
145coeffs=CoefficientList[poly,z];
146TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
147poly=Denominator[approx][[1]];
148coeffs=CoefficientList[poly,z];
149TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
150```
151**/
152#[inline]
153pub(crate) fn k0_small_dd(x: f64) -> DoubleDouble {
154    let dx = DoubleDouble::from_exact_mult(x, x);
155    const P: [(u64, u64); 6] = [
156        (0x3c1be095d044e896, 0x3fbdadb014541eb2),
157        (0x3c7321baa1d0a2d9, 0x3fd1b9f19bc9019a),
158        (0xbc33ce33a244e5bd, 0x3f94ec39f8744183),
159        (0x3bd7008dfc623255, 0x3f3d85175b25727d),
160        (0xbb4aa2a1c4905d30, 0x3ed007a860ef3235),
161        (0xbae8daa77abd6f7f, 0x3e4839e32c19f31a),
162    ];
163
164    let ps_num = f_polyeval3(
165        dx.hi,
166        f64::from_bits(0x3f3d85175b25727d),
167        f64::from_bits(0x3ed007a860ef3235),
168        f64::from_bits(0x3e4839e32c19f31a),
169    );
170
171    let mut p_num = DoubleDouble::mul_f64_add(dx, ps_num, DoubleDouble::from_bit_pair(P[2]));
172    p_num = DoubleDouble::mul_add(dx, p_num, DoubleDouble::from_bit_pair(P[1]));
173    p_num = DoubleDouble::mul_add(dx, p_num, DoubleDouble::from_bit_pair(P[0]));
174
175    const Q: [(u64, u64); 3] = [
176        (0x0000000000000000, 0x3ff0000000000000),
177        (0xbc2a82a292acdc83, 0xbf91be3a25c968d6),
178        (0xbb9d2c37183a6496, 0x3f23bac6961619d8),
179    ];
180
181    let ps_den = f_polyeval3(
182        dx.hi,
183        f64::from_bits(0xbeac315b81faa1bf),
184        f64::from_bits(0x3e2ab2d2fbae0863),
185        f64::from_bits(0xbd9be23550f83df7),
186    );
187    let mut p_den = DoubleDouble::mul_f64_add(dx, ps_den, DoubleDouble::from_bit_pair(Q[2]));
188    p_den = DoubleDouble::mul_add(dx, p_den, DoubleDouble::from_bit_pair(Q[1]));
189    p_den = DoubleDouble::mul_add_f64(dx, p_den, f64::from_bits(0x3ff0000000000000));
190
191    let prod = DoubleDouble::div(p_num, p_den);
192
193    let vi_log = fast_log_d_to_dd(x);
194    let vi = i0_0_to_1_fast(x);
195    let r = DoubleDouble::mul_add(vi_log, -vi, prod);
196    let err = r.hi * f64::from_bits(0x3c00000000000000); // 2^-63
197    let ub = r.hi + (r.lo + err);
198    let lb = r.hi + (r.lo - err);
199    if ub == lb {
200        return r;
201    }
202
203    k0_small_hard(x, vi)
204}
205
206/**
207K0(x) + log(x) * I0(x) = P(x^2)
208hence
209K0(x) = P(x^2) - log(x)*I0(x)
210
211Generated in Wolfram Mathematica:
212```text
213<<FunctionApproximations`
214ClearAll["Global`*"]
215f[x_]:=BesselK[0,x]+Log[x]BesselI[0,x]
216g[z_]:=f[Sqrt[z]]
217{err,approx}=MiniMaxApproximation[g[z],{z,{0.000000001,1},5,5},WorkingPrecision->60]
218poly=Numerator[approx][[1]];
219coeffs=CoefficientList[poly,z];
220TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
221poly=Denominator[approx][[1]];
222coeffs=CoefficientList[poly,z];
223TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
224```
225**/
226#[cold]
227#[inline(never)]
228fn k0_small_hard(x: f64, vi: DoubleDouble) -> DoubleDouble {
229    let dx = DoubleDouble::from_exact_mult(x, x);
230    const P: [(u64, u64); 6] = [
231        (0x3c1be095d044e896, 0x3fbdadb014541eb2),
232        (0x3c7321baa1d0a2d9, 0x3fd1b9f19bc9019a),
233        (0xbc33ce33a244e5bd, 0x3f94ec39f8744183),
234        (0x3bd7008dfc623255, 0x3f3d85175b25727d),
235        (0xbb4aa2a1c4905d30, 0x3ed007a860ef3235),
236        (0xbae8daa77abd6f7f, 0x3e4839e32c19f31a),
237    ];
238
239    let mut p_num = DoubleDouble::mul_add(
240        dx,
241        DoubleDouble::from_bit_pair(P[5]),
242        DoubleDouble::from_bit_pair(P[4]),
243    );
244    p_num = DoubleDouble::mul_add(dx, p_num, DoubleDouble::from_bit_pair(P[3]));
245    p_num = DoubleDouble::mul_add(dx, p_num, DoubleDouble::from_bit_pair(P[2]));
246    p_num = DoubleDouble::mul_add(dx, p_num, DoubleDouble::from_bit_pair(P[1]));
247    p_num = DoubleDouble::mul_add(dx, p_num, DoubleDouble::from_bit_pair(P[0]));
248
249    const Q: [(u64, u64); 6] = [
250        (0x0000000000000000, 0x3ff0000000000000),
251        (0xbc2a82a292acdc83, 0xbf91be3a25c968d6),
252        (0xbb9d2c37183a6496, 0x3f23bac6961619d8),
253        (0xbb32032e14c6c2b2, 0xbeac315b81faa1bf),
254        (0x3aa1a1dc04bfba96, 0x3e2ab2d2fbae0863),
255        (0x3a3e0f678099fcff, 0xbd9be23550f83df7),
256    ];
257
258    let mut p_den = DoubleDouble::mul_add(
259        dx,
260        DoubleDouble::from_bit_pair(Q[5]),
261        DoubleDouble::from_bit_pair(Q[4]),
262    );
263    p_den = DoubleDouble::mul_add(dx, p_den, DoubleDouble::from_bit_pair(Q[3]));
264    p_den = DoubleDouble::mul_add(dx, p_den, DoubleDouble::from_bit_pair(Q[2]));
265    p_den = DoubleDouble::mul_add(dx, p_den, DoubleDouble::from_bit_pair(Q[1]));
266    p_den = DoubleDouble::mul_add_f64(dx, p_den, f64::from_bits(0x3ff0000000000000));
267
268    let prod = DoubleDouble::div(p_num, p_den);
269
270    let v_log = log_dd(x);
271
272    DoubleDouble::mul_add(v_log, -vi, prod)
273}
274
275/**
276Generated in Wolfram
277
278Computes sqrt(x)*exp(x)*K0(x)=Pn(1/x)/Qm(1/x)
279hence
280K0(x) = Pn(1/x)/Qm(1/x) / (sqrt(x) * exp(x))
281
282```text
283<<FunctionApproximations`
284ClearAll["Global`*"]
285f[x_]:=Sqrt[x] Exp[x] BesselK[0,x]
286g[z_]:=f[1/z]
287{err, approx}=MiniMaxApproximation[g[z],{z,{0.0000000000001,1},11,11},WorkingPrecision->60]
288poly=Numerator[approx][[1]];
289coeffs=CoefficientList[poly,z];
290TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
291poly=Denominator[approx][[1]];
292coeffs=CoefficientList[poly,z];
293TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
294```
295**/
296#[inline]
297fn k0_asympt(x: f64) -> f64 {
298    let recip = DoubleDouble::from_quick_recip(x);
299    let e = i0_exp(x * 0.5);
300    let r_sqrt = DoubleDouble::from_sqrt(x);
301
302    const P: [(u64, u64); 12] = [
303        (0xbc9a6a11d237114e, 0x3ff40d931ff62706),
304        (0x3cdd614ddf4929e5, 0x4040645168c3e483),
305        (0xbd1ecf9ea0af6ab2, 0x40757419a703a2ab),
306        (0xbd3da3551fb27770, 0x409d4e65365522a2),
307        (0xbd564d58855d1a46, 0x40b6dd32f5a199d9),
308        (0xbd6cf055ca963a8e, 0x40c4fd2368f19618),
309        (0x3d4b6cdfbdc058df, 0x40c68faa11ebcd59),
310        (0x3d5b4ce4665bfa46, 0x40bb6fbe08e0a8ea),
311        (0xbd4316909063be15, 0x40a1953103a5be31),
312        (0x3d12f3f8edf41af0, 0x4074d2cb001e175c),
313        (0xbcd7bba36540264f, 0x40316cffcad5f8f9),
314        (0xbc6bf28dfdd5d37d, 0x3fc2f487fe78b8d7),
315    ];
316
317    let x2 = DoubleDouble::quick_mult(recip, recip);
318    let x4 = DoubleDouble::quick_mult(x2, x2);
319    let x8 = DoubleDouble::quick_mult(x4, x4);
320
321    let e0 = DoubleDouble::mul_add(
322        recip,
323        DoubleDouble::from_bit_pair(P[1]),
324        DoubleDouble::from_bit_pair(P[0]),
325    );
326    let e1 = DoubleDouble::mul_add(
327        recip,
328        DoubleDouble::from_bit_pair(P[3]),
329        DoubleDouble::from_bit_pair(P[2]),
330    );
331    let e2 = DoubleDouble::mul_add(
332        recip,
333        DoubleDouble::from_bit_pair(P[5]),
334        DoubleDouble::from_bit_pair(P[4]),
335    );
336    let e3 = DoubleDouble::mul_add(
337        recip,
338        DoubleDouble::from_bit_pair(P[7]),
339        DoubleDouble::from_bit_pair(P[6]),
340    );
341    let e4 = DoubleDouble::mul_add(
342        recip,
343        DoubleDouble::from_bit_pair(P[9]),
344        DoubleDouble::from_bit_pair(P[8]),
345    );
346    let e5 = DoubleDouble::mul_add(
347        recip,
348        DoubleDouble::from_bit_pair(P[11]),
349        DoubleDouble::from_bit_pair(P[10]),
350    );
351
352    let f0 = DoubleDouble::mul_add(x2, e1, e0);
353    let f1 = DoubleDouble::mul_add(x2, e3, e2);
354    let f2 = DoubleDouble::mul_add(x2, e5, e4);
355
356    let g0 = DoubleDouble::mul_add(x4, f1, f0);
357
358    let p_num = DoubleDouble::mul_add(x8, f2, g0);
359
360    const Q: [(u64, u64); 12] = [
361        (0x0000000000000000, 0x3ff0000000000000),
362        (0xbcb9e8a5b17e696a, 0x403a485acd64d64a),
363        (0x3cd2e2e9c87f71f7, 0x4071518092320ecb),
364        (0xbd0d05bdb9431a2f, 0x4097e57e4c22c08e),
365        (0x3d5207068ab19ba9, 0x40b2ebadb2db62f9),
366        (0xbd64e37674083471, 0x40c1c0e4e9d6493d),
367        (0x3d3efb7a9a62b020, 0x40c3b94e8d62cdc7),
368        (0x3d47d6ce80a2114b, 0x40b93c2fd39e868e),
369        (0xbd1dfda61f525861, 0x40a1877a53a7f8d8),
370        (0x3d1236ff523dfcfa, 0x4077c3a10c2827de),
371        (0xbcc889997c9b0fe7, 0x4039a1d80b11c4a1),
372        (0x3c7ded0e8d73dddc, 0x3fdafe4913722123),
373    ];
374
375    let e0 = DoubleDouble::mul_add_f64(
376        recip,
377        DoubleDouble::from_bit_pair(Q[1]),
378        f64::from_bits(0x3ff0000000000000),
379    );
380    let e1 = DoubleDouble::mul_add(
381        recip,
382        DoubleDouble::from_bit_pair(Q[3]),
383        DoubleDouble::from_bit_pair(Q[2]),
384    );
385    let e2 = DoubleDouble::mul_add(
386        recip,
387        DoubleDouble::from_bit_pair(Q[5]),
388        DoubleDouble::from_bit_pair(Q[4]),
389    );
390    let e3 = DoubleDouble::mul_add(
391        recip,
392        DoubleDouble::from_bit_pair(Q[7]),
393        DoubleDouble::from_bit_pair(Q[6]),
394    );
395    let e4 = DoubleDouble::mul_add(
396        recip,
397        DoubleDouble::from_bit_pair(Q[9]),
398        DoubleDouble::from_bit_pair(Q[8]),
399    );
400    let e5 = DoubleDouble::mul_add(
401        recip,
402        DoubleDouble::from_bit_pair(Q[11]),
403        DoubleDouble::from_bit_pair(Q[10]),
404    );
405
406    let f0 = DoubleDouble::mul_add(x2, e1, e0);
407    let f1 = DoubleDouble::mul_add(x2, e3, e2);
408    let f2 = DoubleDouble::mul_add(x2, e5, e4);
409
410    let g0 = DoubleDouble::mul_add(x4, f1, f0);
411
412    let p_den = DoubleDouble::mul_add(x8, f2, g0);
413
414    let z = DoubleDouble::div(p_num, p_den);
415
416    let r = DoubleDouble::div(z, e * r_sqrt * e);
417
418    let err = r.hi * f64::from_bits(0x3c10000000000000); // 2^-62
419    let ub = r.hi + (r.lo + err);
420    let lb = r.hi + (r.lo - err);
421    if ub != lb {
422        return k0_asympt_hard(x);
423    }
424    r.to_f64()
425}
426
427/**
428Generated in Wolfram
429
430Computes sqrt(x)*exp(x)*K0(x)=Pn(1/x)/Qm(1/x)
431hence
432K0(x) = Pn(1/x)/Qm(1/x) / (sqrt(x) * exp(x))
433
434```text
435<<FunctionApproximations`
436ClearAll["Global`*"]
437f[x_]:=Sqrt[x] Exp[x] BesselK[0,x]
438g[z_]:=f[1/z]
439{err, approx}=MiniMaxApproximation[g[z],{z,{0.0000000000001,1},14,14},WorkingPrecision->90]
440poly=Numerator[approx][[1]];
441coeffs=CoefficientList[poly,z];
442TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
443poly=Denominator[approx][[1]];
444coeffs=CoefficientList[poly,z];
445TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
446```
447**/
448#[inline(never)]
449#[cold]
450fn k0_asympt_hard(x: f64) -> f64 {
451    static P: [DyadicFloat128; 15] = [
452        DyadicFloat128 {
453            sign: DyadicSign::Pos,
454            exponent: -127,
455            mantissa: 0xa06c98ff_b1382cb2_be520f51_a7b8f970_u128,
456        },
457        DyadicFloat128 {
458            sign: DyadicSign::Pos,
459            exponent: -122,
460            mantissa: 0xc84d8d0c_7faeef84_e56abccc_3d70f8a2_u128,
461        },
462        DyadicFloat128 {
463            sign: DyadicSign::Pos,
464            exponent: -118,
465            mantissa: 0xd1a71096_3da22280_35768c9e_0d3ddf42_u128,
466        },
467        DyadicFloat128 {
468            sign: DyadicSign::Pos,
469            exponent: -115,
470            mantissa: 0xf202e474_3698aabb_05688da0_ea1a088d_u128,
471        },
472        DyadicFloat128 {
473            sign: DyadicSign::Pos,
474            exponent: -112,
475            mantissa: 0xaaa01830_8138af4d_1137b2dd_11a216f5_u128,
476        },
477        DyadicFloat128 {
478            sign: DyadicSign::Pos,
479            exponent: -110,
480            mantissa: 0x99e0649f_320bca1a_c7adadb3_f5d8498e_u128,
481        },
482        DyadicFloat128 {
483            sign: DyadicSign::Pos,
484            exponent: -109,
485            mantissa: 0xb4d81657_de1baf00_918cbc76_c6974e96_u128,
486        },
487        DyadicFloat128 {
488            sign: DyadicSign::Pos,
489            exponent: -108,
490            mantissa: 0x8a9a28c8_a61c2c7a_12416d56_51c0b3d3_u128,
491        },
492        DyadicFloat128 {
493            sign: DyadicSign::Pos,
494            exponent: -108,
495            mantissa: 0x88a079f1_d9bd4582_6353316c_3aeb9dc9_u128,
496        },
497        DyadicFloat128 {
498            sign: DyadicSign::Pos,
499            exponent: -109,
500            mantissa: 0xa82e10eb_9dc6225a_ef6223e7_54aa254d_u128,
501        },
502        DyadicFloat128 {
503            sign: DyadicSign::Pos,
504            exponent: -111,
505            mantissa: 0xf5fc07fe_6b652e8a_0b9e74ba_d0c56118_u128,
506        },
507        DyadicFloat128 {
508            sign: DyadicSign::Pos,
509            exponent: -113,
510            mantissa: 0xc5288444_c7354b24_4a4e1663_86488928_u128,
511        },
512        DyadicFloat128 {
513            sign: DyadicSign::Pos,
514            exponent: -116,
515            mantissa: 0x96d3d226_a220ae6e_d6cca1ae_40f01e27_u128,
516        },
517        DyadicFloat128 {
518            sign: DyadicSign::Pos,
519            exponent: -121,
520            mantissa: 0xa7ab931b_499c4964_499c1091_4ab9673d_u128,
521        },
522        DyadicFloat128 {
523            sign: DyadicSign::Pos,
524            exponent: -129,
525            mantissa: 0xf8373d1a_9ff3f9c6_e5cfbe0a_85ccc131_u128,
526        },
527    ];
528
529    static Q: [DyadicFloat128; 15] = [
530        DyadicFloat128 {
531            sign: DyadicSign::Pos,
532            exponent: -127,
533            mantissa: 0x80000000_00000000_00000000_00000000_u128,
534        },
535        DyadicFloat128 {
536            sign: DyadicSign::Pos,
537            exponent: -122,
538            mantissa: 0xa05190f4_dcf0d35c_277e0f21_0635c538_u128,
539        },
540        DyadicFloat128 {
541            sign: DyadicSign::Pos,
542            exponent: -118,
543            mantissa: 0xa8837381_94c38992_86c0995d_5e5fa474_u128,
544        },
545        DyadicFloat128 {
546            sign: DyadicSign::Pos,
547            exponent: -115,
548            mantissa: 0xc3a4f279_9297e905_f59cc065_75959de8_u128,
549        },
550        DyadicFloat128 {
551            sign: DyadicSign::Pos,
552            exponent: -112,
553            mantissa: 0x8b05ade4_03432e06_881ce37d_a907216d_u128,
554        },
555        DyadicFloat128 {
556            sign: DyadicSign::Pos,
557            exponent: -111,
558            mantissa: 0xfd77f85e_35626f21_355ae728_01b78cbe_u128,
559        },
560        DyadicFloat128 {
561            sign: DyadicSign::Pos,
562            exponent: -109,
563            mantissa: 0x972ed117_254970eb_661121dc_a4462d2f_u128,
564        },
565        DyadicFloat128 {
566            sign: DyadicSign::Pos,
567            exponent: -109,
568            mantissa: 0xec9d204a_9294ab57_2ef500d5_59d701b7_u128,
569        },
570        DyadicFloat128 {
571            sign: DyadicSign::Pos,
572            exponent: -109,
573            mantissa: 0xf033522d_cae45860_53a01453_c56da895_u128,
574        },
575        DyadicFloat128 {
576            sign: DyadicSign::Pos,
577            exponent: -109,
578            mantissa: 0x9a33640c_9896ead5_1ce040d7_b36544f3_u128,
579        },
580        DyadicFloat128 {
581            sign: DyadicSign::Pos,
582            exponent: -111,
583            mantissa: 0xefe714fa_49da0166_fdf8bc68_57b77fa0_u128,
584        },
585        DyadicFloat128 {
586            sign: DyadicSign::Pos,
587            exponent: -113,
588            mantissa: 0xd323b84c_214196b0_e25b8931_930fea0d_u128,
589        },
590        DyadicFloat128 {
591            sign: DyadicSign::Pos,
592            exponent: -116,
593            mantissa: 0xbbb5240b_346642d8_010383cb_1e8a607e_u128,
594        },
595        DyadicFloat128 {
596            sign: DyadicSign::Pos,
597            exponent: -120,
598            mantissa: 0x88dcfa2a_f9f7d2ab_dd017994_8fae7e87_u128,
599        },
600        DyadicFloat128 {
601            sign: DyadicSign::Pos,
602            exponent: -127,
603            mantissa: 0xc891477c_526e0f5e_74c4ae9f_9d8732b5_u128,
604        },
605    ];
606
607    let recip = DyadicFloat128::accurate_reciprocal(x);
608    let e = rational128_exp(x);
609    let r_sqrt = bessel_rsqrt_hard(x, recip);
610
611    let mut p0 = P[14];
612    for i in (0..14).rev() {
613        p0 = recip * p0 + P[i];
614    }
615
616    let mut q = Q[14];
617    for i in (0..14).rev() {
618        q = recip * q + Q[i];
619    }
620
621    let v = p0 * q.reciprocal();
622    let r = v * e.reciprocal() * r_sqrt;
623    r.fast_as_f64()
624}
625
626#[cfg(test)]
627mod tests {
628    use super::*;
629
630    #[test]
631    fn test_k0() {
632        assert_eq!(f_k0(0.11), 2.3332678776741127);
633        assert_eq!(f_k0(0.643), 0.7241025575342853);
634        assert_eq!(f_k0(0.964), 0.4433737413379138);
635        assert_eq!(f_k0(2.964), 0.03621679838808167);
636        assert_eq!(f_k0(423.43), 7.784461905543397e-186);
637        assert_eq!(f_k0(0.), f64::INFINITY);
638        assert_eq!(f_k0(-0.), f64::INFINITY);
639        assert!(f_k0(-0.5).is_nan());
640        assert!(f_k0(f64::NEG_INFINITY).is_nan());
641        assert_eq!(f_k0(f64::INFINITY), 0.);
642    }
643}
pxfm/bessel/k0.rs

pxfm/bessel/
k0.rs
