pxfm/bessel/
k1.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::common::f_fmla;
32use crate::double_double::DoubleDouble;
33use crate::dyadic_float::{DyadicFloat128, DyadicSign};
34use crate::exponents::rational128_exp;
35use crate::logs::{fast_log_d_to_dd, log_dd};
36use crate::polyeval::f_polyeval3;
37
38/// Modified Bessel of the second kind of order 1
39///
40/// Max ULP 0.5
41pub fn f_k1(x: f64) -> f64 {
42    let ix = x.to_bits();
43
44    if ix >= 0x7ffu64 << 52 || ix == 0 {
45        // |x| == NaN, x == inf, |x| == 0, x < 0
46        if ix.wrapping_shl(1) == 0 {
47            // |x| == 0
48            return f64::INFINITY;
49        }
50        if x.is_infinite() {
51            return if x.is_sign_positive() { 0. } else { f64::NAN };
52        }
53        return x + f64::NAN; // x == NaN
54    }
55
56    let xb = x.to_bits();
57
58    if xb >= 0x4086140538aa7d38u64 {
59        // 706.5025494880165
60        return 0.;
61    }
62
63    if xb <= 0x3ff0000000000000 {
64        // x <= 1
65        return k1_small(x).to_f64();
66    }
67
68    k1_asympt(x)
69}
70
71// Generated by Wolfram Mathematica:
72// <<FunctionApproximations`
73// ClearAll["Global`*"]
74// f[x_]:=(BesselI[1,x]*2/x-1-1/2(x/2)^2)/(x/2)^4
75// g[z_]:=f[2 Sqrt[z]]
76// {err,approx}=MiniMaxApproximation[g[z],{z,{0.000000001,1},5,5},WorkingPrecision->60]
77// poly=Numerator[approx][[1]];
78// coeffs=CoefficientList[poly,z];
79// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
80// poly=Denominator[approx][[1]];
81// coeffs=CoefficientList[poly,z];
82// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
83#[inline]
84fn i1_fast(x: f64) -> DoubleDouble {
85    let half_x = x * 0.5; // this is exact
86    let eval_x = DoubleDouble::from_exact_mult(half_x, half_x);
87
88    const P: [(u64, u64); 3] = [
89        (0x3c5555555553c008, 0x3fb5555555555555),
90        (0x3c06f1014b703de8, 0x3f6dfda17d0a2cef),
91        (0xbbc2594d655d84db, 0x3f21b2c299108f7b),
92    ];
93
94    let ps_num = f_polyeval3(
95        eval_x.hi,
96        f64::from_bits(0x3ec37625c178f5e2),
97        f64::from_bits(0x3e5843215f0d5088),
98        f64::from_bits(0x3dd97f1f45f47244),
99    );
100    let mut p_num = DoubleDouble::mul_f64_add(eval_x, ps_num, DoubleDouble::from_bit_pair(P[2]));
101    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[1]));
102    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[0]));
103
104    const Q: [(u64, u64); 3] = [
105        (0x0000000000000000, 0x3ff0000000000000),
106        (0xbc32ebd3ac0e6253, 0xbfa42c718ce308f7),
107        (0xbbe1626e81e3c1bc, 0x3f482772320eab0e),
108    ];
109
110    let ps_den = f_polyeval3(
111        eval_x.hi,
112        f64::from_bits(0xbee169811ef4f4a1),
113        f64::from_bits(0x3e6ebdab5dbe02a5),
114        f64::from_bits(0xbdeb1dbb29fec52a),
115    );
116    let mut p_den = DoubleDouble::mul_f64_add(eval_x, ps_den, DoubleDouble::from_bit_pair(Q[2]));
117    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[1]));
118    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[0]));
119    let p = DoubleDouble::div(p_num, p_den);
120
121    let eval_sqr = DoubleDouble::quick_mult(eval_x, eval_x);
122
123    let mut z = DoubleDouble::mul_f64_add_f64(eval_x, 0.5, 1.);
124    z = DoubleDouble::mul_add(p, eval_sqr, z);
125
126    let x_over_05 = DoubleDouble::from_exact_mult(x, 0.5);
127
128    DoubleDouble::quick_mult(z, x_over_05)
129}
130
131/**
132Rational approximant for
133f(x) := BesselK(1, x) - Log(x)*BesselI(1, x) - 1/x
134
135Generated by Wolfram Mathematica:
136```text
137<<FunctionApproximations`
138ClearAll["Global`*"]
139f[x_]:=(BesselK[1,x]-Log[x]BesselI[1,x]-1/x)/x
140g[z_]:=f[Sqrt[z]]
141{err,approx}=MiniMaxApproximation[g[z],{z,{0.0000000000001,1},5,4},WorkingPrecision->60]
142poly=Numerator[approx][[1]];
143coeffs=CoefficientList[poly,z];
144TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
145poly=Denominator[approx][[1]];
146coeffs=CoefficientList[poly,z];
147TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
148```
149**/
150#[inline]
151pub(crate) fn k1_small(x: f64) -> DoubleDouble {
152    let rcp = DoubleDouble::from_quick_recip(x);
153    let x2 = DoubleDouble::from_exact_mult(x, x);
154
155    const P: [(u64, u64); 6] = [
156        (0xbc7037c12b888927, 0xbfd3b5b6028a83d6),
157        (0x3c39dba459d023e5, 0xbfb4bac288cfe0cd),
158        (0x3be0575395050120, 0xbf6c4a1abe9061df),
159        (0x3b755df8e375b3d4, 0xbf0c850679678599),
160        (0xbb097e0ec926785f, 0xbe98c4a9b608ae1f),
161        (0xbaa029f31c786e81, 0xbe104efe2246ee51),
162    ];
163
164    let ps_num = f_polyeval3(
165        x2.hi,
166        f64::from_bits(0xbf0c850679678599),
167        f64::from_bits(0xbe98c4a9b608ae1f),
168        f64::from_bits(0xbe104efe2246ee51),
169    );
170
171    let mut p_num = DoubleDouble::mul_f64_add(x2, ps_num, DoubleDouble::from_bit_pair(P[2]));
172    p_num = DoubleDouble::mul_add(x2, p_num, DoubleDouble::from_bit_pair(P[1]));
173    p_num = DoubleDouble::mul_add(x2, p_num, DoubleDouble::from_bit_pair(P[0]));
174
175    const Q: [(u64, u64); 5] = [
176        (0x0000000000000000, 0x3ff0000000000000),
177        (0x3c19f62e592f3e71, 0xbf8d3bd595449ca9),
178        (0xbba8472b975a12d7, 0x3f194de71babe24a),
179        (0xbb2eec4b611c19b5, 0xbe994a5dbec84e4d),
180        (0x3a9bae2028402903, 0x3e0981ded64a954b),
181    ];
182
183    let ps_den = f_fmla(
184        x2.hi,
185        f64::from_bits(0x3e0981ded64a954b),
186        f64::from_bits(0xbe994a5dbec84e4d),
187    );
188
189    let mut p_den = DoubleDouble::mul_f64_add(x2, ps_den, DoubleDouble::from_bit_pair(Q[2]));
190    p_den = DoubleDouble::mul_add(x2, p_den, DoubleDouble::from_bit_pair(Q[1]));
191    p_den = DoubleDouble::mul_add_f64(x2, p_den, f64::from_bits(0x3ff0000000000000));
192
193    let p = DoubleDouble::div(p_num, p_den);
194
195    let lg = fast_log_d_to_dd(x);
196    let v_i = i1_fast(x);
197    let z = DoubleDouble::mul_add(v_i, lg, rcp);
198    let r = DoubleDouble::mul_f64_add(p, x, z);
199    let err = f_fmla(
200        r.hi,
201        f64::from_bits(0x3c20000000000000), // 2^-61
202        f64::from_bits(0x3a80000000000000), // 2^-87
203    );
204    let ub = r.hi + (r.lo + err);
205    let lb = r.hi + (r.lo - err);
206    if ub == lb {
207        return r;
208    }
209    k1_small_hard(x)
210}
211
212/**
213Rational approximant for
214f(x) := BesselK(1, x) - Log(x)*BesselI(1, x) - 1/x
215
216Generated by Wolfram Mathematica:
217```text
218<<FunctionApproximations`
219ClearAll["Global`*"]
220f[x_]:=(BesselK[1,x]-Log[x]BesselI[1,x]-1/x)/x
221g[z_]:=f[Sqrt[z]]
222{err,approx}=MiniMaxApproximation[g[z],{z,{0.0000000000001,1},5,4},WorkingPrecision->60]
223poly=Numerator[approx][[1]];
224coeffs=CoefficientList[poly,z];
225TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
226poly=Denominator[approx][[1]];
227coeffs=CoefficientList[poly,z];
228TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
229```
230**/
231#[cold]
232#[inline(never)]
233fn k1_small_hard(x: f64) -> DoubleDouble {
234    let rcp = DoubleDouble::from_quick_recip(x);
235    let x2 = DoubleDouble::from_exact_mult(x, x);
236
237    const P: [(u64, u64); 6] = [
238        (0xbc7037c12b888927, 0xbfd3b5b6028a83d6),
239        (0x3c39dba459d023e5, 0xbfb4bac288cfe0cd),
240        (0x3be0575395050120, 0xbf6c4a1abe9061df),
241        (0x3b755df8e375b3d4, 0xbf0c850679678599),
242        (0xbb097e0ec926785f, 0xbe98c4a9b608ae1f),
243        (0xbaa029f31c786e81, 0xbe104efe2246ee51),
244    ];
245
246    let mut p_num = DoubleDouble::mul_add(
247        x2,
248        DoubleDouble::from_bit_pair(P[5]),
249        DoubleDouble::from_bit_pair(P[4]),
250    );
251    p_num = DoubleDouble::mul_add(x2, p_num, DoubleDouble::from_bit_pair(P[3]));
252    p_num = DoubleDouble::mul_add(x2, p_num, DoubleDouble::from_bit_pair(P[2]));
253    p_num = DoubleDouble::mul_add(x2, p_num, DoubleDouble::from_bit_pair(P[1]));
254    p_num = DoubleDouble::mul_add(x2, p_num, DoubleDouble::from_bit_pair(P[0]));
255
256    const Q: [(u64, u64); 5] = [
257        (0x0000000000000000, 0x3ff0000000000000),
258        (0x3c19f62e592f3e71, 0xbf8d3bd595449ca9),
259        (0xbba8472b975a12d7, 0x3f194de71babe24a),
260        (0xbb2eec4b611c19b5, 0xbe994a5dbec84e4d),
261        (0x3a9bae2028402903, 0x3e0981ded64a954b),
262    ];
263
264    let mut p_den = DoubleDouble::mul_add(
265        x2,
266        DoubleDouble::from_bit_pair(Q[4]),
267        DoubleDouble::from_bit_pair(Q[3]),
268    );
269    p_den = DoubleDouble::mul_add(x2, p_den, DoubleDouble::from_bit_pair(Q[2]));
270    p_den = DoubleDouble::mul_add(x2, p_den, DoubleDouble::from_bit_pair(Q[1]));
271    p_den = DoubleDouble::mul_add(x2, p_den, DoubleDouble::from_bit_pair(Q[0]));
272
273    let p = DoubleDouble::div(p_num, p_den);
274
275    let lg = log_dd(x);
276    let v_i = i1_fast(x);
277    let z = DoubleDouble::mul_add(v_i, lg, rcp);
278    DoubleDouble::mul_f64_add(p, x, z)
279}
280
281/**
282Generated by Wolfram Mathematica:
283```text
284<<FunctionApproximations`
285ClearAll["Global`*"]
286f[x_]:=Sqrt[x] Exp[x] BesselK[1,x]
287g[z_]:=f[1/z]
288{err,approx}=MiniMaxApproximation[g[z],{z,{0.0000000000001,1},11,11},WorkingPrecision->60]
289poly=Numerator[approx][[1]];
290coeffs=CoefficientList[poly,z];
291TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
292poly=Denominator[approx][[1]];
293coeffs=CoefficientList[poly,z];
294TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
295```
296**/
297#[inline]
298fn k1_asympt(x: f64) -> f64 {
299    let recip = DoubleDouble::from_quick_recip(x);
300    let e = i0_exp(x * 0.5);
301    let r_sqrt = DoubleDouble::from_sqrt(x);
302
303    const P: [(u64, u64); 12] = [
304        (0xbc9a6a0690becb3b, 0x3ff40d931ff62706),
305        (0xbce573e1bbf2f0b7, 0x40402cebfab5721d),
306        (0x3d11a739b7c11e7b, 0x4074f58abc0cfbf1),
307        (0xbd2682a09ded0116, 0x409c8315f8facef2),
308        (0xbd3a19e91a120168, 0x40b65f7a4caed8b9),
309        (0x3d449c3d2b834543, 0x40c4fe41fdb4e7b8),
310        (0xbd6bdd415ac7f7e1, 0x40c7aa402d035d03),
311        (0x3d528412ff0d6b24, 0x40bf68faddd7d850),
312        (0xbd48f4bb3f61dac6, 0x40a75f5650249952),
313        (0xbd1dc534b275e309, 0x4081bddd259c0582),
314        (0xbcce5103350bd226, 0x4046c7a049014484),
315        (0x3c8935f8acd6c1d0, 0x3fef7524082b1859),
316    ];
317
318    let x2 = DoubleDouble::quick_mult(recip, recip);
319    let x4 = DoubleDouble::quick_mult(x2, x2);
320    let x8 = DoubleDouble::quick_mult(x4, x4);
321
322    let e0 = DoubleDouble::mul_add(
323        recip,
324        DoubleDouble::from_bit_pair(P[1]),
325        DoubleDouble::from_bit_pair(P[0]),
326    );
327    let e1 = DoubleDouble::mul_add(
328        recip,
329        DoubleDouble::from_bit_pair(P[3]),
330        DoubleDouble::from_bit_pair(P[2]),
331    );
332    let e2 = DoubleDouble::mul_add(
333        recip,
334        DoubleDouble::from_bit_pair(P[5]),
335        DoubleDouble::from_bit_pair(P[4]),
336    );
337    let e3 = DoubleDouble::mul_add(
338        recip,
339        DoubleDouble::from_bit_pair(P[7]),
340        DoubleDouble::from_bit_pair(P[6]),
341    );
342    let e4 = DoubleDouble::mul_add(
343        recip,
344        DoubleDouble::from_bit_pair(P[9]),
345        DoubleDouble::from_bit_pair(P[8]),
346    );
347    let e5 = DoubleDouble::mul_add(
348        recip,
349        DoubleDouble::from_bit_pair(P[11]),
350        DoubleDouble::from_bit_pair(P[10]),
351    );
352
353    let f0 = DoubleDouble::mul_add(x2, e1, e0);
354    let f1 = DoubleDouble::mul_add(x2, e3, e2);
355    let f2 = DoubleDouble::mul_add(x2, e5, e4);
356
357    let g0 = DoubleDouble::mul_add(x4, f1, f0);
358
359    let p_num = DoubleDouble::mul_add(x8, f2, g0);
360
361    const Q: [(u64, u64); 12] = [
362        (0x0000000000000000, 0x3ff0000000000000),
363        (0x3cc0d2508437b3f4, 0x40396ff483adec14),
364        (0xbd130a9c9f8a5338, 0x4070225588d8c15d),
365        (0xbceceba8fa0e65a2, 0x4095481f6684e3bb),
366        (0x3d4099f3c178fd2a, 0x40afedc8a778bf42),
367        (0xbd3a7e6a6276a3e7, 0x40bc0c060112692e),
368        (0x3d11538c155b16d8, 0x40bcb12bd1101782),
369        (0xbd5f7b04cdea2c61, 0x40b07fa363202e10),
370        (0xbce444ed035b66c6, 0x4093d6fe8f44f838),
371        (0xbcf6f88fb942b610, 0x4065c99fa44030c3),
372        (0xbcbd1d2aedee5bc9, 0x40207ffabeb00eea),
373        (0xbc39a0c8091102c9, 0x3facff3d892cd57a),
374    ];
375
376    let e0 = DoubleDouble::mul_add_f64(
377        recip,
378        DoubleDouble::from_bit_pair(Q[1]),
379        f64::from_bits(0x3ff0000000000000),
380    );
381    let e1 = DoubleDouble::mul_add(
382        recip,
383        DoubleDouble::from_bit_pair(Q[3]),
384        DoubleDouble::from_bit_pair(Q[2]),
385    );
386    let e2 = DoubleDouble::mul_add(
387        recip,
388        DoubleDouble::from_bit_pair(Q[5]),
389        DoubleDouble::from_bit_pair(Q[4]),
390    );
391    let e3 = DoubleDouble::mul_add(
392        recip,
393        DoubleDouble::from_bit_pair(Q[7]),
394        DoubleDouble::from_bit_pair(Q[6]),
395    );
396    let e4 = DoubleDouble::mul_add(
397        recip,
398        DoubleDouble::from_bit_pair(Q[9]),
399        DoubleDouble::from_bit_pair(Q[8]),
400    );
401    let e5 = DoubleDouble::mul_add(
402        recip,
403        DoubleDouble::from_bit_pair(Q[11]),
404        DoubleDouble::from_bit_pair(Q[10]),
405    );
406
407    let f0 = DoubleDouble::mul_add(x2, e1, e0);
408    let f1 = DoubleDouble::mul_add(x2, e3, e2);
409    let f2 = DoubleDouble::mul_add(x2, e5, e4);
410
411    let g0 = DoubleDouble::mul_add(x4, f1, f0);
412
413    let p_den = DoubleDouble::mul_add(x8, f2, g0);
414
415    let z = DoubleDouble::div(p_num, p_den);
416
417    let mut r_e = DoubleDouble::quick_mult(e, r_sqrt);
418    r_e = DoubleDouble::from_exact_add(r_e.hi, r_e.lo);
419    r_e = DoubleDouble::quick_mult(r_e, e);
420    r_e = DoubleDouble::from_exact_add(r_e.hi, r_e.lo);
421
422    let r = DoubleDouble::div(z, r_e);
423
424    let err = r.hi * f64::from_bits(0x3c10000000000000); // 2^-61
425    let ub = r.hi + (r.lo + err);
426    let lb = r.hi + (r.lo - err);
427    if ub != lb {
428        return k1_asympt_hard(x);
429    }
430    r.to_f64()
431}
432
433/**
434Generated by Wolfram Mathematica:
435```text
436<<FunctionApproximations`
437ClearAll["Global`*"]
438f[x_]:=Sqrt[x] Exp[x] BesselK[1,x]
439g[z_]:=f[1/z]
440{err,approx}=MiniMaxApproximation[g[z],{z,{0.0000000000001,1},14,14},WorkingPrecision->70]
441poly=Numerator[approx][[1]];
442coeffs=CoefficientList[poly,z];
443TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
444poly=Denominator[approx][[1]];
445coeffs=CoefficientList[poly,z];
446TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
447```
448**/
449#[cold]
450#[inline(never)]
451fn k1_asympt_hard(x: f64) -> f64 {
452    static P: [DyadicFloat128; 15] = [
453        DyadicFloat128 {
454            sign: DyadicSign::Pos,
455            exponent: -127,
456            mantissa: 0xa06c98ff_b1382cb2_be5210ac_f26f25d1_u128,
457        },
458        DyadicFloat128 {
459            sign: DyadicSign::Pos,
460            exponent: -122,
461            mantissa: 0xc5f546cb_659a39d0_fafbd188_36ca05b9_u128,
462        },
463        DyadicFloat128 {
464            sign: DyadicSign::Pos,
465            exponent: -118,
466            mantissa: 0xcd0b7cfa_de158d26_7084bbe9_f1bdb66d_u128,
467        },
468        DyadicFloat128 {
469            sign: DyadicSign::Pos,
470            exponent: -115,
471            mantissa: 0xeac7be2f_957d1260_8849508a_2a5a8972_u128,
472        },
473        DyadicFloat128 {
474            sign: DyadicSign::Pos,
475            exponent: -112,
476            mantissa: 0xa4d14fec_fecc6444_4c7b0287_dad71a86_u128,
477        },
478        DyadicFloat128 {
479            sign: DyadicSign::Pos,
480            exponent: -110,
481            mantissa: 0x94e3180c_01df9932_ad2acd8b_bab59c05_u128,
482        },
483        DyadicFloat128 {
484            sign: DyadicSign::Pos,
485            exponent: -109,
486            mantissa: 0xb0de10f8_74918442_94a96368_8eaa4d0d_u128,
487        },
488        DyadicFloat128 {
489            sign: DyadicSign::Pos,
490            exponent: -108,
491            mantissa: 0x8adfea76_d6dbe5d9_46bfaf83_9341f4b5_u128,
492        },
493        DyadicFloat128 {
494            sign: DyadicSign::Pos,
495            exponent: -108,
496            mantissa: 0x8f0a4337_b69b602c_cf187222_f3a3379f_u128,
497        },
498        DyadicFloat128 {
499            sign: DyadicSign::Pos,
500            exponent: -109,
501            mantissa: 0xbd4c3ebf_c2db0fad_1b425641_cc470043_u128,
502        },
503        DyadicFloat128 {
504            sign: DyadicSign::Pos,
505            exponent: -110,
506            mantissa: 0x9b14d29f_9b97e3c8_c1a7b9d0_787f0ddb_u128,
507        },
508        DyadicFloat128 {
509            sign: DyadicSign::Pos,
510            exponent: -112,
511            mantissa: 0x93e670d2_07a553ef_a90d4895_cf1b5011_u128,
512        },
513        DyadicFloat128 {
514            sign: DyadicSign::Pos,
515            exponent: -115,
516            mantissa: 0x93e0ee0a_cb4d8910_6b4d3e37_f4f9df49_u128,
517        },
518        DyadicFloat128 {
519            sign: DyadicSign::Pos,
520            exponent: -120,
521            mantissa: 0xff0ce10d_5585abd1_e8a53a12_65131ad4_u128,
522        },
523        DyadicFloat128 {
524            sign: DyadicSign::Pos,
525            exponent: -126,
526            mantissa: 0xf020536d_822cbe51_c8de095a_03367c83_u128,
527        },
528    ];
529
530    static Q: [DyadicFloat128; 15] = [
531        DyadicFloat128 {
532            sign: DyadicSign::Pos,
533            exponent: -127,
534            mantissa: 0x80000000_00000000_00000000_00000000_u128,
535        },
536        DyadicFloat128 {
537            sign: DyadicSign::Pos,
538            exponent: -122,
539            mantissa: 0x9c729dd5_4828a918_42807f58_d485a511_u128,
540        },
541        DyadicFloat128 {
542            sign: DyadicSign::Pos,
543            exponent: -118,
544            mantissa: 0x9ff6f631_0794001d_433ab0c5_d4c682a9_u128,
545        },
546        DyadicFloat128 {
547            sign: DyadicSign::Pos,
548            exponent: -115,
549            mantissa: 0xb3f81e8b_1e0e85a6_3928342e_c83088a1_u128,
550        },
551        DyadicFloat128 {
552            sign: DyadicSign::Pos,
553            exponent: -113,
554            mantissa: 0xf6b1c203_a60d4294_239ad045_2c67c224_u128,
555        },
556        DyadicFloat128 {
557            sign: DyadicSign::Pos,
558            exponent: -111,
559            mantissa: 0xd7a98b14_7a499762_abde5c38_3a5b40e4_u128,
560        },
561        DyadicFloat128 {
562            sign: DyadicSign::Pos,
563            exponent: -110,
564            mantissa: 0xf4eb8b77_a2cdc686_afd1273f_d464c8b7_u128,
565        },
566        DyadicFloat128 {
567            sign: DyadicSign::Pos,
568            exponent: -109,
569            mantissa: 0xb4c1e12a_93ee86fc_930c6f94_cfa6ac3a_u128,
570        },
571        DyadicFloat128 {
572            sign: DyadicSign::Pos,
573            exponent: -109,
574            mantissa: 0xaaeaab88_32b776b7_fdd76b0f_24349f41_u128,
575        },
576        DyadicFloat128 {
577            sign: DyadicSign::Pos,
578            exponent: -110,
579            mantissa: 0xc8ec9d61_5bf2ee9b_878b4962_4a5cee85_u128,
580        },
581        DyadicFloat128 {
582            sign: DyadicSign::Pos,
583            exponent: -111,
584            mantissa: 0x8b97bab0_3351673f_22f10d40_fd1c9ff3_u128,
585        },
586        DyadicFloat128 {
587            sign: DyadicSign::Pos,
588            exponent: -114,
589            mantissa: 0xd31cb80a_bf8cbedc_b0dcf7e7_c599f79e_u128,
590        },
591        DyadicFloat128 {
592            sign: DyadicSign::Pos,
593            exponent: -117,
594            mantissa: 0x96b354c8_69197193_ea4f608f_81943988_u128,
595        },
596        DyadicFloat128 {
597            sign: DyadicSign::Pos,
598            exponent: -122,
599            mantissa: 0x989af1bb_e48b5c44_7cd09746_f15e935a_u128,
600        },
601        DyadicFloat128 {
602            sign: DyadicSign::Pos,
603            exponent: -130,
604            mantissa: 0xb7b51326_23c29ed5_8d3dcf5a_79bd9a4f_u128,
605        },
606    ];
607
608    let recip = DyadicFloat128::accurate_reciprocal(x);
609    let e = rational128_exp(x);
610    let r_sqrt = bessel_rsqrt_hard(x, recip);
611
612    let mut p0 = P[14];
613    for i in (0..14).rev() {
614        p0 = recip * p0 + P[i];
615    }
616
617    let mut q0 = Q[14];
618    for i in (0..14).rev() {
619        q0 = recip * q0 + Q[i];
620    }
621
622    let v = p0 * q0.reciprocal();
623    let r = v * (e.reciprocal() * r_sqrt);
624    r.fast_as_f64()
625}
626
627#[cfg(test)]
628mod tests {
629    use super::*;
630    #[test]
631    fn test_k1() {
632        assert_eq!(f_k1(0.643), 1.184534109892725);
633        assert_eq!(f_k1(0.964), 0.6402280656771248);
634        assert_eq!(f_k1(2.964), 0.04192888446074039);
635        assert_eq!(f_k1(8.43), 9.824733212831289e-5);
636        assert_eq!(f_k1(16.43), 2.3142404075259965e-8);
637        assert_eq!(f_k1(423.43), 7.793648638470207e-186);
638        assert_eq!(f_k1(0.), f64::INFINITY);
639        assert_eq!(f_k1(-0.), f64::INFINITY);
640        assert!(f_k1(-0.5).is_nan());
641        assert!(f_k1(f64::NEG_INFINITY).is_nan());
642        assert_eq!(f_k1(f64::INFINITY), 0.);
643    }
644}
pxfm/bessel/k1.rs

pxfm/bessel/
k1.rs
