/*
 * 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.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.OutOfRangeException;
import org.apache.commons.math3.util.ContinuedFraction;
import org.apache.commons.math3.util.FastMath;


This is a utility class that provides computation methods related to the
Beta family of functions.

 Implementation of logBeta(double, double) is based on the algorithms described in 

Didonato and Morris
    (1986), Computation of the Incomplete Gamma Function Ratios
    and their Inverse, TOMS 12(4), 377-393,
Didonato and Morris
    (1992), Algorithm 708: Significant Digit Computation of the
    Incomplete Beta Function Ratios, TOMS 18(3), 360-373,

and implemented in the
NSWC Library of Mathematical Functions,
available
here.
This library is "approved for public release", and the
Copyright guidance
indicates that unless otherwise stated in the code, all FORTRAN functions in
this library are license free. Since no such notice appears in the code these
functions can safely be ported to Commons-Math.

/**
 * <p>
 * This is a utility class that provides computation methods related to the
 * Beta family of functions.
 * </p>
 * <p>
 * Implementation of {@link #logBeta(double, double)} is based on the
 * algorithms described in
 * <ul>
 * <li><a href="http://dx.doi.org/10.1145/22721.23109">Didonato and Morris
 *     (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios
 *     and their Inverse</em>, TOMS 12(4), 377-393,</li>
 * <li><a href="http://dx.doi.org/10.1145/131766.131776">Didonato and Morris
 *     (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the
 *     Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li>
 * </ul>
 * and implemented in the
 * <a href="http://www.dtic.mil/docs/citations/ADA476840">NSWC Library of Mathematical Functions</a>,
 * available
 * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html">here</a>.
 * This library is "approved for public release", and the
 * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf">Copyright guidance</a>
 * indicates that unless otherwise stated in the code, all FORTRAN functions in
 * this library are license free. Since no such notice appears in the code these
 * functions can safely be ported to Commons-Math.
 * </p>
 *
 *
 */
public class Beta {
    
Maximum allowed numerical error. /** Maximum allowed numerical error. */
    private static final double DEFAULT_EPSILON = 1E-14;

    
The constant value of ½log 2π. /** The constant value of ½log 2π. */
    private static final double HALF_LOG_TWO_PI = .9189385332046727;

    

The coefficients of the series expansion of the Δ function. This function
is defined as follows

Δ(x) = log Γ(x) - (x - 0.5) log a + a - 0.5 log 2π,

see equation (23) in Didonato and Morris (1992). The series expansion,
which applies for x ≥ 10, reads

                14
               ====
            1  \                2 n
    Δ(x) = ---  >    d  (10 / x)
            x  /      n
               ====
               n = 0
/**
     * <p>
     * The coefficients of the series expansion of the Δ function. This function
     * is defined as follows
     * </p>
     * <center>Δ(x) = log Γ(x) - (x - 0.5) log a + a - 0.5 log 2π,</center>
     * <p>
     * see equation (23) in Didonato and Morris (1992). The series expansion,
     * which applies for x ≥ 10, reads
     * </p>
     * <pre>
     *                 14
     *                ====
     *             1  \                2 n
     *     Δ(x) = ---  >    d  (10 / x)
     *             x  /      n
     *                ====
     *                n = 0
     * <pre>
     */
    private static final double[] DELTA = {
        .833333333333333333333333333333E-01,
        -.277777777777777777777777752282E-04,
        .793650793650793650791732130419E-07,
        -.595238095238095232389839236182E-09,
        .841750841750832853294451671990E-11,
        -.191752691751854612334149171243E-12,
        .641025640510325475730918472625E-14,
        -.295506514125338232839867823991E-15,
        .179643716359402238723287696452E-16,
        -.139228964661627791231203060395E-17,
        .133802855014020915603275339093E-18,
        -.154246009867966094273710216533E-19,
        .197701992980957427278370133333E-20,
        -.234065664793997056856992426667E-21,
        .171348014966398575409015466667E-22
    };

    
Default constructor.  Prohibit instantiation.
/**
     * Default constructor.  Prohibit instantiation.
     */
    private Beta() {}

    
Returns the

regularized beta function I(x, a, b).
Params:
x – Value.
a – Parameter a.
b – Parameter b.
Throws:
MaxCountExceededException – 
if the algorithm fails to converge.
Returns:
the regularized beta function I(x, a, b)./**
     * Returns the
     * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
     * regularized beta function</a> I(x, a, b).
     *
     * @param x Value.
     * @param a Parameter {@code a}.
     * @param b Parameter {@code b}.
     * @return the regularized beta function I(x, a, b).
     * @throws org.apache.commons.math3.exception.MaxCountExceededException
     * if the algorithm fails to converge.
     */
    public static double regularizedBeta(double x, double a, double b) {
        return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
    }

    
Returns the

regularized beta function I(x, a, b).
Params:
x – Value.
a – Parameter a.
b – Parameter b.
epsilon – When the absolute value of the nth item in the
series is less than epsilon the approximation ceases to calculate
further elements in the series.
Throws:
MaxCountExceededException – 
if the algorithm fails to converge.
Returns:
the regularized beta function I(x, a, b)/**
     * Returns the
     * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
     * regularized beta function</a> I(x, a, b).
     *
     * @param x Value.
     * @param a Parameter {@code a}.
     * @param b Parameter {@code b}.
     * @param epsilon When the absolute value of the nth item in the
     * series is less than epsilon the approximation ceases to calculate
     * further elements in the series.
     * @return the regularized beta function I(x, a, b)
     * @throws org.apache.commons.math3.exception.MaxCountExceededException
     * if the algorithm fails to converge.
     */
    public static double regularizedBeta(double x,
                                         double a, double b,
                                         double epsilon) {
        return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE);
    }

    
Returns the regularized beta function I(x, a, b).
Params:
x – the value.
a – Parameter a.
b – Parameter b.
maxIterations – Maximum number of "iterations" to complete.
Throws:
MaxCountExceededException – 
if the algorithm fails to converge.
Returns:
the regularized beta function I(x, a, b)/**
     * Returns the regularized beta function I(x, a, b).
     *
     * @param x the value.
     * @param a Parameter {@code a}.
     * @param b Parameter {@code b}.
     * @param maxIterations Maximum number of "iterations" to complete.
     * @return the regularized beta function I(x, a, b)
     * @throws org.apache.commons.math3.exception.MaxCountExceededException
     * if the algorithm fails to converge.
     */
    public static double regularizedBeta(double x,
                                         double a, double b,
                                         int maxIterations) {
        return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations);
    }

    
Returns the regularized beta function I(x, a, b).
The implementation of this method is based on:



Regularized Beta Function.


Regularized Beta Function.

Params:
x – the value.
a – Parameter a.
b – Parameter b.
epsilon – When the absolute value of the nth item in the
series is less than epsilon the approximation ceases to calculate
further elements in the series.
maxIterations – Maximum number of "iterations" to complete.
Throws:
MaxCountExceededException – 
if the algorithm fails to converge.
Returns:
the regularized beta function I(x, a, b)/**
     * Returns the regularized beta function I(x, a, b).
     *
     * The implementation of this method is based on:
     * <ul>
     * <li>
     * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
     * Regularized Beta Function</a>.</li>
     * <li>
     * <a href="http://functions.wolfram.com/06.21.10.0001.01">
     * Regularized Beta Function</a>.</li>
     * </ul>
     *
     * @param x the value.
     * @param a Parameter {@code a}.
     * @param b Parameter {@code b}.
     * @param epsilon When the absolute value of the nth item in the
     * series is less than epsilon the approximation ceases to calculate
     * further elements in the series.
     * @param maxIterations Maximum number of "iterations" to complete.
     * @return the regularized beta function I(x, a, b)
     * @throws org.apache.commons.math3.exception.MaxCountExceededException
     * if the algorithm fails to converge.
     */
    public static double regularizedBeta(double x,
                                         final double a, final double b,
                                         double epsilon, int maxIterations) {
        double ret;

        if (Double.isNaN(x) ||
            Double.isNaN(a) ||
            Double.isNaN(b) ||
            x < 0 ||
            x > 1 ||
            a <= 0 ||
            b <= 0) {
            ret = Double.NaN;
        } else if (x > (a + 1) / (2 + b + a) &&
                   1 - x <= (b + 1) / (2 + b + a)) {
            ret = 1 - regularizedBeta(1 - x, b, a, epsilon, maxIterations);
        } else {
            ContinuedFraction fraction = new ContinuedFraction() {

                
{@inheritDoc} /** {@inheritDoc} */
                @Override
                protected double getB(int n, double x) {
                    double ret;
                    double m;
                    if (n % 2 == 0) { // even
                        m = n / 2.0;
                        ret = (m * (b - m) * x) /
                            ((a + (2 * m) - 1) * (a + (2 * m)));
                    } else {
                        m = (n - 1.0) / 2.0;
                        ret = -((a + m) * (a + b + m) * x) /
                                ((a + (2 * m)) * (a + (2 * m) + 1.0));
                    }
                    return ret;
                }

                
{@inheritDoc} /** {@inheritDoc} */
                @Override
                protected double getA(int n, double x) {
                    return 1.0;
                }
            };
            ret = FastMath.exp((a * FastMath.log(x)) + (b * FastMath.log1p(-x)) -
                FastMath.log(a) - logBeta(a, b)) *
                1.0 / fraction.evaluate(x, epsilon, maxIterations);
        }

        return ret;
    }

    
Returns the natural logarithm of the beta function B(a, b).
The implementation of this method is based on:


Beta Function, equation (1).

Params:
a – Parameter a.
b – Parameter b.
epsilon – This parameter is ignored.
maxIterations – This parameter is ignored.
Returns:
log(B(a, b)).
Deprecated:
as of version 3.1, this method is deprecated as the computation of the beta function is no longer iterative; it will be removed in version 4.0. Current implementation of this method internally calls logBeta(double, double)./**
     * Returns the natural logarithm of the beta function B(a, b).
     *
     * The implementation of this method is based on:
     * <ul>
     * <li><a href="http://mathworld.wolfram.com/BetaFunction.html">
     * Beta Function</a>, equation (1).</li>
     * </ul>
     *
     * @param a Parameter {@code a}.
     * @param b Parameter {@code b}.
     * @param epsilon This parameter is ignored.
     * @param maxIterations This parameter is ignored.
     * @return log(B(a, b)).
     * @deprecated as of version 3.1, this method is deprecated as the
     * computation of the beta function is no longer iterative; it will be
     * removed in version 4.0. Current implementation of this method
     * internally calls {@link #logBeta(double, double)}.
     */
    @Deprecated
    public static double logBeta(double a, double b,
                                 double epsilon,
                                 int maxIterations) {

        return logBeta(a, b);
    }


    
Returns the value of log Γ(a + b) for 1 ≤ a, b ≤ 2. Based on the
NSWC Library of Mathematics Subroutines double precision implementation, DGSMLN. In BetaTest.testLogGammaSum(), this private method is accessed through reflection. 
Params:
a – First argument.
b – Second argument.
Throws:
OutOfRangeException – if a or b is lower than 1.0 or greater than 2.0.
Returns:
the value of log(Gamma(a + b))./**
     * Returns the value of log Γ(a + b) for 1 ≤ a, b ≤ 2. Based on the
     * <em>NSWC Library of Mathematics Subroutines</em> double precision
     * implementation, {@code DGSMLN}. In {@code BetaTest.testLogGammaSum()},
     * this private method is accessed through reflection.
     *
     * @param a First argument.
     * @param b Second argument.
     * @return the value of {@code log(Gamma(a + b))}.
     * @throws OutOfRangeException if {@code a} or {@code b} is lower than
     * {@code 1.0} or greater than {@code 2.0}.
     */
    private static double logGammaSum(final double a, final double b)
        throws OutOfRangeException {

        if ((a < 1.0) || (a > 2.0)) {
            throw new OutOfRangeException(a, 1.0, 2.0);
        }
        if ((b < 1.0) || (b > 2.0)) {
            throw new OutOfRangeException(b, 1.0, 2.0);
        }

        final double x = (a - 1.0) + (b - 1.0);
        if (x <= 0.5) {
            return Gamma.logGamma1p(1.0 + x);
        } else if (x <= 1.5) {
            return Gamma.logGamma1p(x) + FastMath.log1p(x);
        } else {
            return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
        }
    }

    
Returns the value of log[Γ(b) / Γ(a + b)] for a ≥ 0 and b ≥ 10. Based on
the NSWC Library of Mathematics Subroutines double precision implementation, DLGDIV. In BetaTest.testLogGammaMinusLogGammaSum(), this private method is accessed through reflection. 
Params:
a – First argument.
b – Second argument.
Throws:
NumberIsTooSmallException – if a < 0.0 or b < 10.0.
Returns:
the value of log(Gamma(b) / Gamma(a + b))./**
     * Returns the value of log[Γ(b) / Γ(a + b)] for a ≥ 0 and b ≥ 10. Based on
     * the <em>NSWC Library of Mathematics Subroutines</em> double precision
     * implementation, {@code DLGDIV}. In
     * {@code BetaTest.testLogGammaMinusLogGammaSum()}, this private method is
     * accessed through reflection.
     *
     * @param a First argument.
     * @param b Second argument.
     * @return the value of {@code log(Gamma(b) / Gamma(a + b))}.
     * @throws NumberIsTooSmallException if {@code a < 0.0} or {@code b < 10.0}.
     */
    private static double logGammaMinusLogGammaSum(final double a,
                                                   final double b)
        throws NumberIsTooSmallException {

        if (a < 0.0) {
            throw new NumberIsTooSmallException(a, 0.0, true);
        }
        if (b < 10.0) {
            throw new NumberIsTooSmallException(b, 10.0, true);
        }

        /*
         * d = a + b - 0.5
         */
        final double d;
        final double w;
        if (a <= b) {
            d = b + (a - 0.5);
            w = deltaMinusDeltaSum(a, b);
        } else {
            d = a + (b - 0.5);
            w = deltaMinusDeltaSum(b, a);
        }

        final double u = d * FastMath.log1p(a / b);
        final double v = a * (FastMath.log(b) - 1.0);

        return u <= v ? (w - u) - v : (w - v) - u;
    }

    
Returns the value of Δ(b) - Δ(a + b), with 0 ≤ a ≤ b and b ≥ 10. Based
on equations (26), (27) and (28) in Didonato and Morris (1992).
Params:
a – First argument.
b – Second argument.
Throws:
OutOfRangeException – if a < 0 or a > b
NumberIsTooSmallException – if b < 10
Returns:
the value of Delta(b) - Delta(a + b)/**
     * Returns the value of Δ(b) - Δ(a + b), with 0 ≤ a ≤ b and b ≥ 10. Based
     * on equations (26), (27) and (28) in Didonato and Morris (1992).
     *
     * @param a First argument.
     * @param b Second argument.
     * @return the value of {@code Delta(b) - Delta(a + b)}
     * @throws OutOfRangeException if {@code a < 0} or {@code a > b}
     * @throws NumberIsTooSmallException if {@code b < 10}
     */
    private static double deltaMinusDeltaSum(final double a,
                                             final double b)
        throws OutOfRangeException, NumberIsTooSmallException {

        if ((a < 0) || (a > b)) {
            throw new OutOfRangeException(a, 0, b);
        }
        if (b < 10) {
            throw new NumberIsTooSmallException(b, 10, true);
        }

        final double h = a / b;
        final double p = h / (1.0 + h);
        final double q = 1.0 / (1.0 + h);
        final double q2 = q * q;
        /*
         * s[i] = 1 + q + ... - q**(2 * i)
         */
        final double[] s = new double[DELTA.length];
        s[0] = 1.0;
        for (int i = 1; i < s.length; i++) {
            s[i] = 1.0 + (q + q2 * s[i - 1]);
        }
        /*
         * w = Delta(b) - Delta(a + b)
         */
        final double sqrtT = 10.0 / b;
        final double t = sqrtT * sqrtT;
        double w = DELTA[DELTA.length - 1] * s[s.length - 1];
        for (int i = DELTA.length - 2; i >= 0; i--) {
            w = t * w + DELTA[i] * s[i];
        }
        return w * p / b;
    }

    
Returns the value of Δ(p) + Δ(q) - Δ(p + q), with p, q ≥ 10. Based on
the NSWC Library of Mathematics Subroutines double precision implementation, DBCORR. In BetaTest.testSumDeltaMinusDeltaSum(), this private method is accessed through reflection. 
Params:
p – First argument.
q – Second argument.
Throws:
NumberIsTooSmallException – if p < 10.0 or q < 10.0.
Returns:
the value of Delta(p) + Delta(q) - Delta(p + q)./**
     * Returns the value of Δ(p) + Δ(q) - Δ(p + q), with p, q ≥ 10. Based on
     * the <em>NSWC Library of Mathematics Subroutines</em> double precision
     * implementation, {@code DBCORR}. In
     * {@code BetaTest.testSumDeltaMinusDeltaSum()}, this private method is
     * accessed through reflection.
     *
     * @param p First argument.
     * @param q Second argument.
     * @return the value of {@code Delta(p) + Delta(q) - Delta(p + q)}.
     * @throws NumberIsTooSmallException if {@code p < 10.0} or {@code q < 10.0}.
     */
    private static double sumDeltaMinusDeltaSum(final double p,
                                                final double q) {

        if (p < 10.0) {
            throw new NumberIsTooSmallException(p, 10.0, true);
        }
        if (q < 10.0) {
            throw new NumberIsTooSmallException(q, 10.0, true);
        }

        final double a = FastMath.min(p, q);
        final double b = FastMath.max(p, q);
        final double sqrtT = 10.0 / a;
        final double t = sqrtT * sqrtT;
        double z = DELTA[DELTA.length - 1];
        for (int i = DELTA.length - 2; i >= 0; i--) {
            z = t * z + DELTA[i];
        }
        return z / a + deltaMinusDeltaSum(a, b);
    }

    
Returns the value of log B(p, q) for 0 ≤ x ≤ 1 and p, q > 0. Based on the
NSWC Library of Mathematics Subroutines implementation, DBETLN. 
Params:
p – First argument.
q – Second argument.
Returns:
the value of log(Beta(p, q)), NaN if p <= 0 or q <= 0./**
     * Returns the value of log B(p, q) for 0 ≤ x ≤ 1 and p, q > 0. Based on the
     * <em>NSWC Library of Mathematics Subroutines</em> implementation,
     * {@code DBETLN}.
     *
     * @param p First argument.
     * @param q Second argument.
     * @return the value of {@code log(Beta(p, q))}, {@code NaN} if
     * {@code p <= 0} or {@code q <= 0}.
     */
    public static double logBeta(final double p, final double q) {
        if (Double.isNaN(p) || Double.isNaN(q) || (p <= 0.0) || (q <= 0.0)) {
            return Double.NaN;
        }

        final double a = FastMath.min(p, q);
        final double b = FastMath.max(p, q);
        if (a >= 10.0) {
            final double w = sumDeltaMinusDeltaSum(a, b);
            final double h = a / b;
            final double c = h / (1.0 + h);
            final double u = -(a - 0.5) * FastMath.log(c);
            final double v = b * FastMath.log1p(h);
            if (u <= v) {
                return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - u) - v;
            } else {
                return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - v) - u;
            }
        } else if (a > 2.0) {
            if (b > 1000.0) {
                final int n = (int) FastMath.floor(a - 1.0);
                double prod = 1.0;
                double ared = a;
                for (int i = 0; i < n; i++) {
                    ared -= 1.0;
                    prod *= ared / (1.0 + ared / b);
                }
                return (FastMath.log(prod) - n * FastMath.log(b)) +
                        (Gamma.logGamma(ared) +
                         logGammaMinusLogGammaSum(ared, b));
            } else {
                double prod1 = 1.0;
                double ared = a;
                while (ared > 2.0) {
                    ared -= 1.0;
                    final double h = ared / b;
                    prod1 *= h / (1.0 + h);
                }
                if (b < 10.0) {
                    double prod2 = 1.0;
                    double bred = b;
                    while (bred > 2.0) {
                        bred -= 1.0;
                        prod2 *= bred / (ared + bred);
                    }
                    return FastMath.log(prod1) +
                           FastMath.log(prod2) +
                           (Gamma.logGamma(ared) +
                           (Gamma.logGamma(bred) -
                            logGammaSum(ared, bred)));
                } else {
                    return FastMath.log(prod1) +
                           Gamma.logGamma(ared) +
                           logGammaMinusLogGammaSum(ared, b);
                }
            }
        } else if (a >= 1.0) {
            if (b > 2.0) {
                if (b < 10.0) {
                    double prod = 1.0;
                    double bred = b;
                    while (bred > 2.0) {
                        bred -= 1.0;
                        prod *= bred / (a + bred);
                    }
                    return FastMath.log(prod) +
                           (Gamma.logGamma(a) +
                            (Gamma.logGamma(bred) -
                             logGammaSum(a, bred)));
                } else {
                    return Gamma.logGamma(a) +
                           logGammaMinusLogGammaSum(a, b);
                }
            } else {
                return Gamma.logGamma(a) +
                       Gamma.logGamma(b) -
                       logGammaSum(a, b);
            }
        } else {
            if (b >= 10.0) {
                return Gamma.logGamma(a) +
                       logGammaMinusLogGammaSum(a, b);
            } else {
                // The following command is the original NSWC implementation.
                // return Gamma.logGamma(a) +
                // (Gamma.logGamma(b) - Gamma.logGamma(a + b));
                // The following command turns out to be more accurate.
                return FastMath.log(Gamma.gamma(a) * Gamma.gamma(b) /
                                    Gamma.gamma(a + b));
            }
        }
    }
}