// ******************************************************************************
   //
   // Title:       Force Field X.
   // Description: Force Field X - Software for Molecular Biophysics.
   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
   //
   // This file is part of Force Field X.
   //
   // Force Field X is free software; you can redistribute it and/or modify it
  // under the terms of the GNU General Public License version 3 as published by
  // the Free Software Foundation.
  //
  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  // details.
  //
  // You should have received a copy of the GNU General Public License along with
  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
  // Place, Suite 330, Boston, MA 02111-1307 USA
  //
  // Linking this library statically or dynamically with other modules is making a
  // combined work based on this library. Thus, the terms and conditions of the
  // GNU General Public License cover the whole combination.
  //
  // As a special exception, the copyright holders of this library give you
  // permission to link this library with independent modules to produce an
  // executable, regardless of the license terms of these independent modules, and
  // to copy and distribute the resulting executable under terms of your choice,
  // provided that you also meet, for each linked independent module, the terms
  // and conditions of the license of that module. An independent module is a
  // module which is not derived from or based on this library. If you modify this
  // library, you may extend this exception to your version of the library, but
  // you are not obligated to do so. If you do not wish to do so, delete this
  // exception statement from your version.
  //
  // ******************************************************************************
  package ffx.numerics.special;
  
  import static org.apache.commons.math3.util.FastMath.abs;
  import static org.apache.commons.math3.util.FastMath.exp;
  import static org.apache.commons.math3.util.FastMath.log;
  import static org.apache.commons.math3.util.FastMath.sqrt;
  
  /**
   * Implementation of the modified Bessel function of the first kind using Chebyshev polynomials.
   * <p>
   * Adapted from the CERN "cern.jet.math.tdouble" package as included with the
   * <a href="http://sourceforge.net/projects/parallelcolt">ParallelColt library</a>.
   */
  public class ModifiedBessel {
  
    private ModifiedBessel() {
      // Prevent instantiation.
    }
  
    /**
     * Chebyshev coefficients for exp(-x) i0(x) in the interval [0,8].
     * <br>
     * lim(x->0){ exp(-x) i0(x) } = 1.
     */
    private static final double[] A_i0 = {
        -4.41534164647933937950E-18,
        3.33079451882223809783E-17,
        -2.43127984654795469359E-16,
        1.71539128555513303061E-15,
        -1.16853328779934516808E-14,
        7.67618549860493561688E-14,
        -4.85644678311192946090E-13,
        2.95505266312963983461E-12,
        -1.72682629144155570723E-11,
        9.67580903537323691224E-11,
        -5.18979560163526290666E-10,
        2.65982372468238665035E-9,
        -1.30002500998624804212E-8,
        6.04699502254191894932E-8,
        -2.67079385394061173391E-7,
        1.11738753912010371815E-6,
        -4.41673835845875056359E-6,
        1.64484480707288970893E-5,
        -5.75419501008210370398E-5,
        1.88502885095841655729E-4,
        -5.76375574538582365885E-4,
        1.63947561694133579842E-3,
        -4.32430999505057594430E-3,
        1.05464603945949983183E-2,
        -2.37374148058994688156E-2,
        4.93052842396707084878E-2,
        -9.49010970480476444210E-2,
        1.71620901522208775349E-1,
        -3.04682672343198398683E-1,
        6.76795274409476084995E-1
    };
    /**
     * Chebyshev coefficients for exp(-x) sqrt(x) i0(x) in the inverted interval [8,infinity].
     * <br>
     * lim(x->inf){ exp(-x) sqrt(x) i0(x) } = 1/sqrt(2pi).
     */
    private static final double[] B_i0 = {
       -7.23318048787475395456E-18,
       -4.83050448594418207126E-18,
       4.46562142029675999901E-17,
       3.46122286769746109310E-17,
       -2.82762398051658348494E-16,
       -3.42548561967721913462E-16,
       1.77256013305652638360E-15,
       3.81168066935262242075E-15,
       -9.55484669882830764870E-15,
       -4.15056934728722208663E-14,
       1.54008621752140982691E-14,
       3.85277838274214270114E-13,
       7.18012445138366623367E-13,
       -1.79417853150680611778E-12,
       -1.32158118404477131188E-11,
       -3.14991652796324136454E-11,
       1.18891471078464383424E-11,
       4.94060238822496958910E-10,
       3.39623202570838634515E-9,
       2.26666899049817806459E-8,
       2.04891858946906374183E-7,
       2.89137052083475648297E-6,
       6.88975834691682398426E-5,
       3.36911647825569408990E-3,
       8.04490411014108831608E-1
   };
   /**
    * Chebyshev coefficients for exp(-x) i1(x) / x in the interval [0,8].
    * <br>
    * lim(x->0){ exp(-x) i1(x) / x } = 1/2.
    */
   private static final double[] A_i1 = {
       2.77791411276104639959E-18,
       -2.11142121435816608115E-17,
       1.55363195773620046921E-16,
       -1.10559694773538630805E-15,
       7.60068429473540693410E-15,
       -5.04218550472791168711E-14,
       3.22379336594557470981E-13,
       -1.98397439776494371520E-12,
       1.17361862988909016308E-11,
       -6.66348972350202774223E-11,
       3.62559028155211703701E-10,
       -1.88724975172282928790E-9,
       9.38153738649577178388E-9,
       -4.44505912879632808065E-8,
       2.00329475355213526229E-7,
       -8.56872026469545474066E-7,
       3.47025130813767847674E-6,
       -1.32731636560394358279E-5,
       4.78156510755005422638E-5,
       -1.61760815825896745588E-4,
       5.12285956168575772895E-4,
       -1.51357245063125314899E-3,
       4.15642294431288815669E-3,
       -1.05640848946261981558E-2,
       2.47264490306265168283E-2,
       -5.29459812080949914269E-2,
       1.02643658689847095384E-1,
       -1.76416518357834055153E-1,
       2.52587186443633654823E-1
   };
   /**
    * Chebyshev coefficients for exp(-x) sqrt(x) i1(x) in the inverted interval [8,infinity].
    * <br>
    * lim(x->inf){ exp(-x) sqrt(x) i1(x) } = 1/sqrt(2pi).
    */
   private static final double[] B_i1 = {
       7.51729631084210481353E-18,
       4.41434832307170791151E-18,
       -4.65030536848935832153E-17,
       -3.20952592199342395980E-17,
       2.96262899764595013876E-16,
       3.30820231092092828324E-16,
       -1.88035477551078244854E-15,
       -3.81440307243700780478E-15,
       1.04202769841288027642E-14,
       4.27244001671195135429E-14,
       -2.10154184277266431302E-14,
       -4.08355111109219731823E-13,
       -7.19855177624590851209E-13,
       2.03562854414708950722E-12,
       1.41258074366137813316E-11,
       3.25260358301548823856E-11,
       -1.89749581235054123450E-11,
       -5.58974346219658380687E-10,
       -3.83538038596423702205E-9,
       -2.63146884688951950684E-8,
       -2.51223623787020892529E-7,
       -3.88256480887769039346E-6,
       -1.10588938762623716291E-4,
       -9.76109749136146840777E-3,
       7.78576235018280120474E-1
   };
 
   /**
    * Modified zero-order Bessel function.
    * <p>
    * The function is defined as i0(x) = J0( ix ). The range is partitioned into the two intervals
    * [0,8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval.
    *
    * @param x input parameter
    * @return i0(x)
    */
   public static double i0(double x) {
     // Handle special cases
     if (Double.isNaN(x)) {
       return Double.NaN;
     }
     if (Double.isInfinite(x)) {
       return 0.0;
     }
 
     double y;
     if (x < 0) {
       x = -x;
     }
     if (x <= 8.0) {
       y = (x * 0.5) - 2.0;
       return (eToThe(x) * evaluateChebyshev(y, A_i0, 30));
     }
     double ix = 1.0 / x;
     return (eToThe(x) * evaluateChebyshev(32.0 * ix - 2.0, B_i0, 25) * sqrt(ix));
   }
 
   /**
    * Modified 1st-order Bessel function.
    * <p>
    * The function is defined as i1(x) = -i j1( ix ). The range is partitioned into the two intervals
    * [0,8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval.
    *
    * @param x input parameter.
    * @return i1(x).
    */
   public static double i1(double x) {
     // Handle special cases
     if (Double.isNaN(x)) {
       return Double.NaN;
     }
     if (Double.isInfinite(x)) {
       return 0.0;
     }
 
     double y, z;
 
     z = abs(x);
     if (z <= 8.0) {
       y = (z * 0.5) - 2.0;
       z = evaluateChebyshev(y, A_i1, 29) * z * eToThe(z);
     } else {
       double iz = 1.0 / z;
       z = eToThe(z) * evaluateChebyshev(32.0 * iz - 2.0, B_i1, 25) * sqrt(iz);
     }
     if (x < 0.0) {
       z = -z;
     }
     return (z);
   }
 
   /**
    * Compute the ratio of i1(x) to i0(x).
    * <p>
    * This implementation avoids redundant calculations by computing both i1(x) and i0(x)
    * efficiently for different ranges of x.
    *
    * @param x input parameter
    * @return i1(x) / i0(x)
    */
   public static double i1OverI0(double x) {
     // Handle special cases
     if (Double.isNaN(x)) {
       return Double.NaN;
     }
     if (Double.isInfinite(x)) {
       return 0.0; // Both i1(±∞) and i0(±∞) are 0, but i1/i0 approaches 0 as x approaches ±∞
     }
     if (x == 0.0) {
       return 0.0; // i1(0) = 0, i0(0) = 1
     }
 
     // For small values, use the direct calculation to avoid precision loss
     if (abs(x) < 1.0) {
       return i1(x) / i0(x);
     }
 
     double absX = abs(x);
     double result;
 
     if (absX <= 8.0) {
       // For moderate values, calculate both functions with shared operations
       double y = (absX * 0.5) - 2.0;
       double expX = eToThe(absX);
       double i0Val = expX * evaluateChebyshev(y, A_i0, 30);
       double i1Val = evaluateChebyshev(y, A_i1, 29) * absX * expX;
       result = i1Val / i0Val;
     } else {
       // For large values, use the asymptotic behavior
       double ix = 1.0 / absX;
       double expX = eToThe(absX);
       double sqrtIx = sqrt(ix);
       double i0Val = expX * evaluateChebyshev(32.0 * ix - 2.0, B_i0, 25) * sqrtIx;
       double i1Val = expX * evaluateChebyshev(32.0 * ix - 2.0, B_i1, 25) * sqrtIx;
       result = i1Val / i0Val;
     }
 
     // Adjust sign for negative x
     return (x < 0.0) ? -result : result;
   }
 
   /**
    * Compute the natural log(i0(x)).
    * <p>
    * This implementation computes log(i0(x)) directly to avoid overflow for large x values.
    * For large x, i0(x) ~ exp(x)/sqrt(2*pi*x), so log(i0(x)) ~ x - 0.5*log(2*pi*x).
    *
    * @param x input parameter.
    * @return the natural log(i0(x)).
    */
   public static double lnI0(double x) {
     // Handle special cases
     if (Double.isNaN(x)) {
       return Double.NaN;
     }
     if (Double.isInfinite(x)) {
       return Double.NEGATIVE_INFINITY; // log(0) = -∞
     }
 
     // For small values, use the direct calculation to avoid precision loss
     if (abs(x) <= 8.0) {
       double i0Val = i0(x);
       // Avoid log(0) which would return -Infinity
       return i0Val > 0.0 ? log(i0Val) : Double.NEGATIVE_INFINITY;
     }
 
     // For large values, use asymptotic formula to avoid overflow
     double absX = abs(x);
     double ix = 1.0 / absX;
 
     // Compute log(i0(x)) directly using the asymptotic formula
     // log(i0(x)) = x + log(evaluateChebyshev(...) * sqrt(ix))
     double chebyshevTerm = evaluateChebyshev(32.0 * ix - 2.0, B_i0, 25);
 
     // Handle potential underflow in the Chebyshev evaluation
     if (chebyshevTerm <= 0.0) {
       return absX; // Fallback to the dominant term
     }
 
     return absX + log(chebyshevTerm * sqrt(ix));
   }
 
   /**
    * Evaluates Chebyshev polynomials (first kind) at x/2.
    * <p>
    * NOTE: Argument x must first be transformed to the interval (-1,1).
    * <p>
    * NOTE: Coefficients are in reverse; zero-order term is last.
    *
    * @param x            argument to the polynomial.
    * @param coefficients the coefficients of the polynomial.
    * @param N            the number of coefficients.
    * @return the result
    */
   private static double evaluateChebyshev(double x, double[] coefficients, int N) {
     double b0, b1, b2;
 
     int p = 0;
     int i;
 
     b0 = coefficients[p++];
     b1 = 0.0;
     i = N - 1;
 
     do {
       b2 = b1;
       b1 = b0;
       b0 = x * b1 - b2 + coefficients[p++];
     } while (--i > 0);
 
     return (0.5 * (b0 - b2));
   }
 
   /**
    * Computes exp(x) with protection against overflow and underflow.
    * <p>
    * Returns Double.MAX_VALUE in place of Double.POSITIVE_INFINITY and
    * Double.MIN_VALUE in place of Double.NEGATIVE_INFINITY. This prevents
    * arithmetic operations from generating NaN results when overflow occurs.
    *
    * @param x input parameter
    * @return exp(x) with overflow/underflow protection
    */
   private static double eToThe(double x) {
     // Handle NaN input
     if (Double.isNaN(x)) {
       return Double.NaN;
     }
 
     // Handle extreme values directly to avoid unnecessary computation
     if (x > 709.0) { // Approximate threshold where exp(x) overflows to Infinity
       return Double.MAX_VALUE;
     }
     if (x < -745.0) { // Approximate threshold where exp(x) underflows to 0
       return Double.MIN_VALUE;
     }
 
     double res = exp(x);
     if (res == Double.POSITIVE_INFINITY) {
       return Double.MAX_VALUE;
     } else if (res == Double.NEGATIVE_INFINITY || res == 0.0) {
       return Double.MIN_VALUE;
     }
     return res;
   }
 }