Returns the logarithm of the Beta function. - Java java.lang

Description

Returns the logarithm of the Beta function.

Demo Code

/*/*from   www . j a  v  a 2 s.com*/
 * ------------------------------------------------------------------------=
 * Copyright (C) 1997 - 1998 by Visual Numerics, Inc. All rights reserved.
 *
 * Permission to use, copy, modify, and distribute this software is freely
 * granted by Visual Numerics, Inc., provided that the copyright notice
 * above and the following warranty disclaimer are preserved in human
 * readable form.
 *
 * Because this software is licenses free of charge, it is provided
 * "AS IS", with NO WARRANTY.  TO THE EXTENT PERMITTED BY LAW, VNI
 * DISCLAIMS ALL WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
 * TO ITS PERFORMANCE, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 * VNI WILL NOT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING OUT OF THE USE
 * OF OR INABILITY TO USE THIS SOFTWARE, INCLUDING BUT NOT LIMITED TO DIRECT,
 * INDIRECT, SPECIAL, CONSEQUENTIAL, PUNITIVE, AND EXEMPLARY DAMAGES, EVEN
 * IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGES. 
 *
 * ------------------------------------------------------------------------=
 */

public class Main{
    private static final double GAMMA_COEF[] = {
            .8571195590989331421920062399942e-2,
            .4415381324841006757191315771652e-2,
            .5685043681599363378632664588789e-1,
            -.4219835396418560501012500186624e-2,
            .1326808181212460220584006796352e-2,
            -.1893024529798880432523947023886e-3,
            .3606925327441245256578082217225e-4,
            -.6056761904460864218485548290365e-5,
            .1055829546302283344731823509093e-5,
            -.1811967365542384048291855891166e-6,
            .3117724964715322277790254593169e-7,
            -.5354219639019687140874081024347e-8,
            .9193275519859588946887786825940e-9,
            -.1577941280288339761767423273953e-9,
            .2707980622934954543266540433089e-10,
            -.4646818653825730144081661058933e-11,
            .7973350192007419656460767175359e-12,
            -.1368078209830916025799499172309e-12,
            .2347319486563800657233471771688e-13,
            -.4027432614949066932766570534699e-14,
            .6910051747372100912138336975257e-15,
            -.1185584500221992907052387126192e-15,
            .2034148542496373955201026051932e-16,
            -.3490054341717405849274012949108e-17,
            .5987993856485305567135051066026e-18,
            -.1027378057872228074490069778431e-18 };
    private static final double R9LGMC_COEF[] = {
            .166638948045186324720572965082e0,
            -.138494817606756384073298605914e-4,
            .981082564692472942615717154749e-8,
            -.180912947557249419426330626672e-10,
            .622109804189260522712601554342e-13,
            -.339961500541772194430333059967e-15,
            .268318199848269874895753884667e-17 };
    final private static double ALNRCS_COEF[] = {
            .103786935627437698006862677191e1,
            -.133643015049089180987660415531,
            .194082491355205633579261993748e-1,
            -.301075511275357776903765377766e-2,
            .486946147971548500904563665091e-3,
            -.810548818931753560668099430086e-4,
            .137788477995595247829382514961e-4,
            -.238022108943589702513699929149e-5,
            .41640416213865183476391859902e-6,
            -.73595828378075994984266837032e-7,
            .13117611876241674949152294345e-7,
            -.235467093177424251366960923302e-8,
            .425227732760349977756380529626e-9,
            -.771908941348407968261081074933e-10,
            .140757464813590699092153564722e-10,
            -.257690720580246806275370786276e-11,
            .473424066662944218491543950059e-12,
            -.872490126747426417453012632927e-13,
            .161246149027405514657398331191e-13,
            -.298756520156657730067107924168e-14,
            .554807012090828879830413216973e-15,
            -.103246191582715695951413339619e-15,
            .192502392030498511778785032449e-16,
            -.359550734652651500111897078443e-17,
            .672645425378768578921945742268e-18,
            -.126026241687352192520824256376e-18 };
    /**
     * Returns the logarithm of the Beta function.
     * 
     * @param a
     *            A double value.
     * @param b
     *            A double value.
     * @return The natural logarithm of the Beta function.
     */
    static public double logBeta(double a, double b) {
        double corr, ans;
        double p = Math.min(a, b);
        double q = Math.max(a, b);

        if (p <= 0.0) {
            ans = Double.NaN;
        } else if (p >= 10.0) {
            // P and Q are large;
            corr = r9lgmc(p) + r9lgmc(q) - r9lgmc(p + q);
            double temp = dlnrel(-p / (p + q));
            ans = -0.5 * Math.log(q) + 0.918938533204672741780329736406
                    + corr + (p - 0.5) * Math.log(p / (p + q)) + q * temp;
        } else if (q >= 10.0) {
            // P is small, but Q is large
            corr = MathUtils.r9lgmc(q) - r9lgmc(p + q);
            // Check from underflow from r9lgmc
            ans = logGamma(p) + corr + p - p * Math.log(p + q) + (q - 0.5)
                    * dlnrel(-p / (p + q));
        } else {
            // P and Q are small;
            ans = Math.log(gamma(p) * (gamma(q) / gamma(p + q)));
        }
        return ans;
    }
    static double r9lgmc(double x) {
        double ans;

        if (x < 10.0) {
            ans = Double.NaN;
        } else if (x < 9.490626562e+07) {
            // 9.490626562e+07 = 1/Math.sqrt(EPSILON_SMALL)
            double y = 10.0 / x;
            ans = csevl(2.0 * y * y - 1.0, R9LGMC_COEF) / x;
        } else if (x < 1.39118e+11) {
            // 1.39118e+11 = exp(min(log(amach(2) / 12.0), -log(12.0 *
            // amach(1))));
            // See A&S 6.1.41
            ans = 1.0 / (12.0 * x);
        } else {
            ans = 0.0; // underflows
        }
        return ans;
    }
    private static double dlnrel(double x) {
        double ans;

        if (x <= -1.0) {
            ans = Double.NaN;
        } else if (Math.abs(x) <= 0.375) {
            ans = x * (1.0 - x * MathUtils.csevl(x / .375, ALNRCS_COEF));
        } else {
            ans = Math.log(1.0 + x);
        }
        return ans;
    }
    /**
     * Returns the logarithm of the Gamma function of a double.
     * 
     * @param x
     *            A double value.
     * @return The natural logarithm of the Gamma function of x. If x is a
     *         negative integer, the result is NaN.
     */
    static public double logGamma(double x) {
        double ans, sinpiy, y;

        y = Math.abs(x);

        if (y <= 10) {
            ans = Math.log(Math.abs(gamma(x)));
        } else if (x > 0) {
            // A&S 6.1.40
            // 0.9189385332046727 = 0.5*log(2*PI)
            ans = 0.9189385332046727 + (x - 0.5) * Math.log(x) - x
                    + r9lgmc(y);
        } else {
            sinpiy = Math.abs(Math.sin(Math.PI * y));
            if (sinpiy == 0 || Math.round(y) == y) {
                // The argument for the function can not be a negative integer.
                ans = Double.NaN;
            } else {
                ans = 0.22579135264472743236 + (x - 0.5) * Math.log(y) - x
                        - Math.log(sinpiy) - r9lgmc(y);
            }
        }
        return ans;
    }
    /**
     * Returns the Gamma function of a double.
     * 
     * @param x
     *            A double value.
     * @return The Gamma function of x. If x is a negative integer, the result
     *         is NaN.
     */
    static public double gamma(double x) {
        double ans;
        double y = Math.abs(x);

        if (y <= 10.0) {
            /*
             * Compute gamma(x) for |x|<=10. First reduce the interval and find
             * gamma(1+y) for 0 <= y < 1.
             */
            int n = (int) x;
            if (x < 0.0)
                n--;
            y = x - n;
            n--;
            ans = 0.9375 + csevl(2.0 * y - 1.0, GAMMA_COEF);
            if (n == 0) {
            } else if (n < 0) {
                // Compute gamma(x) for x < 1
                n = -n;
                if (x == 0.0) {
                    ans = Double.NaN;
                } else if (y < 1.0 / Double.MAX_VALUE) {
                    ans = Double.POSITIVE_INFINITY;
                } else {
                    double xn = n - 2;
                    if (x < 0.0 && x + xn == 0.0) {
                        ans = Double.NaN;
                    } else {
                        for (int i = 0; i < n; i++) {
                            ans /= x + i;
                        }
                    }
                }
            } else { // gamma(x) for x >= 2.0
                for (int i = 1; i <= n; i++) {
                    ans *= y + i;
                }
            }
        } else { // gamma(x) for |x| > 10
            if (x > 171.614) {
                ans = Double.POSITIVE_INFINITY;
            } else if (x < -170.56) {
                ans = 0.0; // underflows
            } else {
                // 0.9189385332046727 = 0.5*log(2*PI)
                ans = Math.exp((y - 0.5) * Math.log(y) - y
                        + 0.9189385332046727 + r9lgmc(y));
                if (x < 0.0) {
                    double sinpiy = Math.sin(Math.PI * y);
                    if (sinpiy == 0 || Math.round(y) == y) {
                        ans = Double.NaN;
                    } else {
                        ans = -Math.PI / (y * sinpiy * ans);
                    }
                }
            }
        }
        return ans;
    }
    static double csevl(double x, double coef[]) {
        double b0, b1, b2, twox;
        int i;
        b1 = 0.0;
        b0 = 0.0;
        b2 = 0.0;
        twox = 2.0 * x;
        for (i = coef.length - 1; i >= 0; i--) {
            b2 = b1;
            b1 = b0;
            b0 = twox * b1 - b2 + coef[i];
        }
        return 0.5 * (b0 - b2);
    }
    /**
     * Rounds value.
     *
     * @param value the number
     * @param decimalplaces number of decimal places
     * @return rounded value as string
     */
    static public String round(double value, double decimalplaces) {

        double mult = Math.pow(10.0, decimalplaces);

        return String.valueOf(Math.round(value * mult) / mult);
    }
}

• Previous
• Next

Related Tutorials

• Compute an approximation to the natural logarithm.
• Compute logFactorial(N).
• Returns the common (base 10) logarithm of a double.
• Returns the logarithm of the Gamma function of a double.
• Gets the gaussian pdf value log.