ENH: add iv_ratioinv

include/xsf/iv_ratio.h

@@ -1,10 +1,11 @@
-// Numerically stable computation of iv(v+1, x) / iv(v, x)
+// Numerically stable computation of iv(v, x) / iv(v - 1, x) and its inverse
 
 #pragma once
 
 #include "cephes/dd_real.h"
 #include "config.h"
 #include "error.h"
+#include "log_exp.h"
 #include "tools.h"
 
 namespace xsf {
@@ -98,7 +99,10 @@ XSF_HOST_DEVICE inline double iv_ratio(double v, double x) {
     if (std::isinf(x)) {
         return 1.0;
     }
-
+    if (v == 0.5) {
+        // Closed-form solution for v = 0.5: iv_ratio(0.5, x) = tanh(x)
+        return std::tanh(x);
+    }
     auto [ret, terms] = _iv_ratio_cf<double>(v, x, false);
     if (terms == 0) { // failed to converge; should not happen
         set_error("iv_ratio", SF_ERROR_NO_RESULT, NULL);
@@ -157,13 +161,117 @@ XSF_HOST_DEVICE inline double iv_ratio_c(double v, double x) {
     } else {
         // The previous branch (v > 0.5) also works for v == 0.5, but
         // the closed-form formula "1 - tanh(x)" is more efficient.
-        double t = std::exp(-2 * x);
-        return (2 * t) / (1 + t);
+        // Equivalently: iv_ratio_c(0.5, x) = 1 - tanh(x) = 2 * expit(-2*x)
+        return 2 * expit(-2 * x);
     }
 }
 
 XSF_HOST_DEVICE inline float iv_ratio_c(float v, float x) {
     return iv_ratio_c(static_cast<double>(v), static_cast<double>(x));
 }
 
+XSF_HOST_DEVICE inline double iv_ratioinv(double v, double r) {
+    if (std::isnan(v) || std::isnan(r)) {
+        return std::numeric_limits<double>::quiet_NaN();
+    }
+    if (!std::isfinite(v) || v < 0.5 || r < 0.0 || r > 1.0) {
+        // iv_ratioinv is only defined for v >= 0.5
+        set_error("iv_ratioinv", SF_ERROR_DOMAIN, NULL);
+        return std::numeric_limits<double>::quiet_NaN();
+    }
+    if (r == 0.0) {
+        return 0.0;
+    }
+    if (r == 1.0) {
+        return std::numeric_limits<double>::infinity();
+    }
+    if (v == 0.5) {
+        // Closed-form solution for v = 0.5: iv_ratio(0.5, x) = tanh(x)
+        // Since tanh(x) = 2*expit(2*x) - 1 = r, we have expit(2*x) = (1+r)/2
+        // For r > 0.5, use logit((1+r)/2) for better stability near r=1
+        // For r <= 0.5, atanh is sufficiently stable
+        if (r > 0.5) {
+            return 0.5 * logit((1.0 + r) * 0.5);
+        } else {
+            return std::atanh(r);
+        }
+    }
+    // Algorithm description: use Chandrupatla's method to find the root of
+    // f(x) = iv_ratio(v, x) - r (or f(x) = 1 - iv_ratio_c(v, x) - r if r > 0.5).
+    // Bounds from Amos, 1973. "Computation of Modified Bessel Function And Their Ratios".
+    // That paper defines the ratio as iv(v+1, x)/iv(v, x), but we define it as iv(v, x)/iv(v-1, x).
+    // Substitute v -> v-1 in eq. (9) and eq. (16) yields the following bounds for x:
+    //   x/(v-1/2 + sqrt(x^2 + (v+1/2)^2)) <= r <= x/(v-1 + sqrt(x^2 + (v+1)^2))
+    //
+    // Inverting (inequality direction flips):
+    //   lower x (from upper r-bound): r*[(v-1) + sqrt((v-1)^2 + 4v*(1-r^2))] / (1-r^2)
+    //   upper x (from lower r-bound): r*[(v-0.5) + sqrt((v-0.5)^2 + 2v*(1-r^2))] / (1-r^2)
+
+    double r_c = 1.0 - r;
+    double one_minus_r_sq = r_c * (1.0 + r); // 1 - r^2
+
+    double vm1 = v - 1.0;
+    double lower_bound = r * (vm1 + std::sqrt(vm1 * vm1 + 4.0 * v * one_minus_r_sq)) / one_minus_r_sq;
+
+    double vm05 = v - 0.5;
+    double upper_bound = r * (vm05 + std::sqrt(vm05 * vm05 + 2.0 * v * one_minus_r_sq)) / one_minus_r_sq;
+
+    // For small r both bounds converge to 2v*r (== true x). Ensure bracket has width.
+    if (upper_bound <= lower_bound) {
+        upper_bound = 2.0 * lower_bound;
+    }
+
+    auto func = [v, r, r_c](double x) {
+        if (r <= 0.5) {
+            return iv_ratio(v, x) - r;
+        }
+        return r_c - iv_ratio_c(v, x);
+    };
+
+    double xl = lower_bound;
+    double xr = upper_bound;
+    double f_xl = func(xl);
+    double f_xr = func(xr);
+
+    if (f_xl * f_xr > 0) {
+        // The Amos bounds are theoretically correct but may yield an invalid bracket
+        // due to precision issues in iv_ratio or iv_ratio_c. This happens especially for very small or large r.
+        // Fallback: use bracket_root_for_cdf_inversion to find a valid bracket.
+        // Bracketing parameters taken from gdtrib, only difference: our function is increasing
+        auto [b_xl, b_xr, b_f_xl, b_f_xr, bracket_status] = detail::bracket_root_for_cdf_inversion(
+            func, 1.0, std::numeric_limits<double>::min(), std::numeric_limits<double>::max(), -0.875, 7.0, 0.125, 8,
+            true, 342
+        );
+        if (bracket_status == 1) {
+            set_error("iv_ratioinv", SF_ERROR_UNDERFLOW, NULL);
+            return 0.0;
+        }
+        if (bracket_status == 2) {
+            set_error("iv_ratioinv", SF_ERROR_OVERFLOW, NULL);
+            return std::numeric_limits<double>::infinity();
+        }
+        if (bracket_status >= 3) {
+            set_error("iv_ratioinv", SF_ERROR_OTHER, "Computational Error");
+            return std::numeric_limits<double>::quiet_NaN();
+        }
+        xl = b_xl;
+        xr = b_xr;
+        f_xl = b_f_xl;
+        f_xr = b_f_xr;
+    }
+
+    auto [result, root_status] =
+        detail::find_root_chandrupatla(func, xl, xr, f_xl, f_xr, std::numeric_limits<double>::epsilon(), 1e-308, 1000);
+    if (root_status) {
+        // Root finding failed. This should never happen.
+        set_error("iv_ratioinv", SF_ERROR_OTHER, "Computational Error");
+        return std::numeric_limits<double>::quiet_NaN();
+    }
+    return result;
+}
+
+XSF_HOST_DEVICE inline float iv_ratioinv(float v, float r) {
+    return static_cast<float>(iv_ratioinv(static_cast<double>(v), static_cast<double>(r)));
+}
+
 } // namespace xsf

tests/xsf_tests/test_iv_ratioinv.cpp

@@ -0,0 +1,99 @@
+#include "../testing_utils.h"
+#include "xsf/config.h"
+#include <xsf/iv_ratio.h>
+
+TEST_CASE("iv_ratioinv round-trip single-double", "[iv_ratioinv][xsf_tests]") {
+    const std::vector<double> xs = logspace<double>(-10, 7, 100);
+    const std::vector<double> vs = logspace<double>(0, 6, 100);
+    for (double x : xs) {
+        for (double v : vs) {
+            const double y_forward = xsf::iv_ratio(v, x);
+            const double y = (y_forward < 0.5) ? y_forward : (1.0 - xsf::iv_ratio_c(v, x));
+            if ((y == 0.) || (y == 1.)) {
+                continue;
+            }
+            const auto result = xsf::iv_ratioinv(v, y);
+            const auto relative_error = xsf::extended_relative_error(result, x);
+            double rtol = 1e-9;
+            CAPTURE(v, x, y, result, rtol, relative_error);
+            REQUIRE(relative_error <= rtol);
+
+            const float v_f = static_cast<float>(v);
+            const float x_f = static_cast<float>(x);
+            const float y_forward_f = xsf::iv_ratio(v_f, x_f);
+            const float y_f = (y_forward_f < 0.5f) ? y_forward_f : (1.0f - xsf::iv_ratio_c(v_f, x_f));
+            if ((y_f < 0.0001f) || (y_f > 0.9999f)) {
+                // Skip the tails in single precision because precision is not sufficient
+                continue;
+            }
+            const auto result_f = xsf::iv_ratioinv(v_f, y_f);
+            const auto relative_error_f = xsf::extended_relative_error(result_f, x_f);
+            float rtol_f = 5e-4f;
+            CAPTURE(v, x, y_f, result_f, rtol_f, relative_error_f);
+            REQUIRE(relative_error_f <= rtol_f);
+        }
+    }
+}
+
+TEST_CASE("iv_ratioinv arbitrary precision", "[iv_ratioinv][xsf_tests]") {
+    // Test with arbitrary precision values for v, x, and iv_ratio(v,x)
+    // Reference values for iv_ratio were computed with the Python library mpmath.
+    // from mpmath import mp
+    // mp.dps = 1000
+    // def iv_ratio(v, x):
+    //     v = mp.mpf(v)
+    //     x = mp.mpf(x)
+    //     result = mp.besseli(v, x) / mp.besseli(v - 1, x)
+    //     return float(result)
+
+    using test_case = std::tuple<double, double, double>;
+
+    auto [v, x, iv_ratio_ref] = GENERATE(
+        test_case{10, 100000, 0.9999050040375403}, test_case{1, 0.1, 0.04993760398793892},
+        test_case{0.5, 0.001, 0.0009999996666668}, test_case{1000, 500, 0.23607911616885813},
+        test_case{100, 1e-15, 5e-18}, test_case{1e5, 5, 2.4999999984375157e-05},
+        test_case{1e5, 1e5, 0.41421399130458997}, test_case{1e5, 1e10, 0.999990000099999}
+    );
+    double x_result = xsf::iv_ratioinv(v, iv_ratio_ref);
+    const auto rel_error = xsf::extended_relative_error(x_result, x);
+
+    CAPTURE(v, x, iv_ratio_ref, x_result, rel_error);
+    REQUIRE(rel_error <= 1e-10);
+}
+
+TEST_CASE("iv_ratioinv nan propagation", "[iv_ratioinv][xsf_tests]") {
+    const double nan = std::numeric_limits<double>::quiet_NaN();
+    REQUIRE(std::isnan(xsf::iv_ratioinv(nan, 0.5)));
+    REQUIRE(std::isnan(xsf::iv_ratioinv(1.0, nan)));
+}
+
+TEST_CASE("iv_ratioinv domain error", "[iv_ratioinv][xsf_tests]") {
+    const double inf = std::numeric_limits<double>::infinity();
+    REQUIRE(std::isnan(xsf::iv_ratioinv(0.4, 0.5)));  // v < 0.5
+    REQUIRE(std::isnan(xsf::iv_ratioinv(1.0, -1.0))); // y < 0
+    REQUIRE(std::isnan(xsf::iv_ratioinv(1.0, 2.0)));  // y > 1
+    REQUIRE(std::isnan(xsf::iv_ratioinv(inf, 0.5)));  // v = inf
+}
+
+TEST_CASE("iv_ratioinv domain boundary", "[iv_ratioinv][xsf_tests]") {
+    const double inf = std::numeric_limits<double>::infinity();
+    REQUIRE(xsf::iv_ratioinv(1.0, 0.0) == 0.0); // y = 0
+    REQUIRE(xsf::iv_ratioinv(1.0, 1.0) == inf); // y = 1
+}
+
+TEST_CASE("iv_ratioinv v = 0.5 roundtrip", "[iv_ratioinv][xsf_tests]") {
+    const std::vector<double> xs = logspace<double>(-50, 2, 500);
+    double threshold = 1.0 - 1e-12; // threshold to avoid numerical issues near y=1
+    for (double x : xs) {
+        const double y_forward = xsf::iv_ratio(0.5, x);
+        const double y = (y_forward < 0.5) ? y_forward : (1.0 - xsf::iv_ratio_c(0.5, x));
+        if ((y == 0.) || (y > threshold)) {
+            continue;
+        }
+        const double x_result = xsf::iv_ratioinv(0.5, y);
+        const auto relative_error = xsf::extended_relative_error(x_result, x);
+        double rtol = 1e-7;
+        CAPTURE(x, y, x_result, rtol, relative_error);
+        REQUIRE(relative_error <= rtol);
+    }
+}