Merge pull request #1456.

fixed modified Bessel function for nu=x=y=0 case.
lballabio · Aug 16, 2022 · 9c2c62e · 9c2c62e
2 parents 0caf57a + db1e85b
commit 9c2c62e
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 2 deletions.
diff --git a/ql/math/modifiedbessel.cpp b/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);
     }

diff --git a/test-suite/functions.cpp b/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);