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 // Assumptions:
87 // Never called with a < 0; a > b; or b < 10
88
89 final double h = a / b;
90 final double p = h / (1 + h);
91 final double q = 1 / (1 + h);
92 final double q2 = q * q;
93 /*
94 * s[i] = 1 + q + ... - q**(2 * i)
95 */
96 final double[] s = new double[DELTA.length];
97 s[0] = 1;
98 for (int i = 1; i < s.length; i++) {
99 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 > 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 }