Query cyl_bessel_k for cpp_double_double

[ckormanyos]

Hi John (@jzmaddock) I'm down to only 2 remaining functions having troubles with the candidate cpp_double_double.

On cyl_bessel_k(-1000, T(700.0)) I have issues that I can't resolve alone. The result is not bad, having $O(10^{-31}$, but the intermediate steps do some underflowing.

This calculation breaks the min_exponent10/max_exponent10 near-symmetry, so the proposed (big) change in precision.hpp has motivated by cpp_double_double having min_exponent10/max_exponent10 of $-291/+308$.

The real problem, however came when starting the recursion for $K_{1000}(700)$. The value of $K_0(700)$ came back subnormal and the value of $K_1(0)$ was zero. This messed up the recursion start and lead to the wrong answer, regardless of how exp(-x) was calculated (halved or not).

I know there would theoreticlaly be a way to do the recursion with another scaling, but that would be a very large change in Math requiring start points for scaled/unscaled $K_0(x)$ and $K_1(x)$.

Do you see a way that the unusual cpp_double_double can actually with modest impact on Math calculate $K_{1000}(700)$?

Cc: @mborland

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.83%. Comparing base (52b029f) to head (2a7bef3).

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #1228      +/-   ##
===========================================
- Coverage    93.83%   93.83%   -0.01%     
===========================================
  Files          657      657              
  Lines        55243    55231      -12     
===========================================
- Hits         51838    51826      -12     
  Misses        3405     3405              

Files with missing lines	Coverage Δ
include/boost/math/tools/precision.hpp	95.74% <ø> (-0.87%)	⬇️
test/test_bessel_k.hpp	100.00% <ø> (ø)

---

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 52b029f...2a7bef3. Read the comment docs.

[jzmaddock]

The real problem, however came when starting the recursion for K 1000 ( 700 ) . The value of K 0 ( 700 ) came back subnormal and the value of K 1 ( 0 ) was zero. This messed up the recursion start and lead to the wrong answer, regardless of how exp(-x) was calculated (halved or not).

I know there would theoreticlaly be a way to do the recursion with another scaling, but that would be a very large change in Math requiring start points for scaled/unscaled K 0 ( x ) and K 1 ( x ) .

Do you see a way that the unusual cpp_double_double can actually with modest impact on Math calculate K 1000 ( 700 ) ?

This is a Math lib bug, and I've filed an issue with reduced test case here: #1229.

It's not trivial to fix though, and I'm not sure how important the use case is, so my gut feeling is to push this to the long grass for the minute.

[ckormanyos]

It's not trivial to fix though, and I'm not sure how important the use case is, so my gut feeling is to push this to the long grass for the minute.

Thank you John. That seems like a wise idea since scaling the recursion or other heavy-handed activities seem like a bit much.

[ckormanyos]

Hi John,

I thought it might be easiest to keep the test point (the arrays and all the indexing stay the same), but change the order and argument. The test case appears in three individual areas.

I'm liking

$$ K_{-129}\left(200\right)~{\approx}~3.61744436315860678558682169223740584132967454950379795115566^{-71}{\text{,}} $$

for this test case, with $K_{0}(200)$ and $K_{1}(200)$ having $O(-88)$, and only about 2-3 decimal digits of cancellation.

[ckormanyos]

If you get a chance, John (@jzmaddock), no hurry though, I'd like your opinion on my proposed changes to precision.hpp in this PR also.

Cc: @mborland

[ckormanyos]

Remember to also compute the derivative.

Note to self:

N[D[BesselK[-129, x], x] /. {x -> 200}, 62]

$$ {\approx}-4.3110345255133348027545113739271337415489367194240537230455182^{-71} $$

[jzmaddock]

Hi @ckormanyos and Merry Christmas! Yes those changes are all good: if you wanted to be extra picky about it, you could add full specializations of log_max_value/log_min_value (in your MP code) given that those could presumably be noexcept and constexpr then? But yes, please do go ahead and merge this PR when ready!

[ckormanyos]

Thanks John (@jzmaddock) Merry Christmas also! Nothing like a few special functions of applied mathematics to raise the holida spirit!

Yes those changes are all good:

Great!

if you wanted to be extra picky about it, you could add full specializations of log_max_value/log_min_value (in your MP code) given that those could presumably be noexcept and constexpr then?

I was considering adding these within the code of the multi-precision backend. Your comment now completes my motivation to do this. I'll do it in the next round. It is too time consuming to compute the logarithms in actual run-time, so a constexpr pure number would, in fact, be superior.

But yes, please do go ahead and merge this PR when ready!

DONE!