fixed case nu=x=y=0

Author: klausspanderen

ql/math/modifiedbessel.cpp

@@ -101,6 +101,9 @@ namespace QuantLib {
 
     std::complex<Real> modifiedBesselFunction_i(Real nu,
                                                 const std::complex<Real> &z) {
+        if (z.imag() == 0.0 && z.real() >= 0.0)
+            return std::complex<Real>(modifiedBesselFunction_i(nu, z.real()));
+
         return modifiedBesselFunction_i_impl<
             std::complex<Real>, Unweighted>(nu, z);
     }
@@ -111,6 +114,9 @@ namespace QuantLib {
 
     std::complex<Real> modifiedBesselFunction_k(Real nu,
                                                 const std::complex<Real> &z) {
+        if (z.imag() == 0.0 && z.real() >= 0.0)
+            return std::complex<Real>(modifiedBesselFunction_k(nu, z.real()));
+
         return modifiedBesselFunction_k_impl<
             std::complex<Real>, Unweighted>(nu, z);
     }
@@ -124,6 +130,11 @@ namespace QuantLib {
 
     std::complex<Real> modifiedBesselFunction_i_exponentiallyWeighted(
         Real nu, const std::complex<Real> &z) {
+
+        if (z.imag() == 0.0 && z.real() >= 0.0)
+            return std::complex<Real>(
+                modifiedBesselFunction_i_exponentiallyWeighted(nu, z.real()));
+
         return modifiedBesselFunction_i_impl<
             std::complex<Real>, ExponentiallyWeighted>(nu, z);
     }
@@ -135,6 +146,11 @@ namespace QuantLib {
 
     std::complex<Real> modifiedBesselFunction_k_exponentiallyWeighted(
         Real nu, const std::complex<Real> &z) {
+
+        if (z.imag() == 0.0 && z.real() >= 0.0)
+            return std::complex<Real>(
+                modifiedBesselFunction_k_exponentiallyWeighted(nu, z.real()));
+
         return modifiedBesselFunction_k_impl<
             std::complex<Real>, ExponentiallyWeighted>(nu, z);
     }

test-suite/functions.cpp

@@ -259,8 +259,6 @@ void FunctionsTest::testWeightedModifiedBesselFunctions() {
     for (Real nu = -5.0; nu <= 5.0; nu += 0.5) {
         for (Real x = -5.0; x <= 5.0; x += 0.5) {
             for (Real y = -5.0; y <= 5.0; y += 0.5) {
-                if (close_enough(x, 0.0) && close_enough(y, 0.0))
-                    continue;
                 std::complex<Real> z(x, y);
                 std::complex<Real> calculated_i =
                     modifiedBesselFunction_i_exponentiallyWeighted(nu, z);