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

package org.apache.commons.math3.special;

import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.exception.ConvergenceException;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.MathArrays;

This class provides computation methods related to Bessel
functions of the first kind. Detailed descriptions of these functions are
available in Wikipedia, Abrabowitz and
Stegun (Ch. 9-11), and DLMF (Ch. 10).

This implementation is based on the rjbesl Fortran routine at
Netlib.

From the Fortran code: 

This program is based on a program written by David J. Sookne (2) that
computes values of the Bessel functions J or I of real argument and integer
order. Modifications include the restriction of the computation to the J
Bessel function of non-negative real argument, the extension of the
computation to arbitrary positive order, and the elimination of most
underflow.

References:

"A Note on Backward Recurrence Algorithms," Olver, F. W. J., and Sookne,
D. J., Math. Comp. 26, 1972, pp 941-947.
"Bessel Functions of Real Argument and Integer Order," Sookne, D. J., NBS
Jour. of Res. B. 77B, 1973, pp 125-132.
 

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

    // ---------------------------------------------------------------------
    // Mathematical constants
    // ---------------------------------------------------------------------

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

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

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

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

    // ---------------------------------------------------------------------
    // Machine-dependent parameters
    // ---------------------------------------------------------------------

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

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

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

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

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

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

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

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

    
Create a new BesselJ with the given order.
Params:
order – order of the function computed when using value(double)./**
     * Create a new BesselJ with the given order.
     *
     * @param order order of the function computed when using {@link #value(double)}.
     */
    public BesselJ(double order) {
        this.order = order;
    }

    
Returns the value of the constructed Bessel function of the first kind,
for the passed argument.
Params:
x – Argument
Throws:
MathIllegalArgumentException – if x is too large relative to order
ConvergenceException – if the algorithm fails to converge
Returns:
Value of the Bessel function at x/**
     * Returns the value of the constructed Bessel function of the first kind,
     * for the passed argument.
     *
     * @param x Argument
     * @return Value of the Bessel function at x
     * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order}
     * @throws ConvergenceException if the algorithm fails to converge
     */
    public double value(double x)
        throws MathIllegalArgumentException, ConvergenceException {
        return BesselJ.value(order, x);
    }

    
Returns the first Bessel function, \(J_{order}(x)\).
Params:
order – Order of the Bessel function
x – Argument
Throws:
MathIllegalArgumentException – if x is too large relative to order
ConvergenceException – if the algorithm fails to converge
Returns:
Value of the Bessel function of the first kind, \(J_{order}(x)\)/**
     * Returns the first Bessel function, \(J_{order}(x)\).
     *
     * @param order Order of the Bessel function
     * @param x Argument
     * @return Value of the Bessel function of the first kind, \(J_{order}(x)\)
     * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order}
     * @throws ConvergenceException if the algorithm fails to converge
     */
    public static double value(double order, double x)
        throws MathIllegalArgumentException, ConvergenceException {
        final int n = (int) order;
        final double alpha = order - n;
        final int nb = n + 1;
        final BesselJResult res = rjBesl(x, alpha, nb);

        if (res.nVals >= nb) {
            return res.vals[n];
        } else if (res.nVals < 0) {
            throw new MathIllegalArgumentException(LocalizedFormats.BESSEL_FUNCTION_BAD_ARGUMENT,order, x);
        } else if (FastMath.abs(res.vals[res.nVals - 1]) < 1e-100) {
            return res.vals[n]; // underflow; return value (will be zero)
        }
        throw new ConvergenceException(LocalizedFormats.BESSEL_FUNCTION_FAILED_CONVERGENCE, order, x);
    }

    
Encapsulates the results returned by BesselJ.rjBesl(double, double, int). 
 getVals() returns the computed function values. getnVals() is the number of values among those returned by getnVals() that can be considered accurate. 


nVals < 0: An argument is out of range. For example, nb <= 0, alpha
< 0 or > 1, or x is too large. In this case, b(0) is set to zero, the
remainder of the b-vector is not calculated, and nVals is set to
MIN(nb,0) - 1 so that nVals != nb.
nb > nVals > 0: Not all requested function values could be calculated
accurately. This usually occurs because nb is much larger than abs(x). In
this case, b(n) is calculated to the desired accuracy for n < nVals, but
precision is lost for nVals < n <= nb. If b(n) does not vanish for n >
nVals (because it is too small to be represented), and b(n)/b(nVals) =
\(10^{-k}\), then only the first NSIG-k significant figures of b(n) can be
trusted.

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

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

        
Valid value count /** Valid value count */
        private final int nVals;

        
Create a new BesselJResult with the given values and valid value count.
Params:
b – values
n – count of valid values/**
         * Create a new BesselJResult with the given values and valid value count.
         *
         * @param b values
         * @param n count of valid values
         */
        public BesselJResult(double[] b, int n) {
            vals = MathArrays.copyOf(b, b.length);
            nVals = n;
        }

        
Returns:
the computed function values/**
         * @return the computed function values
         */
        public double[] getVals() {
            return MathArrays.copyOf(vals, vals.length);
        }

        
Returns:
the number of valid function values (normally the same as the length of the array returned by getnVals())/**
         * @return the number of valid function values (normally the same as the
         *         length of the array returned by {@link #getnVals()})
         */
        public int getnVals() {
            return nVals;
        }
    }

    
Calculates Bessel functions \(J_{n+alpha}(x)\) for
non-negative argument x, and non-negative order n + alpha.

Before using the output vector, the user should check that
nVals = nb, i.e., all orders have been calculated to the desired accuracy.
See BesselResult class javadoc for details on return values.

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

        int ncalc = 0;
        double alpem = 0;
        double alp2em = 0;

        // ---------------------------------------------------------------------
        // Check for out of range arguments.
        // ---------------------------------------------------------------------
        final int magx = (int) x;
        if ((nb > 0) && (x >= X_MIN) && (x <= X_MAX) && (alpha >= 0) &&
            (alpha < 1)) {
            // ---------------------------------------------------------------------
            // Initialize result array to zero.
            // ---------------------------------------------------------------------
            ncalc = nb;
            for (int i = 0; i < nb; ++i) {
                b[i] = 0;
            }

            // ---------------------------------------------------------------------
            // Branch to use 2-term ascending series for small X and asymptotic
            // form for large X when NB is not too large.
            // ---------------------------------------------------------------------
            double tempa;
            double tempb;
            double tempc;
            double tover;
            if (x < RTNSIG) {
                // ---------------------------------------------------------------------
                // Two-term ascending series for small X.
                // ---------------------------------------------------------------------
                tempa = 1;
                alpem = 1 + alpha;
                double halfx = 0;
                if (x > ENMTEN) {
                    halfx = 0.5 * x;
                }
                if (alpha != 0) {
                    tempa = FastMath.pow(halfx, alpha) /
                            (alpha * Gamma.gamma(alpha));
                }
                tempb = 0;
                if (x + 1 > 1) {
                    tempb = -halfx * halfx;
                }
                b[0] = tempa + (tempa * tempb / alpem);
                if ((x != 0) && (b[0] == 0)) {
                    ncalc = 0;
                }
                if (nb != 1) {
                    if (x <= 0) {
                        for (int n = 1; n < nb; ++n) {
                            b[n] = 0;
                        }
                    } else {
                        // ---------------------------------------------------------------------
                        // Calculate higher order functions.
                        // ---------------------------------------------------------------------
                        tempc = halfx;
                        tover = tempb != 0 ? ENMTEN / tempb :  2 * ENMTEN / x;
                        for (int n = 1; n < nb; ++n) {
                            tempa /= alpem;
                            alpem += 1;
                            tempa *= tempc;
                            if (tempa <= tover * alpem) {
                                tempa = 0;
                            }
                            b[n] = tempa + (tempa * tempb / alpem);
                            if ((b[n] == 0) && (ncalc > n)) {
                                ncalc = n;
                            }
                        }
                    }
                }
            } else if ((x > 25.0) && (nb <= magx + 1)) {
                // ---------------------------------------------------------------------
                // Asymptotic series for X > 25
                // ---------------------------------------------------------------------
                final double xc = FastMath.sqrt(PI2 / x);
                final double mul = 0.125 / x;
                final double xin = mul * mul;
                int m = 0;
                if (x >= 130.0) {
                    m = 4;
                } else if (x >= 35.0) {
                    m = 8;
                } else {
                    m = 11;
                }

                final double xm = 4.0 * m;
                // ---------------------------------------------------------------------
                // Argument reduction for SIN and COS routines.
                // ---------------------------------------------------------------------
                double t = (double) ((int) ((x / TWOPI) + 0.5));
                final double z = x - t * TOWPI1 - t * TWOPI2 - (alpha + 0.5) / PI2;
                double vsin = FastMath.sin(z);
                double vcos = FastMath.cos(z);
                double gnu = 2 * alpha;
                double capq;
                double capp;
                double s;
                double t1;
                double xk;
                for (int i = 1; i <= 2; i++) {
                    s = (xm - 1 - gnu) * (xm - 1 + gnu) * xin * 0.5;
                    t = (gnu - (xm - 3.0)) * (gnu + (xm - 3.0));
                    capp = (s * t) / FACT[2 * m];
                    t1 = (gnu - (xm + 1)) * (gnu + (xm + 1));
                    capq = (s * t1) / FACT[2 * m + 1];
                    xk = xm;
                    int k = 2 * m;
                    t1 = t;

                    for (int j = 2; j <= m; j++) {
                        xk -= 4.0;
                        s = (xk - 1 - gnu) * (xk - 1 + gnu);
                        t = (gnu - (xk - 3.0)) * (gnu + (xk - 3.0));
                        capp = (capp + 1 / FACT[k - 2]) * s * t * xin;
                        capq = (capq + 1 / FACT[k - 1]) * s * t1 * xin;
                        k -= 2;
                        t1 = t;
                    }

                    capp += 1;
                    capq = (capq + 1) * ((gnu * gnu) - 1) * (0.125 / x);
                    b[i - 1] = xc * (capp * vcos - capq * vsin);
                    if (nb == 1) {
                        return new BesselJResult(MathArrays.copyOf(b, b.length),
                                                 ncalc);
                    }
                    t = vsin;
                    vsin = -vcos;
                    vcos = t;
                    gnu += 2.0;
                }

                // ---------------------------------------------------------------------
                // If NB > 2, compute J(X,ORDER+I) I = 2, NB-1
                // ---------------------------------------------------------------------
                if (nb > 2) {
                    gnu = 2 * alpha + 2.0;
                    for (int j = 2; j < nb; ++j) {
                        b[j] = gnu * b[j - 1] / x - b[j - 2];
                        gnu += 2.0;
                    }
                }
            } else {
                // ---------------------------------------------------------------------
                // Use recurrence to generate results. First initialize the
                // calculation of P*S.
                // ---------------------------------------------------------------------
                final int nbmx = nb - magx;
                int n = magx + 1;
                int nstart = 0;
                int nend = 0;
                double en = 2 * (n + alpha);
                double plast = 1;
                double p = en / x;
                double pold;
                // ---------------------------------------------------------------------
                // Calculate general significance test.
                // ---------------------------------------------------------------------
                double test = 2 * ENSIG;
                boolean readyToInitialize = false;
                if (nbmx >= 3) {
                    // ---------------------------------------------------------------------
                    // Calculate P*S until N = NB-1. Check for possible
                    // overflow.
                    // ---------------------------------------------------------------------
                    tover = ENTEN / ENSIG;
                    nstart = magx + 2;
                    nend = nb - 1;
                    en = 2 * (nstart - 1 + alpha);
                    double psave;
                    double psavel;
                    for (int k = nstart; k <= nend; k++) {
                        n = k;
                        en += 2.0;
                        pold = plast;
                        plast = p;
                        p = (en * plast / x) - pold;
                        if (p > tover) {
                            // ---------------------------------------------------------------------
                            // To avoid overflow, divide P*S by TOVER. Calculate
                            // P*S until
                            // ABS(P) > 1.
                            // ---------------------------------------------------------------------
                            tover = ENTEN;
                            p /= tover;
                            plast /= tover;
                            psave = p;
                            psavel = plast;
                            nstart = n + 1;
                            do {
                                n += 1;
                                en += 2.0;
                                pold = plast;
                                plast = p;
                                p = (en * plast / x) - pold;
                            } while (p <= 1);
                            tempb = en / x;
                            // ---------------------------------------------------------------------
                            // Calculate backward test and find NCALC, the
                            // highest N such that
                            // the test is passed.
                            // ---------------------------------------------------------------------
                            test = pold * plast * (0.5 - 0.5 / (tempb * tempb));
                            test /= ENSIG;
                            p = plast * tover;
                            n -= 1;
                            en -= 2.0;
                            nend = FastMath.min(nb, n);
                            for (int l = nstart; l <= nend; l++) {
                                pold = psavel;
                                psavel = psave;
                                psave = (en * psavel / x) - pold;
                                if (psave * psavel > test) {
                                    ncalc = l - 1;
                                    readyToInitialize = true;
                                    break;
                                }
                            }
                            ncalc = nend;
                            readyToInitialize = true;
                            break;
                        }
                    }
                    if (!readyToInitialize) {
                        n = nend;
                        en = 2 * (n + alpha);
                        // ---------------------------------------------------------------------
                        // Calculate special significance test for NBMX > 2.
                        // ---------------------------------------------------------------------
                        test = FastMath.max(test, FastMath.sqrt(plast * ENSIG) *
                                                  FastMath.sqrt(2 * p));
                    }
                }
                // ---------------------------------------------------------------------
                // Calculate P*S until significance test passes.
                // ---------------------------------------------------------------------
                if (!readyToInitialize) {
                    do {
                        n += 1;
                        en += 2.0;
                        pold = plast;
                        plast = p;
                        p = (en * plast / x) - pold;
                    } while (p < test);
                }
                // ---------------------------------------------------------------------
                // Initialize the backward recursion and the normalization sum.
                // ---------------------------------------------------------------------
                n += 1;
                en += 2.0;
                tempb = 0;
                tempa = 1 / p;
                int m = (2 * n) - 4 * (n / 2);
                double sum = 0;
                double em = (double) (n / 2);
                alpem = em - 1 + alpha;
                alp2em = 2 * em + alpha;
                if (m != 0) {
                    sum = tempa * alpem * alp2em / em;
                }
                nend = n - nb;

                boolean readyToNormalize = false;
                boolean calculatedB0 = false;

                // ---------------------------------------------------------------------
                // Recur backward via difference equation, calculating (but not
                // storing) B(N), until N = NB.
                // ---------------------------------------------------------------------
                for (int l = 1; l <= nend; l++) {
                    n -= 1;
                    en -= 2.0;
                    tempc = tempb;
                    tempb = tempa;
                    tempa = (en * tempb / x) - tempc;
                    m = 2 - m;
                    if (m != 0) {
                        em -= 1;
                        alp2em = 2 * em + alpha;
                        if (n == 1) {
                            break;
                        }
                        alpem = em - 1 + alpha;
                        if (alpem == 0) {
                            alpem = 1;
                        }
                        sum = (sum + tempa * alp2em) * alpem / em;
                    }
                }

                // ---------------------------------------------------------------------
                // Store B(NB).
                // ---------------------------------------------------------------------
                b[n - 1] = tempa;
                if (nend >= 0) {
                    if (nb <= 1) {
                        alp2em = alpha;
                        if (alpha + 1 == 1) {
                            alp2em = 1;
                        }
                        sum += b[0] * alp2em;
                        readyToNormalize = true;
                    } else {
                        // ---------------------------------------------------------------------
                        // Calculate and store B(NB-1).
                        // ---------------------------------------------------------------------
                        n -= 1;
                        en -= 2.0;
                        b[n - 1] = (en * tempa / x) - tempb;
                        if (n == 1) {
                            calculatedB0 = true;
                        } else {
                            m = 2 - m;
                            if (m != 0) {
                                em -= 1;
                                alp2em = 2 * em + alpha;
                                alpem = em - 1 + alpha;
                                if (alpem == 0) {
                                    alpem = 1;
                                }

                                sum = (sum + (b[n - 1] * alp2em)) * alpem / em;
                            }
                        }
                    }
                }
                if (!readyToNormalize && !calculatedB0) {
                    nend = n - 2;
                    if (nend != 0) {
                        // ---------------------------------------------------------------------
                        // Calculate via difference equation and store B(N),
                        // until N = 2.
                        // ---------------------------------------------------------------------

                        for (int l = 1; l <= nend; l++) {
                            n -= 1;
                            en -= 2.0;
                            b[n - 1] = (en * b[n] / x) - b[n + 1];
                            m = 2 - m;
                            if (m != 0) {
                                em -= 1;
                                alp2em = 2 * em + alpha;
                                alpem = em - 1 + alpha;
                                if (alpem == 0) {
                                    alpem = 1;
                                }

                                sum = (sum + b[n - 1] * alp2em) * alpem / em;
                            }
                        }
                    }
                }
                // ---------------------------------------------------------------------
                // Calculate b[0]
                // ---------------------------------------------------------------------
                if (!readyToNormalize) {
                    if (!calculatedB0) {
                        b[0] = 2.0 * (alpha + 1) * b[1] / x - b[2];
                    }
                    em -= 1;
                    alp2em = 2 * em + alpha;
                    if (alp2em == 0) {
                        alp2em = 1;
                    }
                    sum += b[0] * alp2em;
                }
                // ---------------------------------------------------------------------
                // Normalize. Divide all B(N) by sum.
                // ---------------------------------------------------------------------

                if (FastMath.abs(alpha) > 1e-16) {
                    sum *= Gamma.gamma(alpha) * FastMath.pow(x * 0.5, -alpha);
                }
                tempa = ENMTEN;
                if (sum > 1) {
                    tempa *= sum;
                }

                for (n = 0; n < nb; n++) {
                    if (FastMath.abs(b[n]) < tempa) {
                        b[n] = 0;
                    }
                    b[n] /= sum;
                }
            }
            // ---------------------------------------------------------------------
            // Error return -- X, NB, or ALPHA is out of range.
            // ---------------------------------------------------------------------
        } else {
            if (b.length > 0) {
                b[0] = 0;
            }
            ncalc = FastMath.min(nb, 0) - 1;
        }
        return new BesselJResult(MathArrays.copyOf(b, b.length), ncalc);
    }
}