Source code

001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.numbers.gamma;
018
019/**
020 * Computes \( log_e B(p, q) \).
021 * <p>
022 * This class is immutable.
023 * </p>
024 */
025public final class LogBeta {
026    /** The threshold value of 10 where the series expansion of the \( \Delta \) function applies. */
027    private static final double TEN = 10;
028    /** The threshold value of 2 for algorithm switch. */
029    private static final double TWO = 2;
030    /** The threshold value of 1000 for algorithm switch. */
031    private static final double THOUSAND = 1000;
032
033    /** The constant value of ½log 2π. */
034    private static final double HALF_LOG_TWO_PI = 0.9189385332046727;
035
036    /**
037     * The coefficients of the series expansion of the \( \Delta \) function.
038     * This function is defined as follows:
039     * \[
040     *  \Delta(x) = \log \Gamma(x) - (x - \frac{1}{2}) \log a + a - \frac{1}{2} \log 2\pi,
041     * \]
042     * <p>
043     * See equation (23) in Didonato and Morris (1992). The series expansion,
044     * which applies for \( x \geq 10 \), reads
045     * </p>
046     * \[
047     *  \Delta(x) = \frac{1}{x} \sum_{n = 0}^{14} d_n (\frac{10}{x})^{2 n}
048     * \]
049     */
050    private static final double[] DELTA = {
051        .833333333333333333333333333333E-01,
052        -.277777777777777777777777752282E-04,
053        .793650793650793650791732130419E-07,
054        -.595238095238095232389839236182E-09,
055        .841750841750832853294451671990E-11,
056        -.191752691751854612334149171243E-12,
057        .641025640510325475730918472625E-14,
058        -.295506514125338232839867823991E-15,
059        .179643716359402238723287696452E-16,
060        -.139228964661627791231203060395E-17,
061        .133802855014020915603275339093E-18,
062        -.154246009867966094273710216533E-19,
063        .197701992980957427278370133333E-20,
064        -.234065664793997056856992426667E-21,
065        .171348014966398575409015466667E-22
066    };
067
068    /** Private constructor. */
069    private LogBeta() {
070        // intentionally empty.
071    }
072
073    /**
074     * Returns the value of \( \Delta(b) - \Delta(a + b) \),
075     * with \( 0 \leq a \leq b \) and \( b \geq 10 \).
076     * Based on equations (26), (27) and (28) in Didonato and Morris (1992).
077     *
078     * @param a First argument.
079     * @param b Second argument.
080     * @return the value of \( \Delta(b) - \Delta(a + b) \)
081     * @throws IllegalArgumentException if {@code a < 0} or {@code a > b}
082     * @throws IllegalArgumentException if {@code b < 10}
083     */
084    private static double deltaMinusDeltaSum(final double a,
085                                             final double b) {
086        // Assumptions:
087        // Never called with a < 0; a > b; or b < 10
088
089        final double h = a / b;
090        final double p = h / (1 + h);
091        final double q = 1 / (1 + h);
092        final double q2 = q * q;
093        /*
094         * s[i] = 1 + q + ... - q**(2 * i)
095         */
096        final double[] s = new double[DELTA.length];
097        s[0] = 1;
098        for (int i = 1; i < s.length; i++) {
099            s[i] = 1 + (q + q2 * s[i - 1]);
100        }
101        /*
102         * w = Delta(b) - Delta(a + b)
103         */
104        final double sqrtT = 10 / b;
105        final double t = sqrtT * sqrtT;
106        double w = DELTA[DELTA.length - 1] * s[s.length - 1];
107        for (int i = DELTA.length - 2; i >= 0; i--) {
108            w = t * w + DELTA[i] * s[i];
109        }
110        return w * p / b;
111    }
112
113    /**
114     * Returns the value of \( \Delta(p) + \Delta(q) - \Delta(p + q) \),
115     * with \( p, q \geq 10 \).
116     * Based on the <em>NSWC Library of Mathematics Subroutines</em> implementation,
117     * {@code DBCORR}.
118     *
119     * @param p First argument.
120     * @param q Second argument.
121     * @return the value of \( \Delta(p) + \Delta(q) - \Delta(p + q) \).
122     * @throws IllegalArgumentException if {@code p < 10} or {@code q < 10}.
123     */
124    private static double sumDeltaMinusDeltaSum(final double p,
125                                                final double q) {
126
127        if (p < TEN) {
128            throw new GammaException(GammaException.OUT_OF_RANGE, p, TEN, Double.POSITIVE_INFINITY);
129        }
130        if (q < TEN) {
131            throw new GammaException(GammaException.OUT_OF_RANGE, q, TEN, Double.POSITIVE_INFINITY);
132        }
133
134        final double a = Math.min(p, q);
135        final double b = Math.max(p, q);
136        final double sqrtT = 10 / a;
137        final double t = sqrtT * sqrtT;
138        double z = DELTA[DELTA.length - 1];
139        for (int i = DELTA.length - 2; i >= 0; i--) {
140            z = t * z + DELTA[i];
141        }
142        return z / a + deltaMinusDeltaSum(a, b);
143    }
144
145    /**
146     * Returns the value of \( \log B(p, q) \) for \( 0 \leq x \leq 1 \) and \( p, q &gt; 0 \).
147     * Based on the <em>NSWC Library of Mathematics Subroutines</em> implementation,
148     * {@code DBETLN}.
149     *
150     * @param p First argument.
151     * @param q Second argument.
152     * @return the value of \( \log B(p, q) \), or {@code NaN} if
153     * {@code p <= 0} or {@code q <= 0}.
154     */
155    public static double value(double p,
156                               double q) {
157        if (Double.isNaN(p) ||
158            Double.isNaN(q) ||
159            p <= 0 ||
160            q <= 0) {
161            return Double.NaN;
162        }
163
164        final double a = Math.min(p, q);
165        final double b = Math.max(p, q);
166        if (a >= TEN) {
167            final double w = sumDeltaMinusDeltaSum(a, b);
168            final double h = a / b;
169            final double c = h / (1 + h);
170            final double u = -(a - 0.5) * Math.log(c);
171            final double v = b * Math.log1p(h);
172            if (u <= v) {
173                return (((-0.5 * Math.log(b) + HALF_LOG_TWO_PI) + w) - u) - v;
174            }
175            return (((-0.5 * Math.log(b) + HALF_LOG_TWO_PI) + w) - v) - u;
176        } else if (a > TWO) {
177            if (b > THOUSAND) {
178                final int n = (int) Math.floor(a - 1);
179                double prod = 1;
180                double ared = a;
181                for (int i = 0; i < n; i++) {
182                    ared -= 1;
183                    prod *= ared / (1 + ared / b);
184                }
185                return (Math.log(prod) - n * Math.log(b)) +
186                        (LogGamma.value(ared) +
187                         logGammaMinusLogGammaSum(ared, b));
188            }
189            double prod1 = 1;
190            double ared = a;
191            while (ared > TWO) {
192                ared -= 1;
193                final double h = ared / b;
194                prod1 *= h / (1 + h);
195            }
196            if (b < TEN) {
197                double prod2 = 1;
198                double bred = b;
199                while (bred > TWO) {
200                    bred -= 1;
201                    prod2 *= bred / (ared + bred);
202                }
203                return Math.log(prod1) +
204                       Math.log(prod2) +
205                       (LogGamma.value(ared) +
206                       (LogGamma.value(bred) -
207                        LogGammaSum.value(ared, bred)));
208            }
209            return Math.log(prod1) +
210                   LogGamma.value(ared) +
211                   logGammaMinusLogGammaSum(ared, b);
212        } else if (a >= 1) {
213            if (b > TWO) {
214                if (b < TEN) {
215                    double prod = 1;
216                    double bred = b;
217                    while (bred > TWO) {
218                        bred -= 1;
219                        prod *= bred / (a + bred);
220                    }
221                    return Math.log(prod) +
222                           (LogGamma.value(a) +
223                            (LogGamma.value(bred) -
224                             LogGammaSum.value(a, bred)));
225                }
226                return LogGamma.value(a) +
227                    logGammaMinusLogGammaSum(a, b);
228            }
229            return LogGamma.value(a) +
230                   LogGamma.value(b) -
231                   LogGammaSum.value(a, b);
232        } else {
233            if (b >= TEN) {
234                return LogGamma.value(a) +
235                       logGammaMinusLogGammaSum(a, b);
236            }
237            // The original NSWC implementation was
238            //   LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b));
239            // but the following command turned out to be more accurate.
240            // Note: Check for overflow that occurs if a and/or b are tiny.
241            final double beta = Gamma.value(a) * Gamma.value(b) / Gamma.value(a + b);
242            if (Double.isFinite(beta)) {
243                return Math.log(beta);
244            }
245            return LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b));
246        }
247    }
248
249    /**
250     * Returns the value of \( \log ( \Gamma(b) / \Gamma(a + b) ) \)
251     * for \( a \geq 0 \) and \( b \geq 10 \).
252     * Based on the <em>NSWC Library of Mathematics Subroutines</em> implementation,
253     * {@code DLGDIV}.
254     *
255     * <p>This method assumes \( a \leq b \).
256     *
257     * @param a First argument.
258     * @param b Second argument.
259     * @return the value of \( \log(\Gamma(b) / \Gamma(a + b) \).
260     * @throws IllegalArgumentException if {@code a < 0} or {@code b < 10}.
261     */
262    private static double logGammaMinusLogGammaSum(double a,
263                                                   double b) {
264        if (a < 0) {
265            throw new GammaException(GammaException.OUT_OF_RANGE, a, 0, Double.POSITIVE_INFINITY);
266        }
267        if (b < TEN) {
268            throw new GammaException(GammaException.OUT_OF_RANGE, b, TEN, Double.POSITIVE_INFINITY);
269        }
270
271        /*
272         * d = a + b - 0.5
273         */
274        // Assumption: always called with a <= b.
275        final double d = b + (a - 0.5);
276        final double w = deltaMinusDeltaSum(a, b);
277
278        final double u = d * Math.log1p(a / b);
279        final double v = a * (Math.log(b) - 1);
280
281        return u <= v ?
282            (w - u) - v :
283            (w - v) - u;
284    }
285}