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 \Beta(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 Δ 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 Δ function. This function
038     * is defined as follows:
039     * <pre>
040     * Δ(x) = log Γ(x) - (x - 0.5) log a + a - 0.5 log 2π,
041     * </pre>
042     * <p>
043     * See equation (23) in Didonato and Morris (1992). The series expansion,
044     * which applies for x ≥ 10, reads
045     * </p>
046     * <pre>
047     *                 14
048     *                ====
049     *             1  \                2 n
050     *     Δ(x) = ---  >    d  (10 / x)
051     *             x  /      n
052     *                ====
053     *                n = 0
054     * </pre>
055     */
056    private static final double[] DELTA = {
057        .833333333333333333333333333333E-01,
058        -.277777777777777777777777752282E-04,
059        .793650793650793650791732130419E-07,
060        -.595238095238095232389839236182E-09,
061        .841750841750832853294451671990E-11,
062        -.191752691751854612334149171243E-12,
063        .641025640510325475730918472625E-14,
064        -.295506514125338232839867823991E-15,
065        .179643716359402238723287696452E-16,
066        -.139228964661627791231203060395E-17,
067        .133802855014020915603275339093E-18,
068        -.154246009867966094273710216533E-19,
069        .197701992980957427278370133333E-20,
070        -.234065664793997056856992426667E-21,
071        .171348014966398575409015466667E-22
072    };
073
074    /** Private constructor. */
075    private LogBeta() {
076        // intentional empty.
077    }
078
079    /**
080     * Returns the value of Δ(b) - Δ(a + b), with 0 ≤ a ≤ b and b ≥ 10. Based
081     * on equations (26), (27) and (28) in Didonato and Morris (1992).
082     *
083     * @param a First argument.
084     * @param b Second argument.
085     * @return the value of {@code Delta(b) - Delta(a + b)}
086     * @throws IllegalArgumentException if {@code a < 0} or {@code a > b}
087     * @throws IllegalArgumentException if {@code b < 10}
088     */
089    private static double deltaMinusDeltaSum(final double a,
090                                             final double b) {
091        if (a < 0 ||
092            a > b) {
093            throw new GammaException(GammaException.OUT_OF_RANGE, a, 0, b);
094        }
095        if (b < TEN) {
096            throw new GammaException(GammaException.OUT_OF_RANGE, b, TEN, Double.POSITIVE_INFINITY);
097        }
098
099        final double h = a / b;
100        final double p = h / (1 + h);
101        final double q = 1 / (1 + h);
102        final double q2 = q * q;
103        /*
104         * s[i] = 1 + q + ... - q**(2 * i)
105         */
106        final double[] s = new double[DELTA.length];
107        s[0] = 1;
108        for (int i = 1; i < s.length; i++) {
109            s[i] = 1 + (q + q2 * s[i - 1]);
110        }
111        /*
112         * w = Delta(b) - Delta(a + b)
113         */
114        final double sqrtT = 10 / b;
115        final double t = sqrtT * sqrtT;
116        double w = DELTA[DELTA.length - 1] * s[s.length - 1];
117        for (int i = DELTA.length - 2; i >= 0; i--) {
118            w = t * w + DELTA[i] * s[i];
119        }
120        return w * p / b;
121    }
122
123    /**
124     * Returns the value of Δ(p) + Δ(q) - Δ(p + q), with p, q ≥ 10.
125     * Based on the <em>NSWC Library of Mathematics Subroutines</em> implementation,
126     * {@code DBCORR}.
127     *
128     * @param p First argument.
129     * @param q Second argument.
130     * @return the value of {@code Delta(p) + Delta(q) - Delta(p + q)}.
131     * @throws IllegalArgumentException if {@code p < 10} or {@code q < 10}.
132     */
133    private static double sumDeltaMinusDeltaSum(final double p,
134                                                final double q) {
135
136        if (p < TEN) {
137            throw new GammaException(GammaException.OUT_OF_RANGE, p, TEN, Double.POSITIVE_INFINITY);
138        }
139        if (q < TEN) {
140            throw new GammaException(GammaException.OUT_OF_RANGE, q, TEN, Double.POSITIVE_INFINITY);
141        }
142
143        final double a = Math.min(p, q);
144        final double b = Math.max(p, q);
145        final double sqrtT = 10 / a;
146        final double t = sqrtT * sqrtT;
147        double z = DELTA[DELTA.length - 1];
148        for (int i = DELTA.length - 2; i >= 0; i--) {
149            z = t * z + DELTA[i];
150        }
151        return z / a + deltaMinusDeltaSum(a, b);
152    }
153
154    /**
155     * Returns the value of {@code log B(p, q)} for {@code 0 ≤ x ≤ 1} and {@code p, q > 0}.
156     * Based on the <em>NSWC Library of Mathematics Subroutines</em> implementation,
157     * {@code DBETLN}.
158     *
159     * @param p First argument.
160     * @param q Second argument.
161     * @return the value of {@code log(Beta(p, q))}, {@code NaN} if
162     * {@code p <= 0} or {@code q <= 0}.
163     */
164    public static double value(double p,
165                               double q) {
166        if (Double.isNaN(p) ||
167            Double.isNaN(q) ||
168            p <= 0 ||
169            q <= 0) {
170            return Double.NaN;
171        }
172
173        final double a = Math.min(p, q);
174        final double b = Math.max(p, q);
175        if (a >= TEN) {
176            final double w = sumDeltaMinusDeltaSum(a, b);
177            final double h = a / b;
178            final double c = h / (1 + h);
179            final double u = -(a - 0.5) * Math.log(c);
180            final double v = b * Math.log1p(h);
181            if (u <= v) {
182                return (((-0.5 * Math.log(b) + HALF_LOG_TWO_PI) + w) - u) - v;
183            } else {
184                return (((-0.5 * Math.log(b) + HALF_LOG_TWO_PI) + w) - v) - u;
185            }
186        } else if (a > TWO) {
187            if (b > THOUSAND) {
188                final int n = (int) Math.floor(a - 1);
189                double prod = 1;
190                double ared = a;
191                for (int i = 0; i < n; i++) {
192                    ared -= 1;
193                    prod *= ared / (1 + ared / b);
194                }
195                return (Math.log(prod) - n * Math.log(b)) +
196                        (LogGamma.value(ared) +
197                         logGammaMinusLogGammaSum(ared, b));
198            } else {
199                double prod1 = 1;
200                double ared = a;
201                while (ared > 2) {
202                    ared -= 1;
203                    final double h = ared / b;
204                    prod1 *= h / (1 + h);
205                }
206                if (b < TEN) {
207                    double prod2 = 1;
208                    double bred = b;
209                    while (bred > 2) {
210                        bred -= 1;
211                        prod2 *= bred / (ared + bred);
212                    }
213                    return Math.log(prod1) +
214                           Math.log(prod2) +
215                           (LogGamma.value(ared) +
216                           (LogGamma.value(bred) -
217                            LogGammaSum.value(ared, bred)));
218                } else {
219                    return Math.log(prod1) +
220                           LogGamma.value(ared) +
221                           logGammaMinusLogGammaSum(ared, b);
222                }
223            }
224        } else if (a >= 1) {
225            if (b > TWO) {
226                if (b < TEN) {
227                    double prod = 1;
228                    double bred = b;
229                    while (bred > 2) {
230                        bred -= 1;
231                        prod *= bred / (a + bred);
232                    }
233                    return Math.log(prod) +
234                           (LogGamma.value(a) +
235                            (LogGamma.value(bred) -
236                             LogGammaSum.value(a, bred)));
237                } else {
238                    return LogGamma.value(a) +
239                        logGammaMinusLogGammaSum(a, b);
240                }
241            } else {
242                return LogGamma.value(a) +
243                       LogGamma.value(b) -
244                       LogGammaSum.value(a, b);
245            }
246        } else {
247            if (b >= TEN) {
248                return LogGamma.value(a) +
249                       logGammaMinusLogGammaSum(a, b);
250            } else {
251                // The original NSWC implementation was
252                //   LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b));
253                // but the following command turned out to be more accurate.
254                return Math.log(Gamma.value(a) * Gamma.value(b) /
255                                Gamma.value(a + b));
256            }
257        }
258    }
259
260    /**
261     * Returns the value of log[Γ(b) / Γ(a + b)] for a ≥ 0 and b ≥ 10.
262     * Based on the <em>NSWC Library of Mathematics Subroutines</em> implementation,
263     * {@code DLGDIV}.
264     *
265     * @param a First argument.
266     * @param b Second argument.
267     * @return the value of {@code log(Gamma(b) / Gamma(a + b))}.
268     * @throws IllegalArgumentException if {@code a < 0} or {@code b < 10}.
269     */
270    private static double logGammaMinusLogGammaSum(double a,
271                                                   double b) {
272        if (a < 0) {
273            throw new GammaException(GammaException.OUT_OF_RANGE, a, 0, Double.POSITIVE_INFINITY);
274        }
275        if (b < TEN) {
276            throw new GammaException(GammaException.OUT_OF_RANGE, b, TEN, Double.POSITIVE_INFINITY);
277        }
278
279        /*
280         * d = a + b - 0.5
281         */
282        final double d;
283        final double w;
284        if (a <= b) {
285            d = b + (a - 0.5);
286            w = deltaMinusDeltaSum(a, b);
287        } else {
288            d = a + (b - 0.5);
289            w = deltaMinusDeltaSum(b, a);
290        }
291
292        final double u = d * Math.log1p(a / b);
293        final double v = a * (Math.log(b) - 1);
294
295        return u <= v ?
296            (w - u) - v :
297            (w - v) - u;
298    }
299}