Fix modifiedBesselFunction_k at integer order and for large arguments

[Abhayindia]

Fixes #2785.

modifiedBesselFunction_k is built from pi/2 * (I_(-nu) - I_nu) / sin(nu*pi), which is correct analytically but fails in the two regimes described in the issue. This replaces the evaluation scheme in those two regimes and leaves it untouched elsewhere.

Large argument

The issue attributes the x >= 13 zeros to cancellation. It is worth being precise, because the mechanism is stronger than that and explains why the failure is abrupt rather than gradual: for |x| >= 13, modifiedBesselFunction_i_impl takes its asymptotic branch, whose coefficients depend on nu only through nu*nu, and I<Real>::value() is 0.0 so the i*exp(i*nu*pi) term drops out. I_(-nu) and I_nu are therefore computed from identical arithmetic and the difference is exactly zero by construction, not merely ill-conditioned.

Note this is specific to real arguments. For std::complex that term is i rather than 0, the two values genuinely differ, and the subtraction remains usable, which is why the complex overloads never produced the silent zero.

K is now summed directly from its own asymptotic series, which shares the coefficients already built in the I branch but without the alternating sign:

K_nu(x) ~ sqrt(pi/2x) e^-x [ 1 + (mu-1)/(8x) + (mu-1)(mu-9)/(2!(8x)^2) + ... ],  mu = 4nu^2

The existing weight2LargeX covers both overloads: exp(-x) unweighted and exp(-2x) weighted, matching the header's f(x)*exp(-x) contract.

The switch stays at 13, matching the I branch. Lowering it to 10 measurably improves real arguments in [10, 13) (at x = 12, about 3e-12 against 4e-7 for the subtraction) but costs accuracy for complex arguments near the imaginary axis, where this series is less converged than the difference it would replace. That trade seemed worth leaving out of a bug fix.

Integer order

K_n is computed from the integer-order series, A&S 9.6.11, when |nu - round(nu)| is below a threshold. K is even in nu, so negative integers use K_|n|.

The issue notes this needs digamma. In fact only psi at positive integers is required, where psi(1) = -gamma and psi(m) = -gamma + sum_{j<m} 1/j, so no general digamma function is introduced.

The threshold is 1e-9, chosen so the change never makes anything worse. Measured against boost::math::cyl_bessel_k at nu = 1 + delta, the crossover where the integer series becomes the more accurate of the two is about 2.4e-9 at x = 0.5 and larger for larger x; sitting below the smallest crossover means the series is only used where it wins. This fixes the removable singularity and its immediate neighbourhood. It does not address the gradual near-integer loss the issue also measures at nu = 0.001, which needs an evaluation uniform in nu rather than a threshold.

Tests

testModifiedBesselFunctions and testWeightedModifiedBesselFunctions both pass unchanged in substance, but the two weighted checks needed correcting: they computed their expected value by inlining the same I_(-nu) - I_nu quotient, over loops stepping nu from -5 to 5 in halves, so at every integer order they compared the broken formula against itself and passed. The first-kind half of those same checks already compares against the public modifiedBesselFunction_i; the second-kind half now does the same.

This may be relevant to #2432. If a test's reference is the implementation it is testing, a failure there is not evidence about tolerance.

testModifiedBesselFunctionsRegressionCases is new: eight points from boost::math, four at integer order and four past the old cliff. It is verified to fail on master with exactly the values in the issue (-2.04e+16 at nu = 1, 0 at x = 13). Its second-kind tolerance scales by |K| rather than by max(|I|, |K|) as the existing block does, since out there I is around 1e11 while K is around 1e-14 and a bound taken from the larger accepts anything.

Verification

Reproducer from the issue:

	before	after
k(0, 2)	nan	2.4e-16
k(1, 2)	-2.04e+16	exact
k(0.5, 14)	0	exact

Sweep against boost::math::cyl_bessel_k, 144 real points over nu in [-10, 10] and x in [0.1, 40]: worst relative error 1.5e-07, at x = 8, which is below the switch and unchanged by this PR. Complex overload on the positive real axis: worst 3.5e-11. Full test suite passes.

One limitation worth recording: the asymptotic series needs x large relative to nu, so for nu around 30 at x near 13 it is only good to about 1e-6. That is inherent to the expansion and the existing I branch shares it, but previously those points returned zero.

[coveralls]

coverage: 75.644% (+0.01%) from 75.63% — Abhayindia:fix/2785-bessel-k-evaluation into lballabio:master