| /* |
| * 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.numbers.gamma; |
| |
| /** |
| * Computes \( log_e \Beta(p, q) \). |
| * <p> |
| * This class is immutable. |
| * </p> |
| */ |
| public final class LogBeta { |
| /** The threshold value of 10 where the series expansion of the Δ function applies. */ |
| private static final double TEN = 10; |
| /** The threshold value of 2 for algorithm switch. */ |
| private static final double TWO = 2; |
| /** The threshold value of 1000 for algorithm switch. */ |
| private static final double THOUSAND = 1000; |
| |
| /** The constant value of ½log 2π. */ |
| private static final double HALF_LOG_TWO_PI = 0.9189385332046727; |
| |
| /** |
| * The coefficients of the series expansion of the Δ function. This function |
| * is defined as follows: |
| * <pre> |
| * Δ(x) = log Γ(x) - (x - 0.5) log a + a - 0.5 log 2π, |
| * </pre> |
| * <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 |
| }; |
| |
| /** Private constructor. */ |
| private LogBeta() { |
| // intentional empty. |
| } |
| |
| /** |
| * 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 IllegalArgumentException if {@code a < 0} or {@code a > b} |
| * @throws IllegalArgumentException if {@code b < 10} |
| */ |
| private static double deltaMinusDeltaSum(final double a, |
| final double b) { |
| if (a < 0 || |
| a > b) { |
| throw new GammaException(GammaException.OUT_OF_RANGE, a, 0, b); |
| } |
| if (b < TEN) { |
| throw new GammaException(GammaException.OUT_OF_RANGE, b, TEN, Double.POSITIVE_INFINITY); |
| } |
| |
| final double h = a / b; |
| final double p = h / (1 + h); |
| final double q = 1 / (1 + h); |
| final double q2 = q * q; |
| /* |
| * s[i] = 1 + q + ... - q**(2 * i) |
| */ |
| final double[] s = new double[DELTA.length]; |
| s[0] = 1; |
| for (int i = 1; i < s.length; i++) { |
| s[i] = 1 + (q + q2 * s[i - 1]); |
| } |
| /* |
| * w = Delta(b) - Delta(a + b) |
| */ |
| final double sqrtT = 10 / 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 <em>NSWC Library of Mathematics Subroutines</em> implementation, |
| * {@code DBCORR}. |
| * |
| * @param p First argument. |
| * @param q Second argument. |
| * @return the value of {@code Delta(p) + Delta(q) - Delta(p + q)}. |
| * @throws IllegalArgumentException if {@code p < 10} or {@code q < 10}. |
| */ |
| private static double sumDeltaMinusDeltaSum(final double p, |
| final double q) { |
| |
| if (p < TEN) { |
| throw new GammaException(GammaException.OUT_OF_RANGE, p, TEN, Double.POSITIVE_INFINITY); |
| } |
| if (q < TEN) { |
| throw new GammaException(GammaException.OUT_OF_RANGE, q, TEN, Double.POSITIVE_INFINITY); |
| } |
| |
| final double a = Math.min(p, q); |
| final double b = Math.max(p, q); |
| final double sqrtT = 10 / 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 {@code log B(p, q)} for {@code 0 ≤ x ≤ 1} and {@code 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 value(double p, |
| double q) { |
| if (Double.isNaN(p) || |
| Double.isNaN(q) || |
| p <= 0 || |
| q <= 0) { |
| return Double.NaN; |
| } |
| |
| final double a = Math.min(p, q); |
| final double b = Math.max(p, q); |
| if (a >= TEN) { |
| final double w = sumDeltaMinusDeltaSum(a, b); |
| final double h = a / b; |
| final double c = h / (1 + h); |
| final double u = -(a - 0.5) * Math.log(c); |
| final double v = b * Math.log1p(h); |
| if (u <= v) { |
| return (((-0.5 * Math.log(b) + HALF_LOG_TWO_PI) + w) - u) - v; |
| } else { |
| return (((-0.5 * Math.log(b) + HALF_LOG_TWO_PI) + w) - v) - u; |
| } |
| } else if (a > TWO) { |
| if (b > THOUSAND) { |
| final int n = (int) Math.floor(a - 1); |
| double prod = 1; |
| double ared = a; |
| for (int i = 0; i < n; i++) { |
| ared -= 1; |
| prod *= ared / (1 + ared / b); |
| } |
| return (Math.log(prod) - n * Math.log(b)) + |
| (LogGamma.value(ared) + |
| logGammaMinusLogGammaSum(ared, b)); |
| } else { |
| double prod1 = 1; |
| double ared = a; |
| while (ared > 2) { |
| ared -= 1; |
| final double h = ared / b; |
| prod1 *= h / (1 + h); |
| } |
| if (b < TEN) { |
| double prod2 = 1; |
| double bred = b; |
| while (bred > 2) { |
| bred -= 1; |
| prod2 *= bred / (ared + bred); |
| } |
| return Math.log(prod1) + |
| Math.log(prod2) + |
| (LogGamma.value(ared) + |
| (LogGamma.value(bred) - |
| LogGammaSum.value(ared, bred))); |
| } else { |
| return Math.log(prod1) + |
| LogGamma.value(ared) + |
| logGammaMinusLogGammaSum(ared, b); |
| } |
| } |
| } else if (a >= 1) { |
| if (b > TWO) { |
| if (b < TEN) { |
| double prod = 1; |
| double bred = b; |
| while (bred > 2) { |
| bred -= 1; |
| prod *= bred / (a + bred); |
| } |
| return Math.log(prod) + |
| (LogGamma.value(a) + |
| (LogGamma.value(bred) - |
| LogGammaSum.value(a, bred))); |
| } else { |
| return LogGamma.value(a) + |
| logGammaMinusLogGammaSum(a, b); |
| } |
| } else { |
| return LogGamma.value(a) + |
| LogGamma.value(b) - |
| LogGammaSum.value(a, b); |
| } |
| } else { |
| if (b >= TEN) { |
| return LogGamma.value(a) + |
| logGammaMinusLogGammaSum(a, b); |
| } else { |
| // The original NSWC implementation was |
| // LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b)); |
| // but the following command turned out to be more accurate. |
| return Math.log(Gamma.value(a) * Gamma.value(b) / |
| Gamma.value(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> implementation, |
| * {@code DLGDIV}. |
| * |
| * @param a First argument. |
| * @param b Second argument. |
| * @return the value of {@code log(Gamma(b) / Gamma(a + b))}. |
| * @throws IllegalArgumentException if {@code a < 0} or {@code b < 10}. |
| */ |
| private static double logGammaMinusLogGammaSum(double a, |
| double b) { |
| if (a < 0) { |
| throw new GammaException(GammaException.OUT_OF_RANGE, a, 0, Double.POSITIVE_INFINITY); |
| } |
| if (b < TEN) { |
| throw new GammaException(GammaException.OUT_OF_RANGE, b, TEN, Double.POSITIVE_INFINITY); |
| } |
| |
| /* |
| * 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 * Math.log1p(a / b); |
| final double v = a * (Math.log(b) - 1); |
| |
| return u <= v ? |
| (w - u) - v : |
| (w - v) - u; |
| } |
| } |
