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 if (a < 0 || 087 a > b) { 088 throw new GammaException(GammaException.OUT_OF_RANGE, a, 0, b); 089 } 090 if (b < TEN) { 091 throw new GammaException(GammaException.OUT_OF_RANGE, b, TEN, Double.POSITIVE_INFINITY); 092 } 093 094 final double h = a / b; 095 final double p = h / (1 + h); 096 final double q = 1 / (1 + h); 097 final double q2 = q * q; 098 /* 099 * 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}