I am fairly new to Julia and have been working with ArbNumerics.jl BesselK function as part of evaluating a complicated response function that can be expanded in a series of BesselK functions. The arguments it ends up taking sometimes cause it to underflow at 64-bit precision. Although the function decays to very small values, these cancel with other parts of my equation that are large to give finite values between 0-1000. So, as long as I extend the precision of the BesselK functions, everything comes out nicely. However, I have noticed that sometimes ArbNumerics implementation gives the wrong value for certain arguments for a range of precisions. For example:
setextrabits(0) setworkingprecision(ArbFloat, bits = p) bessel = ArbNumerics.besselk(ArbFloat(2.0), ArbFloat(57.966855579194714))
is one term in my code that returns the wrong value for precisions p between 116 and 175 bits which ultimately gives unphysical results.
I admit I do not fully understand what is going on here. How do I deal with this numerically volatile expression? Is there a way to correct my usage of ArbNumerics.jl, or any suggestions to tame these extreme floats?