1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.numbers.gamma;
18
19 /**
20 * Computes \( log_e B(p, q) \).
21 * <p>
22 * This class is immutable.
23 * </p>
24 */
25 public final class LogBeta {
26 /** The threshold value of 10 where the series expansion of the \( \Delta \) function applies. */
27 private static final double TEN = 10;
28 /** The threshold value of 2 for algorithm switch. */
29 private static final double TWO = 2;
30 /** The threshold value of 1000 for algorithm switch. */
31 private static final double THOUSAND = 1000;
32
33 /** The constant value of ½log 2π. */
34 private static final double HALF_LOG_TWO_PI = 0.9189385332046727;
35
36 /**
37 * The coefficients of the series expansion of the \( \Delta \) function.
38 * This function is defined as follows:
39 * \[
40 * \Delta(x) = \log \Gamma(x) - (x - \frac{1}{2}) \log a + a - \frac{1}{2} \log 2\pi,
41 * \]
42 * <p>
43 * See equation (23) in Didonato and Morris (1992). The series expansion,
44 * which applies for \( x \geq 10 \), reads
45 * </p>
46 * \[
47 * \Delta(x) = \frac{1}{x} \sum_{n = 0}^{14} d_n (\frac{10}{x})^{2 n}
48 * \]
49 */
50 private static final double[] DELTA = {
51 .833333333333333333333333333333E-01,
52 -.277777777777777777777777752282E-04,
53 .793650793650793650791732130419E-07,
54 -.595238095238095232389839236182E-09,
55 .841750841750832853294451671990E-11,
56 -.191752691751854612334149171243E-12,
57 .641025640510325475730918472625E-14,
58 -.295506514125338232839867823991E-15,
59 .179643716359402238723287696452E-16,
60 -.139228964661627791231203060395E-17,
61 .133802855014020915603275339093E-18,
62 -.154246009867966094273710216533E-19,
63 .197701992980957427278370133333E-20,
64 -.234065664793997056856992426667E-21,
65 .171348014966398575409015466667E-22
66 };
67
68 /** Private constructor. */
69 private LogBeta() {
70 // intentionally empty.
71 }
72
73 /**
74 * Returns the value of \( \Delta(b) - \Delta(a + b) \),
75 * with \( 0 \leq a \leq b \) and \( b \geq 10 \).
76 * Based on equations (26), (27) and (28) in Didonato and Morris (1992).
77 *
78 * @param a First argument.
79 * @param b Second argument.
80 * @return the value of \( \Delta(b) - \Delta(a + b) \)
81 * @throws IllegalArgumentException if {@code a < 0} or {@code a > b}
82 * @throws IllegalArgumentException if {@code b < 10}
83 */
84 private static double deltaMinusDeltaSum(final double a,
85 final double b) {
86 if (a < 0 ||
87 a > b) {
88 throw new GammaException(GammaException.OUT_OF_RANGE, a, 0, b);
89 }
90 if (b < TEN) {
91 throw new GammaException(GammaException.OUT_OF_RANGE, b, TEN, Double.POSITIVE_INFINITY);
92 }
93
94 final double h = a / b;
95 final double p = h / (1 + h);
96 final double q = 1 / (1 + h);
97 final double q2 = q * q;
98 /*
99 * s[i] = 1 + q + ... - q**(2 * i)
100 */
101 final double[] s = new double[DELTA.length];
102 s[0] = 1;
103 for (int i = 1; i < s.length; i++) {
104 s[i] = 1 + (q + q2 * s[i - 1]);
105 }
106 /*
107 * w = Delta(b) - Delta(a + b)
108 */
109 final double sqrtT = 10 / b;
110 final double t = sqrtT * sqrtT;
111 double w = DELTA[DELTA.length - 1] * s[s.length - 1];
112 for (int i = DELTA.length - 2; i >= 0; i--) {
113 w = t * w + DELTA[i] * s[i];
114 }
115 return w * p / b;
116 }
117
118 /**
119 * Returns the value of \( \Delta(p) + \Delta(q) - \Delta(p + q) \),
120 * with \( p, q \geq 10 \).
121 * Based on the <em>NSWC Library of Mathematics Subroutines</em> implementation,
122 * {@code DBCORR}.
123 *
124 * @param p First argument.
125 * @param q Second argument.
126 * @return the value of \( \Delta(p) + \Delta(q) - \Delta(p + q) \).
127 * @throws IllegalArgumentException if {@code p < 10} or {@code q < 10}.
128 */
129 private static double sumDeltaMinusDeltaSum(final double p,
130 final double q) {
131
132 if (p < TEN) {
133 throw new GammaException(GammaException.OUT_OF_RANGE, p, TEN, Double.POSITIVE_INFINITY);
134 }
135 if (q < TEN) {
136 throw new GammaException(GammaException.OUT_OF_RANGE, q, TEN, Double.POSITIVE_INFINITY);
137 }
138
139 final double a = Math.min(p, q);
140 final double b = Math.max(p, q);
141 final double sqrtT = 10 / a;
142 final double t = sqrtT * sqrtT;
143 double z = DELTA[DELTA.length - 1];
144 for (int i = DELTA.length - 2; i >= 0; i--) {
145 z = t * z + DELTA[i];
146 }
147 return z / a + deltaMinusDeltaSum(a, b);
148 }
149
150 /**
151 * Returns the value of \( \log B(p, q) \) for \( 0 \leq x \leq 1 \) and \( p, q > 0 \).
152 * Based on the <em>NSWC Library of Mathematics Subroutines</em> implementation,
153 * {@code DBETLN}.
154 *
155 * @param p First argument.
156 * @param q Second argument.
157 * @return the value of \( \log B(p, q) \), or {@code NaN} if
158 * {@code p <= 0} or {@code q <= 0}.
159 */
160 public static double value(double p,
161 double q) {
162 if (Double.isNaN(p) ||
163 Double.isNaN(q) ||
164 p <= 0 ||
165 q <= 0) {
166 return Double.NaN;
167 }
168
169 final double a = Math.min(p, q);
170 final double b = Math.max(p, q);
171 if (a >= TEN) {
172 final double w = sumDeltaMinusDeltaSum(a, b);
173 final double h = a / b;
174 final double c = h / (1 + h);
175 final double u = -(a - 0.5) * Math.log(c);
176 final double v = b * Math.log1p(h);
177 if (u <= v) {
178 return (((-0.5 * Math.log(b) + HALF_LOG_TWO_PI) + w) - u) - v;
179 }
180 return (((-0.5 * Math.log(b) + HALF_LOG_TWO_PI) + w) - v) - u;
181 } else if (a > TWO) {
182 if (b > THOUSAND) {
183 final int n = (int) Math.floor(a - 1);
184 double prod = 1;
185 double ared = a;
186 for (int i = 0; i < n; i++) {
187 ared -= 1;
188 prod *= ared / (1 + ared / b);
189 }
190 return (Math.log(prod) - n * Math.log(b)) +
191 (LogGamma.value(ared) +
192 logGammaMinusLogGammaSum(ared, b));
193 }
194 double prod1 = 1;
195 double ared = a;
196 while (ared > 2) {
197 ared -= 1;
198 final double h = ared / b;
199 prod1 *= h / (1 + h);
200 }
201 if (b < TEN) {
202 double prod2 = 1;
203 double bred = b;
204 while (bred > 2) {
205 bred -= 1;
206 prod2 *= bred / (ared + bred);
207 }
208 return Math.log(prod1) +
209 Math.log(prod2) +
210 (LogGamma.value(ared) +
211 (LogGamma.value(bred) -
212 LogGammaSum.value(ared, bred)));
213 }
214 return Math.log(prod1) +
215 LogGamma.value(ared) +
216 logGammaMinusLogGammaSum(ared, b);
217 } else if (a >= 1) {
218 if (b > TWO) {
219 if (b < TEN) {
220 double prod = 1;
221 double bred = b;
222 while (bred > 2) {
223 bred -= 1;
224 prod *= bred / (a + bred);
225 }
226 return Math.log(prod) +
227 (LogGamma.value(a) +
228 (LogGamma.value(bred) -
229 LogGammaSum.value(a, bred)));
230 }
231 return LogGamma.value(a) +
232 logGammaMinusLogGammaSum(a, b);
233 }
234 return LogGamma.value(a) +
235 LogGamma.value(b) -
236 LogGammaSum.value(a, b);
237 } else {
238 if (b >= TEN) {
239 return LogGamma.value(a) +
240 logGammaMinusLogGammaSum(a, b);
241 }
242 // The original NSWC implementation was
243 // LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b));
244 // but the following command turned out to be more accurate.
245 // Note: Check for overflow that occurs if a and/or b are tiny.
246 final double beta = Gamma.value(a) * Gamma.value(b) / Gamma.value(a + b);
247 if (Double.isFinite(beta)) {
248 return Math.log(beta);
249 }
250 return LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b));
251 }
252 }
253
254 /**
255 * Returns the value of \( \log ( \Gamma(b) / \Gamma(a + b) ) \)
256 * for \( a \geq 0 \) and \( b \geq 10 \).
257 * Based on the <em>NSWC Library of Mathematics Subroutines</em> implementation,
258 * {@code DLGDIV}.
259 *
260 * @param a First argument.
261 * @param b Second argument.
262 * @return the value of \( \log(\Gamma(b) / \Gamma(a + b) \).
263 * @throws IllegalArgumentException if {@code a < 0} or {@code b < 10}.
264 */
265 private static double logGammaMinusLogGammaSum(double a,
266 double b) {
267 if (a < 0) {
268 throw new GammaException(GammaException.OUT_OF_RANGE, a, 0, Double.POSITIVE_INFINITY);
269 }
270 if (b < TEN) {
271 throw new GammaException(GammaException.OUT_OF_RANGE, b, TEN, Double.POSITIVE_INFINITY);
272 }
273
274 /*
275 * d = a + b - 0.5
276 */
277 final double d;
278 final double w;
279 if (a <= b) {
280 d = b + (a - 0.5);
281 w = deltaMinusDeltaSum(a, b);
282 } else {
283 d = a + (b - 0.5);
284 w = deltaMinusDeltaSum(b, a);
285 }
286
287 final double u = d * Math.log1p(a / b);
288 final double v = a * (Math.log(b) - 1);
289
290 return u <= v ?
291 (w - u) - v :
292 (w - v) - u;
293 }
294 }