BesselJ.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. package org.apache.commons.math4.legacy.special;

  20. import java.util.Arrays;
  21. import org.apache.commons.numbers.gamma.Gamma;
  22. import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
  23. import org.apache.commons.math4.legacy.exception.ConvergenceException;
  24. import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
  25. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  26. import org.apache.commons.math4.core.jdkmath.JdkMath;

  28. /**
  29.  * This class provides computation methods related to Bessel
  30.  * functions of the first kind. Detailed descriptions of these functions are
  31.  * available in <a
  32.  * href="http://en.wikipedia.org/wiki/Bessel_function">Wikipedia</a>, <a
  33.  * href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">Abrabowitz and
  34.  * Stegun</a> (Ch. 9-11), and <a href="http://dlmf.nist.gov/">DLMF</a> (Ch. 10).
  35.  * <p>
  36.  * This implementation is based on the rjbesl Fortran routine at
  37.  * <a href="http://www.netlib.org/specfun/rjbesl">Netlib</a>.</p>
  38.  * <p>
  39.  * From the Fortran code: </p>
  40.  * <p>
  41.  * This program is based on a program written by David J. Sookne (2) that
  42.  * computes values of the Bessel functions J or I of real argument and integer
  43.  * order. Modifications include the restriction of the computation to the J
  44.  * Bessel function of non-negative real argument, the extension of the
  45.  * computation to arbitrary positive order, and the elimination of most
  46.  * underflow.</p>
  47.  * <p>
  48.  * References:</p>
  49.  * <ul>
  50.  * <li>"A Note on Backward Recurrence Algorithms," Olver, F. W. J., and Sookne,
  51.  * D. J., Math. Comp. 26, 1972, pp 941-947.</li>
  52.  * <li>"Bessel Functions of Real Argument and Integer Order," Sookne, D. J., NBS
  53.  * Jour. of Res. B. 77B, 1973, pp 125-132.</li>
  54.  * </ul>
  55.  * @since 3.4
  56.  */
  57. public class BesselJ
  58.     implements UnivariateFunction {

  60.     // ---------------------------------------------------------------------
  61.     // Mathematical constants
  62.     // ---------------------------------------------------------------------

  64.     /** -2 / pi. */
  65.     private static final double PI2 = 0.636619772367581343075535;

  67.     /** first few significant digits of 2pi. */
  68.     private static final double TOWPI1 = 6.28125;

  70.     /** 2pi - TWOPI1 to working precision. */
  71.     private static final double TWOPI2 = 1.935307179586476925286767e-3;

  73.     /** TOWPI1 + TWOPI2. */
  74.     private static final double TWOPI = TOWPI1 + TWOPI2;

  76.     // ---------------------------------------------------------------------
  77.     // Machine-dependent parameters
  78.     // ---------------------------------------------------------------------

  80.     /**
  81.      * 10.0^K, where K is the largest integer such that ENTEN is
  82.      * machine-representable in working precision.
  83.      */
  84.     private static final double ENTEN = 1.0e308;

  86.     /**
  87.      * Decimal significance desired. Should be set to (INT(log_{10}(2) * (it)+1)).
  88.      * Setting NSIG lower will result in decreased accuracy while setting
  89.      * NSIG higher will increase CPU time without increasing accuracy.
  90.      * The truncation error is limited to a relative error of
  91.      * T=.5(10^(-NSIG)).
  92.      */
  93.     private static final double ENSIG = 1.0e16;

  95.     /** 10.0 ** (-K) for the smallest integer K such that K >= NSIG/4. */
  96.     private static final double RTNSIG = 1.0e-4;

  98.     /** Smallest ABS(X) such that X/4 does not underflow. */
  99.     private static final double ENMTEN = 8.90e-308;

  101.     /** Minimum acceptable value for x. */
  102.     private static final double X_MIN = 0.0;

  104.     /**
  105.      * Upper limit on the magnitude of x. If abs(x) = n, then at least
  106.      * n iterations of the backward recursion will be executed. The value of
  107.      * 10.0 ** 4 is used on every machine.
  108.      */
  109.     private static final double X_MAX = 1.0e4;

  111.     /** First 25 factorials as doubles. */
  112.     private static final double[] FACT = {
  113.         1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
  114.         3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
  115.         1.307674368e12, 2.0922789888e13, 3.55687428096e14, 6.402373705728e15,
  116.         1.21645100408832e17, 2.43290200817664e18, 5.109094217170944e19,
  117.         1.12400072777760768e21, 2.585201673888497664e22,
  118.         6.2044840173323943936e23
  119.     };

  121.     /** Order of the function computed when {@link #value(double)} is used. */
  122.     private final double order;

  124.     /**
  125.      * Create a new BesselJ with the given order.
  126.      *
  127.      * @param order order of the function computed when using {@link #value(double)}.
  128.      */
  129.     public BesselJ(double order) {
  130.         this.order = order;
  131.     }

  133.     /**
  134.      * Returns the value of the constructed Bessel function of the first kind,
  135.      * for the passed argument.
  136.      *
  137.      * @param x Argument
  138.      * @return Value of the Bessel function at x
  139.      * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order}
  140.      * @throws ConvergenceException if the algorithm fails to converge
  141.      */
  142.     @Override
  143.     public double value(double x)
  144.         throws MathIllegalArgumentException, ConvergenceException {
  145.         return BesselJ.value(order, x);
  146.     }

  148.     /**
  149.      * Returns the first Bessel function, \(J_{order}(x)\).
  150.      *
  151.      * @param order Order of the Bessel function
  152.      * @param x Argument
  153.      * @return Value of the Bessel function of the first kind, \(J_{order}(x)\)
  154.      * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order}
  155.      * @throws ConvergenceException if the algorithm fails to converge
  156.      */
  157.     public static double value(double order, double x)
  158.         throws MathIllegalArgumentException, ConvergenceException {
  159.         final int n = (int) order;
  160.         final double alpha = order - n;
  161.         final int nb = n + 1;
  162.         final BesselJResult res = rjBesl(x, alpha, nb);

  164.         if (res.nVals >= nb) {
  165.             return res.vals[n];
  166.         } else if (res.nVals < 0) {
  167.             throw new MathIllegalArgumentException(LocalizedFormats.BESSEL_FUNCTION_BAD_ARGUMENT,order, x);
  168.         } else if (JdkMath.abs(res.vals[res.nVals - 1]) < 1e-100) {
  169.             return res.vals[n]; // underflow; return value (will be zero)
  170.         }
  171.         throw new ConvergenceException(LocalizedFormats.BESSEL_FUNCTION_FAILED_CONVERGENCE, order, x);
  172.     }

  174.     /**
  175.      * Encapsulates the results returned by {@link BesselJ#rjBesl(double, double, int)}.
  176.      * <p>
  177.      * {@link #getVals()} returns the computed function values.
  178.      * {@link #getnVals()} is the number of values among those returned by {@link #getnVals()}
  179.      * that can be considered accurate.
  180.      * </p><ul>
  181.      * <li>{@code nVals < 0}: An argument is out of range. For example, {@code nb <= 0},
  182.      * {@code alpha < 0 or > 1}, or x is too large. In this case, b(0) is set to zero, the
  183.      * remainder of the b-vector is not calculated, and nVals is set to
  184.      * MIN(nb,0) - 1 so that nVals != nb.</li>
  185.      * <li>{@code nb > nVals > 0}: Not all requested function values could be calculated
  186.      * accurately. This usually occurs because nb is much larger than abs(x). In
  187.      * this case, b(n) is calculated to the desired accuracy for {@code n < nVals}, but
  188.      * precision is lost for {@code nVals < n <= nb}. If b(n) does not vanish for
  189.      * {@code n > nVals} (because it is too small to be represented), and b(n)/b(nVals) =
  190.      * \(10^{-k}\), then only the first NSIG-k significant figures of b(n) can be
  191.      * trusted.</li></ul>
  192.      */
  193.     public static class BesselJResult {

  195.         /** Bessel function values. */
  196.         private final double[] vals;

  198.         /** Valid value count. */
  199.         private final int nVals;

  201.         /**
  202.          * Create a new BesselJResult with the given values and valid value count.
  203.          *
  204.          * @param b values
  205.          * @param n count of valid values
  206.          */
  207.         public BesselJResult(double[] b, int n) {
  208.             vals = Arrays.copyOf(b, b.length);
  209.             nVals = n;
  210.         }

  212.         /**
  213.          * @return the computed function values
  214.          */
  215.         public double[] getVals() {
  216.             return Arrays.copyOf(vals, vals.length);
  217.         }

  219.         /**
  220.          * @return the number of valid function values (normally the same as the
  221.          *         length of the array returned by {@link #getnVals()})
  222.          */
  223.         public int getnVals() {
  224.             return nVals;
  225.         }
  226.     }

  228.     /**
  229.      * Calculates Bessel functions \(J_{n+alpha}(x)\) for
  230.      * non-negative argument x, and non-negative order n + alpha.
  231.      * <p>
  232.      * Before using the output vector, the user should check that
  233.      * nVals = nb, i.e., all orders have been calculated to the desired accuracy.
  234.      * See BesselResult class javadoc for details on return values.
  235.      * </p>
  236.      * @param x non-negative real argument for which J's are to be calculated
  237.      * @param alpha fractional part of order for which J's or exponentially
  238.      * scaled J's (\(J\cdot e^{x}\)) are to be calculated. {@code 0 <= alpha < 1.0}
  239.      * @param nb integer number of functions to be calculated, {@code nb > 0}. The first
  240.      * function calculated is of order alpha, and the last is of order
  241.      * {@code nb - 1 + alpha}.
  242.      * @return BesselJResult a vector of the functions
  243.      * \(J_{alpha}(x)\) through \(J_{nb-1+alpha}(x)\), or the corresponding exponentially
  244.      * scaled functions and an integer output variable indicating possible errors
  245.      */
  246.     public static BesselJResult rjBesl(double x, double alpha, int nb) {
  247.         final double[] b = new double[nb];

  249.         int ncalc = 0;
  250.         double alpem = 0;
  251.         double alp2em = 0;

  253.         // ---------------------------------------------------------------------
  254.         // Check for out of range arguments.
  255.         // ---------------------------------------------------------------------
  256.         final int magx = (int) x;
  257.         if (nb > 0 && x >= X_MIN && x <= X_MAX && alpha >= 0 && alpha < 1) {
  258.             // ---------------------------------------------------------------------
  259.             // Initialize result array to zero.
  260.             // ---------------------------------------------------------------------
  261.             ncalc = nb;
  262.             for (int i = 0; i < nb; ++i) {
  263.                 b[i] = 0;
  264.             }

  266.             // ---------------------------------------------------------------------
  267.             // Branch to use 2-term ascending series for small X and asymptotic
  268.             // form for large X when NB is not too large.
  269.             // ---------------------------------------------------------------------
  270.             double tempa;
  271.             double tempb;
  272.             double tempc;
  273.             double tover;
  274.             if (x < RTNSIG) {
  275.                 // ---------------------------------------------------------------------
  276.                 // Two-term ascending series for small X.
  277.                 // ---------------------------------------------------------------------
  278.                 tempa = 1;
  279.                 alpem = 1 + alpha;
  280.                 double halfx = 0;
  281.                 if (x > ENMTEN) {
  282.                     halfx = 0.5 * x;
  283.                 }
  284.                 if (alpha != 0) {
  285.                     tempa = JdkMath.pow(halfx, alpha) /
  286.                             (alpha * Gamma.value(alpha));
  287.                 }
  288.                 tempb = 0;
  289.                 if (x + 1 > 1) {
  290.                     tempb = -halfx * halfx;
  291.                 }
  292.                 b[0] = tempa + (tempa * tempb / alpem);
  293.                 if (x != 0 && b[0] == 0) {
  294.                     ncalc = 0;
  295.                 }
  296.                 if (nb != 1) {
  297.                     if (x <= 0) {
  298.                         for (int n = 1; n < nb; ++n) {
  299.                             b[n] = 0;
  300.                         }
  301.                     } else {
  302.                         // ---------------------------------------------------------------------
  303.                         // Calculate higher order functions.
  304.                         // ---------------------------------------------------------------------
  305.                         tempc = halfx;
  306.                         tover = tempb != 0 ? ENMTEN / tempb :  2 * ENMTEN / x;
  307.                         for (int n = 1; n < nb; ++n) {
  308.                             tempa /= alpem;
  309.                             alpem += 1;
  310.                             tempa *= tempc;
  311.                             if (tempa <= tover * alpem) {
  312.                                 tempa = 0;
  313.                             }
  314.                             b[n] = tempa + (tempa * tempb / alpem);
  315.                             if (b[n] == 0 && ncalc > n) {
  316.                                 ncalc = n;
  317.                             }
  318.                         }
  319.                     }
  320.                 }
  321.             } else if (x > 25.0 && nb <= magx + 1) {
  322.                 // ---------------------------------------------------------------------
  323.                 // Asymptotic series for X > 25
  324.                 // ---------------------------------------------------------------------
  325.                 final double xc = JdkMath.sqrt(PI2 / x);
  326.                 final double mul = 0.125 / x;
  327.                 final double xin = mul * mul;
  328.                 int m = 0;
  329.                 if (x >= 130.0) {
  330.                     m = 4;
  331.                 } else if (x >= 35.0) {
  332.                     m = 8;
  333.                 } else {
  334.                     m = 11;
  335.                 }

  337.                 final double xm = 4.0 * m;
  338.                 // ---------------------------------------------------------------------
  339.                 // Argument reduction for SIN and COS routines.
  340.                 // ---------------------------------------------------------------------
  341.                 double t = (double) ((int) ((x / TWOPI) + 0.5));
  342.                 final double z = x - t * TOWPI1 - t * TWOPI2 - (alpha + 0.5) / PI2;
  343.                 double vsin = JdkMath.sin(z);
  344.                 double vcos = JdkMath.cos(z);
  345.                 double gnu = 2 * alpha;
  346.                 double capq;
  347.                 double capp;
  348.                 double s;
  349.                 double t1;
  350.                 double xk;
  351.                 for (int i = 1; i <= 2; i++) {
  352.                     s = (xm - 1 - gnu) * (xm - 1 + gnu) * xin * 0.5;
  353.                     t = (gnu - (xm - 3.0)) * (gnu + (xm - 3.0));
  354.                     capp = (s * t) / FACT[2 * m];
  355.                     t1 = (gnu - (xm + 1)) * (gnu + (xm + 1));
  356.                     capq = (s * t1) / FACT[2 * m + 1];
  357.                     xk = xm;
  358.                     int k = 2 * m;
  359.                     t1 = t;

  361.                     for (int j = 2; j <= m; j++) {
  362.                         xk -= 4.0;
  363.                         s = (xk - 1 - gnu) * (xk - 1 + gnu);
  364.                         t = (gnu - (xk - 3.0)) * (gnu + (xk - 3.0));
  365.                         capp = (capp + 1 / FACT[k - 2]) * s * t * xin;
  366.                         capq = (capq + 1 / FACT[k - 1]) * s * t1 * xin;
  367.                         k -= 2;
  368.                         t1 = t;
  369.                     }

  371.                     capp += 1;
  372.                     capq = (capq + 1) * ((gnu * gnu) - 1) * (0.125 / x);
  373.                     b[i - 1] = xc * (capp * vcos - capq * vsin);
  374.                     if (nb == 1) {
  375.                         return new BesselJResult(Arrays.copyOf(b, b.length),
  376.                                                  ncalc);
  377.                     }
  378.                     t = vsin;
  379.                     vsin = -vcos;
  380.                     vcos = t;
  381.                     gnu += 2.0;
  382.                 }

  384.                 // ---------------------------------------------------------------------
  385.                 // If NB > 2, compute J(X,ORDER+I) I = 2, NB-1
  386.                 // ---------------------------------------------------------------------
  387.                 if (nb > 2) {
  388.                     gnu = 2 * alpha + 2.0;
  389.                     for (int j = 2; j < nb; ++j) {
  390.                         b[j] = gnu * b[j - 1] / x - b[j - 2];
  391.                         gnu += 2.0;
  392.                     }
  393.                 }
  394.             } else {
  395.                 // ---------------------------------------------------------------------
  396.                 // Use recurrence to generate results. First initialize the
  397.                 // calculation of P*S.
  398.                 // ---------------------------------------------------------------------
  399.                 final int nbmx = nb - magx;
  400.                 int n = magx + 1;
  401.                 int nstart = 0;
  402.                 int nend = 0;
  403.                 double en = 2 * (n + alpha);
  404.                 double plast = 1;
  405.                 double p = en / x;
  406.                 double pold;
  407.                 // ---------------------------------------------------------------------
  408.                 // Calculate general significance test.
  409.                 // ---------------------------------------------------------------------
  410.                 double test = 2 * ENSIG;
  411.                 boolean readyToInitialize = false;
  412.                 if (nbmx >= 3) {
  413.                     // ---------------------------------------------------------------------
  414.                     // Calculate P*S until N = NB-1. Check for possible
  415.                     // overflow.
  416.                     // ---------------------------------------------------------------------
  417.                     tover = ENTEN / ENSIG;
  418.                     nstart = magx + 2;
  419.                     nend = nb - 1;
  420.                     en = 2 * (nstart - 1 + alpha);
  421.                     double psave;
  422.                     double psavel;
  423.                     for (int k = nstart; k <= nend; k++) {
  424.                         n = k;
  425.                         en += 2.0;
  426.                         pold = plast;
  427.                         plast = p;
  428.                         p = (en * plast / x) - pold;
  429.                         if (p > tover) {
  430.                             // ---------------------------------------------------------------------
  431.                             // To avoid overflow, divide P*S by TOVER. Calculate
  432.                             // P*S until
  433.                             // ABS(P) > 1.
  434.                             // ---------------------------------------------------------------------
  435.                             tover = ENTEN;
  436.                             p /= tover;
  437.                             plast /= tover;
  438.                             psave = p;
  439.                             psavel = plast;
  440.                             nstart = n + 1;
  441.                             do {
  442.                                 n += 1;
  443.                                 en += 2.0;
  444.                                 pold = plast;
  445.                                 plast = p;
  446.                                 p = (en * plast / x) - pold;
  447.                             } while (p <= 1);
  448.                             tempb = en / x;
  449.                             // ---------------------------------------------------------------------
  450.                             // Calculate backward test and find NCALC, the
  451.                             // highest N such that
  452.                             // the test is passed.
  453.                             // ---------------------------------------------------------------------
  454.                             test = pold * plast * (0.5 - 0.5 / (tempb * tempb));
  455.                             test /= ENSIG;
  456.                             p = plast * tover;
  457.                             n -= 1;
  458.                             en -= 2.0;
  459.                             nend = JdkMath.min(nb, n);
  460.                             for (int l = nstart; l <= nend; l++) {
  461.                                 pold = psavel;
  462.                                 psavel = psave;
  463.                                 psave = (en * psavel / x) - pold;
  464.                                 if (psave * psavel > test) {
  465.                                     ncalc = l - 1;
  466.                                     readyToInitialize = true;
  467.                                     break;
  468.                                 }
  469.                             }
  470.                             ncalc = nend;
  471.                             readyToInitialize = true;
  472.                             break;
  473.                         }
  474.                     }
  475.                     if (!readyToInitialize) {
  476.                         n = nend;
  477.                         en = 2 * (n + alpha);
  478.                         // ---------------------------------------------------------------------
  479.                         // Calculate special significance test for NBMX > 2.
  480.                         // ---------------------------------------------------------------------
  481.                         test = JdkMath.max(test, JdkMath.sqrt(plast * ENSIG) *
  482.                                                   JdkMath.sqrt(2 * p));
  483.                     }
  484.                 }
  485.                 // ---------------------------------------------------------------------
  486.                 // Calculate P*S until significance test passes.
  487.                 // ---------------------------------------------------------------------
  488.                 if (!readyToInitialize) {
  489.                     do {
  490.                         n += 1;
  491.                         en += 2.0;
  492.                         pold = plast;
  493.                         plast = p;
  494.                         p = (en * plast / x) - pold;
  495.                     } while (p < test);
  496.                 }
  497.                 // ---------------------------------------------------------------------
  498.                 // Initialize the backward recursion and the normalization sum.
  499.                 // ---------------------------------------------------------------------
  500.                 n += 1;
  501.                 en += 2.0;
  502.                 tempb = 0;
  503.                 tempa = 1 / p;
  504.                 int m = (2 * n) - 4 * (n / 2);
  505.                 double sum = 0;
  506.                 double em = (double) (n / 2);
  507.                 alpem = em - 1 + alpha;
  508.                 alp2em = 2 * em + alpha;
  509.                 if (m != 0) {
  510.                     sum = tempa * alpem * alp2em / em;
  511.                 }
  512.                 nend = n - nb;

  514.                 boolean readyToNormalize = false;
  515.                 boolean calculatedB0 = false;

  517.                 // ---------------------------------------------------------------------
  518.                 // Recur backward via difference equation, calculating (but not
  519.                 // storing) B(N), until N = NB.
  520.                 // ---------------------------------------------------------------------
  521.                 for (int l = 1; l <= nend; l++) {
  522.                     n -= 1;
  523.                     en -= 2.0;
  524.                     tempc = tempb;
  525.                     tempb = tempa;
  526.                     tempa = (en * tempb / x) - tempc;
  527.                     m = 2 - m;
  528.                     if (m != 0) {
  529.                         em -= 1;
  530.                         alp2em = 2 * em + alpha;
  531.                         if (n == 1) {
  532.                             break;
  533.                         }
  534.                         alpem = em - 1 + alpha;
  535.                         if (alpem == 0) {
  536.                             alpem = 1;
  537.                         }
  538.                         sum = (sum + tempa * alp2em) * alpem / em;
  539.                     }
  540.                 }

  542.                 // ---------------------------------------------------------------------
  543.                 // Store B(NB).
  544.                 // ---------------------------------------------------------------------
  545.                 b[n - 1] = tempa;
  546.                 if (nend >= 0) {
  547.                     if (nb <= 1) {
  548.                         alp2em = alpha;
  549.                         if (alpha + 1 == 1) {
  550.                             alp2em = 1;
  551.                         }
  552.                         sum += b[0] * alp2em;
  553.                         readyToNormalize = true;
  554.                     } else {
  555.                         // ---------------------------------------------------------------------
  556.                         // Calculate and store B(NB-1).
  557.                         // ---------------------------------------------------------------------
  558.                         n -= 1;
  559.                         en -= 2.0;
  560.                         b[n - 1] = (en * tempa / x) - tempb;
  561.                         if (n == 1) {
  562.                             calculatedB0 = true;
  563.                         } else {
  564.                             m = 2 - m;
  565.                             if (m != 0) {
  566.                                 em -= 1;
  567.                                 alp2em = 2 * em + alpha;
  568.                                 alpem = em - 1 + alpha;
  569.                                 if (alpem == 0) {
  570.                                     alpem = 1;
  571.                                 }

  573.                                 sum = (sum + (b[n - 1] * alp2em)) * alpem / em;
  574.                             }
  575.                         }
  576.                     }
  577.                 }
  578.                 if (!readyToNormalize && !calculatedB0) {
  579.                     nend = n - 2;
  580.                     if (nend != 0) {
  581.                         // ---------------------------------------------------------------------
  582.                         // Calculate via difference equation and store B(N),
  583.                         // until N = 2.
  584.                         // ---------------------------------------------------------------------

  586.                         for (int l = 1; l <= nend; l++) {
  587.                             n -= 1;
  588.                             en -= 2.0;
  589.                             b[n - 1] = (en * b[n] / x) - b[n + 1];
  590.                             m = 2 - m;
  591.                             if (m != 0) {
  592.                                 em -= 1;
  593.                                 alp2em = 2 * em + alpha;
  594.                                 alpem = em - 1 + alpha;
  595.                                 if (alpem == 0) {
  596.                                     alpem = 1;
  597.                                 }

  599.                                 sum = (sum + b[n - 1] * alp2em) * alpem / em;
  600.                             }
  601.                         }
  602.                     }
  603.                 }
  604.                 // ---------------------------------------------------------------------
  605.                 // Calculate b[0]
  606.                 // ---------------------------------------------------------------------
  607.                 if (!readyToNormalize) {
  608.                     if (!calculatedB0) {
  609.                         b[0] = 2.0 * (alpha + 1) * b[1] / x - b[2];
  610.                     }
  611.                     em -= 1;
  612.                     alp2em = 2 * em + alpha;
  613.                     if (alp2em == 0) {
  614.                         alp2em = 1;
  615.                     }
  616.                     sum += b[0] * alp2em;
  617.                 }
  618.                 // ---------------------------------------------------------------------
  619.                 // Normalize. Divide all B(N) by sum.
  620.                 // ---------------------------------------------------------------------

  622.                 if (JdkMath.abs(alpha) > 1e-16) {
  623.                     sum *= Gamma.value(alpha) * JdkMath.pow(x * 0.5, -alpha);
  624.                 }
  625.                 tempa = ENMTEN;
  626.                 if (sum > 1) {
  627.                     tempa *= sum;
  628.                 }

  630.                 for (n = 0; n < nb; n++) {
  631.                     if (JdkMath.abs(b[n]) < tempa) {
  632.                         b[n] = 0;
  633.                     }
  634.                     b[n] /= sum;
  635.                 }
  636.             }
  637.             // ---------------------------------------------------------------------
  638.             // Error return -- X, NB, or ALPHA is out of range.
  639.             // ---------------------------------------------------------------------
  640.         } else {
  641.             if (b.length > 0) {
  642.                 b[0] = 0;
  643.             }
  644.             ncalc = JdkMath.min(nb, 0) - 1;
  645.         }
  646.         return new BesselJResult(Arrays.copyOf(b, b.length), ncalc);
  647.     }
  648. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT