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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.special;
23
24 import org.hipparchus.CalculusFieldElement;
25 import org.hipparchus.Field;
26 import org.hipparchus.exception.LocalizedCoreFormats;
27 import org.hipparchus.exception.MathIllegalArgumentException;
28 import org.hipparchus.exception.MathIllegalStateException;
29 import org.hipparchus.util.ContinuedFraction;
30 import org.hipparchus.util.FastMath;
31 import org.hipparchus.util.FieldContinuedFraction;
32
33 /**
34 * <p>
35 * This is a utility class that provides computation methods related to the
36 * Γ (Gamma) family of functions.
37 * </p>
38 * <p>
39 * Implementation of {@link #invGamma1pm1(double)} and
40 * {@link #logGamma1p(double)} is based on the algorithms described in
41 * </p>
42 * <ul>
43 * <li><a href="http://dx.doi.org/10.1145/22721.23109">Didonato and Morris
44 * (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios and
45 * their Inverse</em>, TOMS 12(4), 377-393,</li>
46 * <li><a href="http://dx.doi.org/10.1145/131766.131776">Didonato and Morris
47 * (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the
48 * Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li>
49 * </ul>
50 * <p>
51 * and implemented in the
52 * <a href="http://www.dtic.mil/docs/citations/ADA476840">NSWC Library of Mathematical Functions</a>,
53 * available
54 * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html">here</a>.
55 * This library is "approved for public release", and the
56 * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf">Copyright guidance</a>
57 * indicates that unless otherwise stated in the code, all FORTRAN functions in
58 * this library are license free. Since no such notice appears in the code these
59 * functions can safely be ported to Hipparchus.
60 * </p>
61 *
62 */
63 public class Gamma {
64 /**
65 * <a href="http://en.wikipedia.org/wiki/Euler-Mascheroni_constant">Euler-Mascheroni constant</a>
66 */
67 public static final double GAMMA = 0.577215664901532860606512090082; // NOPMD - the fact the function and the constant have the same name is intentional and comes from mathematics conventions
68
69 /**
70 * The value of the {@code g} constant in the Lanczos approximation, see
71 * {@link #lanczos(double)}.
72 */
73 public static final double LANCZOS_G = 607.0 / 128.0;
74
75 /** Maximum allowed numerical error. */
76 private static final double DEFAULT_EPSILON = 10e-15;
77
78 /** Lanczos coefficients */
79 private static final double[] LANCZOS = {
80 0.99999999999999709182,
81 57.156235665862923517,
82 -59.597960355475491248,
83 14.136097974741747174,
84 -0.49191381609762019978,
85 .33994649984811888699e-4,
86 .46523628927048575665e-4,
87 -.98374475304879564677e-4,
88 .15808870322491248884e-3,
89 -.21026444172410488319e-3,
90 .21743961811521264320e-3,
91 -.16431810653676389022e-3,
92 .84418223983852743293e-4,
93 -.26190838401581408670e-4,
94 .36899182659531622704e-5,
95 };
96
97 /** Avoid repeated computation of log of 2 PI in logGamma */
98 private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(2.0 * FastMath.PI);
99
100 /** The constant value of √(2π). */
101 private static final double SQRT_TWO_PI = 2.506628274631000502;
102
103 // limits for switching algorithm in digamma
104 /** C limit. */
105 private static final double C_LIMIT = 49;
106
107 /** S limit. */
108 private static final double S_LIMIT = 1e-8;
109
110 /*
111 * Constants for the computation of double invGamma1pm1(double).
112 * Copied from DGAM1 in the NSWC library.
113 */
114
115 /** The constant {@code A0} defined in {@code DGAM1}. */
116 private static final double INV_GAMMA1P_M1_A0 = .611609510448141581788E-08;
117
118 /** The constant {@code A1} defined in {@code DGAM1}. */
119 private static final double INV_GAMMA1P_M1_A1 = .624730830116465516210E-08;
120
121 /** The constant {@code B1} defined in {@code DGAM1}. */
122 private static final double INV_GAMMA1P_M1_B1 = .203610414066806987300E+00;
123
124 /** The constant {@code B2} defined in {@code DGAM1}. */
125 private static final double INV_GAMMA1P_M1_B2 = .266205348428949217746E-01;
126
127 /** The constant {@code B3} defined in {@code DGAM1}. */
128 private static final double INV_GAMMA1P_M1_B3 = .493944979382446875238E-03;
129
130 /** The constant {@code B4} defined in {@code DGAM1}. */
131 private static final double INV_GAMMA1P_M1_B4 = -.851419432440314906588E-05;
132
133 /** The constant {@code B5} defined in {@code DGAM1}. */
134 private static final double INV_GAMMA1P_M1_B5 = -.643045481779353022248E-05;
135
136 /** The constant {@code B6} defined in {@code DGAM1}. */
137 private static final double INV_GAMMA1P_M1_B6 = .992641840672773722196E-06;
138
139 /** The constant {@code B7} defined in {@code DGAM1}. */
140 private static final double INV_GAMMA1P_M1_B7 = -.607761895722825260739E-07;
141
142 /** The constant {@code B8} defined in {@code DGAM1}. */
143 private static final double INV_GAMMA1P_M1_B8 = .195755836614639731882E-09;
144
145 /** The constant {@code P0} defined in {@code DGAM1}. */
146 private static final double INV_GAMMA1P_M1_P0 = .6116095104481415817861E-08;
147
148 /** The constant {@code P1} defined in {@code DGAM1}. */
149 private static final double INV_GAMMA1P_M1_P1 = .6871674113067198736152E-08;
150
151 /** The constant {@code P2} defined in {@code DGAM1}. */
152 private static final double INV_GAMMA1P_M1_P2 = .6820161668496170657918E-09;
153
154 /** The constant {@code P3} defined in {@code DGAM1}. */
155 private static final double INV_GAMMA1P_M1_P3 = .4686843322948848031080E-10;
156
157 /** The constant {@code P4} defined in {@code DGAM1}. */
158 private static final double INV_GAMMA1P_M1_P4 = .1572833027710446286995E-11;
159
160 /** The constant {@code P5} defined in {@code DGAM1}. */
161 private static final double INV_GAMMA1P_M1_P5 = -.1249441572276366213222E-12;
162
163 /** The constant {@code P6} defined in {@code DGAM1}. */
164 private static final double INV_GAMMA1P_M1_P6 = .4343529937408594255178E-14;
165
166 /** The constant {@code Q1} defined in {@code DGAM1}. */
167 private static final double INV_GAMMA1P_M1_Q1 = .3056961078365221025009E+00;
168
169 /** The constant {@code Q2} defined in {@code DGAM1}. */
170 private static final double INV_GAMMA1P_M1_Q2 = .5464213086042296536016E-01;
171
172 /** The constant {@code Q3} defined in {@code DGAM1}. */
173 private static final double INV_GAMMA1P_M1_Q3 = .4956830093825887312020E-02;
174
175 /** The constant {@code Q4} defined in {@code DGAM1}. */
176 private static final double INV_GAMMA1P_M1_Q4 = .2692369466186361192876E-03;
177
178 /** The constant {@code C} defined in {@code DGAM1}. */
179 private static final double INV_GAMMA1P_M1_C = -.422784335098467139393487909917598E+00;
180
181 /** The constant {@code C0} defined in {@code DGAM1}. */
182 private static final double INV_GAMMA1P_M1_C0 = .577215664901532860606512090082402E+00;
183
184 /** The constant {@code C1} defined in {@code DGAM1}. */
185 private static final double INV_GAMMA1P_M1_C1 = -.655878071520253881077019515145390E+00;
186
187 /** The constant {@code C2} defined in {@code DGAM1}. */
188 private static final double INV_GAMMA1P_M1_C2 = -.420026350340952355290039348754298E-01;
189
190 /** The constant {@code C3} defined in {@code DGAM1}. */
191 private static final double INV_GAMMA1P_M1_C3 = .166538611382291489501700795102105E+00;
192
193 /** The constant {@code C4} defined in {@code DGAM1}. */
194 private static final double INV_GAMMA1P_M1_C4 = -.421977345555443367482083012891874E-01;
195
196 /** The constant {@code C5} defined in {@code DGAM1}. */
197 private static final double INV_GAMMA1P_M1_C5 = -.962197152787697356211492167234820E-02;
198
199 /** The constant {@code C6} defined in {@code DGAM1}. */
200 private static final double INV_GAMMA1P_M1_C6 = .721894324666309954239501034044657E-02;
201
202 /** The constant {@code C7} defined in {@code DGAM1}. */
203 private static final double INV_GAMMA1P_M1_C7 = -.116516759185906511211397108401839E-02;
204
205 /** The constant {@code C8} defined in {@code DGAM1}. */
206 private static final double INV_GAMMA1P_M1_C8 = -.215241674114950972815729963053648E-03;
207
208 /** The constant {@code C9} defined in {@code DGAM1}. */
209 private static final double INV_GAMMA1P_M1_C9 = .128050282388116186153198626328164E-03;
210
211 /** The constant {@code C10} defined in {@code DGAM1}. */
212 private static final double INV_GAMMA1P_M1_C10 = -.201348547807882386556893914210218E-04;
213
214 /** The constant {@code C11} defined in {@code DGAM1}. */
215 private static final double INV_GAMMA1P_M1_C11 = -.125049348214267065734535947383309E-05;
216
217 /** The constant {@code C12} defined in {@code DGAM1}. */
218 private static final double INV_GAMMA1P_M1_C12 = .113302723198169588237412962033074E-05;
219
220 /** The constant {@code C13} defined in {@code DGAM1}. */
221 private static final double INV_GAMMA1P_M1_C13 = -.205633841697760710345015413002057E-06;
222
223 /**
224 * Default constructor. Prohibit instantiation.
225 */
226 private Gamma() {}
227
228 /**
229 * <p>
230 * Returns the value of log Γ(x) for x > 0.
231 * </p>
232 * <p>
233 * For x ≤ 8, the implementation is based on the double precision
234 * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
235 * {@code DGAMLN}. For x > 8, the implementation is based on
236 * </p>
237 * <ul>
238 * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
239 * Function</a>, equation (28).</li>
240 * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
241 * Lanczos Approximation</a>, equations (1) through (5).</li>
242 * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
243 * the computation of the convergent Lanczos complex Gamma
244 * approximation</a></li>
245 * </ul>
246 *
247 * @param x Argument.
248 * @return the value of {@code log(Gamma(x))}, {@code Double.NaN} if
249 * {@code x <= 0.0}.
250 */
251 public static double logGamma(double x) {
252 double ret;
253
254 if (Double.isNaN(x) || (x <= 0.0)) {
255 ret = Double.NaN;
256 } else if (x < 0.5) {
257 return logGamma1p(x) - FastMath.log(x);
258 } else if (x <= 2.5) {
259 return logGamma1p((x - 0.5) - 0.5);
260 } else if (x <= 8.0) {
261 final int n = (int) FastMath.floor(x - 1.5);
262 double prod = 1.0;
263 for (int i = 1; i <= n; i++) {
264 prod *= x - i;
265 }
266 return logGamma1p(x - (n + 1)) + FastMath.log(prod);
267 } else {
268 double sum = lanczos(x);
269 double tmp = x + LANCZOS_G + .5;
270 ret = ((x + .5) * FastMath.log(tmp)) - tmp +
271 HALF_LOG_2_PI + FastMath.log(sum / x);
272 }
273
274 return ret;
275 }
276
277 /**
278 * <p>
279 * Returns the value of log Γ(x) for x > 0.
280 * </p>
281 * <p>
282 * For x ≤ 8, the implementation is based on the double precision
283 * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
284 * {@code DGAMLN}. For x > 8, the implementation is based on
285 * </p>
286 * <ul>
287 * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
288 * Function</a>, equation (28).</li>
289 * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
290 * Lanczos Approximation</a>, equations (1) through (5).</li>
291 * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
292 * the computation of the convergent Lanczos complex Gamma
293 * approximation</a></li>
294 * </ul>
295 *
296 * @param x Argument.
297 * @param <T> Type of the field elements.
298 * @return the value of {@code log(Gamma(x))}, {@code Double.NaN} if
299 * {@code x <= 0.0}.
300 */
301 public static <T extends CalculusFieldElement<T>> T logGamma(T x) {
302 final Field<T> field = x.getField();
303 T ret;
304
305 if (x.isNaN() || (x.getReal() <= 0.0)) {
306 ret = field.getOne().multiply(Double.NaN);
307 }
308 else if (x.getReal() < 0.5) {
309 return logGamma1p(x).subtract(x.log());
310 }
311 else if (x.getReal() <= 2.5) {
312 return logGamma1p(x.subtract(1));
313 }
314 else if (x.getReal() <= 8.0) {
315 final int n = (int) x.subtract(1.5).floor().getReal();
316 T prod = field.getOne();
317 for (int i = 1; i <= n; i++) {
318 prod = prod.multiply(x.subtract(i));
319 }
320 return logGamma1p(x.subtract(n + 1)).add(prod.log());
321 }
322 else {
323 T sum = lanczos(x);
324 T tmp = x.add(LANCZOS_G + .5);
325 ret = x.add(.5).multiply(tmp.log()).subtract(tmp).add(HALF_LOG_2_PI).add(sum.divide(x).log());
326 }
327
328 return ret;
329 }
330
331 /**
332 * Returns the regularized gamma function P(a, x).
333 *
334 * @param a Parameter.
335 * @param x Value.
336 * @return the regularized gamma function P(a, x).
337 * @throws MathIllegalStateException if the algorithm fails to converge.
338 */
339 public static double regularizedGammaP(double a, double x) {
340 return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
341 }
342
343 /**
344 * Returns the regularized gamma function P(a, x).
345 *
346 * @param a Parameter.
347 * @param x Value.
348 * @param <T> Type of the field elements.
349 * @return the regularized gamma function P(a, x).
350 * @throws MathIllegalStateException if the algorithm fails to converge.
351 */
352 public static <T extends CalculusFieldElement<T>> T regularizedGammaP(T a, T x) {
353 return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
354 }
355
356 /**
357 * Returns the regularized gamma function P(a, x).
358 * <p>
359 * The implementation of this method is based on:
360 * <ul>
361 * <li>
362 * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
363 * Regularized Gamma Function</a>, equation (1)
364 * </li>
365 * <li>
366 * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
367 * Incomplete Gamma Function</a>, equation (4).
368 * </li>
369 * <li>
370 * <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html">
371 * Confluent Hypergeometric Function of the First Kind</a>, equation (1).
372 * </li>
373 * </ul>
374 *
375 * @param a the a parameter.
376 * @param x the value.
377 * @param epsilon When the absolute value of the nth item in the
378 * series is less than epsilon the approximation ceases to calculate
379 * further elements in the series.
380 * @param maxIterations Maximum number of "iterations" to complete.
381 * @return the regularized gamma function P(a, x)
382 * @throws MathIllegalStateException if the algorithm fails to converge.
383 */
384 public static double regularizedGammaP(double a,
385 double x,
386 double epsilon,
387 int maxIterations) {
388 double ret;
389
390 if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
391 ret = Double.NaN;
392 } else if (x == 0.0) {
393 ret = 0.0;
394 } else if (x >= a + 1) {
395 // use regularizedGammaQ because it should converge faster in this
396 // case.
397 ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations);
398 } else {
399 // calculate series
400 double n = 0.0; // current element index
401 double an = 1.0 / a; // n-th element in the series
402 double sum = an; // partial sum
403 while (FastMath.abs(an/sum) > epsilon &&
404 n < maxIterations &&
405 sum < Double.POSITIVE_INFINITY) {
406 // compute next element in the series
407 n += 1.0;
408 an *= x / (a + n);
409
410 // update partial sum
411 sum += an;
412 }
413 if (n >= maxIterations) {
414 throw new MathIllegalStateException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, maxIterations);
415 } else if (Double.isInfinite(sum)) {
416 ret = 1.0;
417 } else {
418 ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * sum;
419 }
420 }
421
422 return ret;
423 }
424
425 /**
426 * Returns the regularized gamma function P(a, x).
427 * <p>
428 * The implementation of this method is based on:
429 * <ul>
430 * <li>
431 * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
432 * Regularized Gamma Function</a>, equation (1)
433 * </li>
434 * <li>
435 * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
436 * Incomplete Gamma Function</a>, equation (4).
437 * </li>
438 * <li>
439 * <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html">
440 * Confluent Hypergeometric Function of the First Kind</a>, equation (1).
441 * </li>
442 * </ul>
443 *
444 * @param a the a parameter.
445 * @param x the value.
446 * @param epsilon When the absolute value of the nth item in the
447 * series is less than epsilon the approximation ceases to calculate
448 * further elements in the series.
449 * @param maxIterations Maximum number of "iterations" to complete.
450 * @param <T> Type of the field elements.
451 * @return the regularized gamma function P(a, x)
452 * @throws MathIllegalStateException if the algorithm fails to converge.
453 */
454 public static <T extends CalculusFieldElement<T>> T regularizedGammaP(T a,
455 T x,
456 double epsilon,
457 int maxIterations) {
458 final Field<T> field = x.getField();
459 final T zero = field.getZero();
460 final T one = field.getOne();
461
462 T ret;
463
464 if (a.isNaN() || x.isNaN() || (a.getReal() <= 0.0) || (x.getReal() < 0.0)) {
465 ret = one.multiply(Double.NaN);
466 }
467 else if (x.getReal() == 0.0) {
468 ret = zero;
469 }
470 else if (x.getReal() >= a.add(1).getReal()) {
471 // use regularizedGammaQ because it should converge faster in this
472 // case.
473 ret = one.subtract(regularizedGammaQ(a, x, epsilon, maxIterations));
474 }
475 else {
476 // calculate series
477 double n = 0.0; // current element index
478 T an = one.divide(a); // n-th element in the series
479 T sum = an; // partial sum
480 while (an.divide(sum).abs().getReal() > epsilon &&
481 n < maxIterations &&
482 sum.getReal() < Double.POSITIVE_INFINITY) {
483 // compute next element in the series
484 n += 1.0;
485 an = an.multiply(x.divide(a.add(n)));
486
487 // update partial sum
488 sum = sum.add(an);
489 }
490 if (n >= maxIterations) {
491 throw new MathIllegalStateException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, maxIterations);
492 }
493 else if (sum.isInfinite()) {
494 ret = one;
495 }
496 else {
497 ret = a.multiply(x.log()).subtract(logGamma(a)).subtract(x).exp().multiply(sum);
498 }
499 }
500
501 return ret;
502 }
503
504 /**
505 * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
506 *
507 * @param a the a parameter.
508 * @param x the value.
509 * @return the regularized gamma function Q(a, x)
510 * @throws MathIllegalStateException if the algorithm fails to converge.
511 */
512 public static double regularizedGammaQ(double a, double x) {
513 return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
514 }
515
516 /**
517 * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
518 *
519 * @param a the a parameter.
520 * @param x the value.
521 * @param <T> Type of the field elements.
522 * @return the regularized gamma function Q(a, x)
523 * @throws MathIllegalStateException if the algorithm fails to converge.
524 */
525 public static <T extends CalculusFieldElement<T>> T regularizedGammaQ(T a, T x) {
526 return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
527 }
528
529 /**
530 * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
531 * <p>
532 * The implementation of this method is based on:
533 * <ul>
534 * <li>
535 * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
536 * Regularized Gamma Function</a>, equation (1).
537 * </li>
538 * <li>
539 * <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/">
540 * Regularized incomplete gamma function: Continued fraction representations
541 * (formula 06.08.10.0003)</a>
542 * </li>
543 * </ul>
544 *
545 * @param a the a parameter.
546 * @param x the value.
547 * @param epsilon When the absolute value of the nth item in the
548 * series is less than epsilon the approximation ceases to calculate
549 * further elements in the series.
550 * @param maxIterations Maximum number of "iterations" to complete.
551 * @return the regularized gamma function P(a, x)
552 * @throws MathIllegalStateException if the algorithm fails to converge.
553 */
554 public static double regularizedGammaQ(final double a,
555 double x,
556 double epsilon,
557 int maxIterations) {
558 double ret;
559
560 if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
561 ret = Double.NaN;
562 } else if (x == 0.0) {
563 ret = 1.0;
564 } else if (x < a + 1.0) {
565 // use regularizedGammaP because it should converge faster in this
566 // case.
567 ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations);
568 } else {
569 // create continued fraction
570 ContinuedFraction cf = new ContinuedFraction() {
571
572 /** {@inheritDoc} */
573 @Override
574 protected double getA(int n, double x) {
575 return ((2.0 * n) + 1.0) - a + x;
576 }
577
578 /** {@inheritDoc} */
579 @Override
580 protected double getB(int n, double x) {
581 return n * (a - n);
582 }
583 };
584
585 ret = 1.0 / cf.evaluate(x, epsilon, maxIterations);
586 ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * ret;
587 }
588
589 return ret;
590 }
591
592 /**
593 * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
594 * <p>
595 * The implementation of this method is based on:
596 * <ul>
597 * <li>
598 * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
599 * Regularized Gamma Function</a>, equation (1).
600 * </li>
601 * <li>
602 * <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/">
603 * Regularized incomplete gamma function: Continued fraction representations
604 * (formula 06.08.10.0003)</a>
605 * </li>
606 * </ul>
607 *
608 * @param a the a parameter.
609 * @param x the value.
610 * @param epsilon When the absolute value of the nth item in the
611 * series is less than epsilon the approximation ceases to calculate
612 * further elements in the series.
613 * @param maxIterations Maximum number of "iterations" to complete.
614 * @param <T> Type fo the field elements.
615 * @return the regularized gamma function P(a, x)
616 * @throws MathIllegalStateException if the algorithm fails to converge.
617 */
618 public static <T extends CalculusFieldElement<T>> T regularizedGammaQ(final T a,
619 T x,
620 double epsilon,
621 int maxIterations) {
622 final Field<T> field = x.getField();
623 final T one = field.getOne();
624
625 T ret;
626
627 if (a.isNaN() || x.isNaN() || a.getReal() <= 0.0 || x.getReal() < 0.0) {
628 ret = field.getOne().multiply(Double.NaN);
629 }
630 else if (x.getReal() == 0.0) {
631 ret = one;
632 }
633 else if (x.getReal() < a.add(1.0).getReal()) {
634 // use regularizedGammaP because it should converge faster in this
635 // case.
636 ret = one.subtract(regularizedGammaP(a, x, epsilon, maxIterations));
637 }
638 else {
639 // create continued fraction
640 FieldContinuedFraction cf = new FieldContinuedFraction() {
641
642 /** {@inheritDoc} */
643 @Override
644 @SuppressWarnings("unchecked")
645 public <C extends CalculusFieldElement<C>> C getA(final int n, final C x) {
646 return x.subtract((C) a).add((2.0 * n) + 1.0);
647 }
648
649 /** {@inheritDoc} */
650 @Override
651 @SuppressWarnings("unchecked")
652 public <C extends CalculusFieldElement<C>> C getB(final int n, final C x) {
653 return (C) a.subtract(n).multiply(n);
654 }
655 };
656
657 ret = one.divide(cf.evaluate(x, epsilon, maxIterations));
658 ret = a.multiply(x.log()).subtract(logGamma(a)).subtract(x).exp().multiply(ret);
659 }
660
661 return ret;
662 }
663
664 /**
665 * <p>Computes the digamma function of x.</p>
666 *
667 * <p>This is an independently written implementation of the algorithm described in
668 * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p>
669 *
670 * <p>Some of the constants have been changed to increase accuracy at the moderate expense
671 * of run-time. The result should be accurate to within 10^-8 absolute tolerance for
672 * x >= 10^-5 and within 10^-8 relative tolerance for x > 0.</p>
673 *
674 * <p>Performance for large negative values of x will be quite expensive (proportional to
675 * |x|). Accuracy for negative values of x should be about 10^-8 absolute for results
676 * less than 10^5 and 10^-8 relative for results larger than that.</p>
677 *
678 * @param x Argument.
679 * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller.
680 * @see <a href="http://en.wikipedia.org/wiki/Digamma_function">Digamma</a>
681 * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf">Bernardo's original article </a>
682 */
683 public static double digamma(double x) {
684 if (Double.isNaN(x) || Double.isInfinite(x)) {
685 return x;
686 }
687
688 if (x > 0 && x <= S_LIMIT) {
689 // use method 5 from Bernardo AS103
690 // accurate to O(x)
691 return -GAMMA - 1 / x;
692 }
693 if (x >= C_LIMIT) {
694 // use method 8 (accurate to O(1/x^8))
695 double inv = 1 / (x * x);
696 // 1 1 1 1 1 5 691 1
697 // log(x) - --- - ------ + ------- - ------- + ------- - ------- + ---------- - -------
698 // 2 x 12 x^2 120 x^4 252 x^6 240 x^8 660 x^10 32760 x^12 12 x^14
699 return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv * (1.0 / 252 + inv *
700 (1.0 / 240 - inv * (5.0 / 660 + inv * (691.0 / 32760 - inv / 12))))));
701 }
702
703 return digamma(x + 1) - 1 / x;
704 }
705
706 /**
707 * <p>Computes the digamma function of x.</p>
708 *
709 * <p>This is an independently written implementation of the algorithm described in
710 * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p>
711 *
712 * <p>Some of the constants have been changed to increase accuracy at the moderate expense
713 * of run-time. The result should be accurate to within 10^-8 absolute tolerance for
714 * x >= 10^-5 and within 10^-8 relative tolerance for x > 0.</p>
715 *
716 * <p>Performance for large negative values of x will be quite expensive (proportional to
717 * |x|). Accuracy for negative values of x should be about 10^-8 absolute for results
718 * less than 10^5 and 10^-8 relative for results larger than that.</p>
719 *
720 * @param x Argument.
721 * @param <T> Type of the field elements.
722 * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller.
723 * @see <a href="http://en.wikipedia.org/wiki/Digamma_function">Digamma</a>
724 * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf">Bernardo's original article </a>
725 */
726 public static <T extends CalculusFieldElement<T>> T digamma(T x) {
727 if (x.isNaN() || x.isInfinite()) {
728 return x;
729 }
730
731 if (x.getReal() > 0 && x.getReal() <= S_LIMIT) {
732 // use method 5 from Bernardo AS103
733 // accurate to O(x)
734 return x.pow(-1).negate().subtract(GAMMA);
735 }
736
737 if (x.getReal() >= C_LIMIT) {
738 // use method 8 (accurate to O(1/x^8))
739 T inv = x.square().reciprocal();
740 // 1 1 1 1 1 5 691 1
741 // log(x) - --- - ------ + ------- - ------- + ------- - ------- + ---------- - -------
742 // 2 x 12 x^2 120 x^4 252 x^6 240 x^8 660 x^10 32760 x^12 12 x^14
743 return x.log().subtract(x.pow(-1).multiply(0.5)).add(
744 inv.multiply(
745 inv.multiply(
746 inv.multiply(
747 inv.multiply(
748 inv.multiply(
749 inv.multiply(inv.divide(-12.)
750 .add(691. / 32760))
751 .subtract(5. / 660))
752 .add(1.0 / 240))
753 .subtract(1.0 / 252))
754 .add(1.0 / 120))
755 .subtract(1.0 / 12)));
756 }
757
758 return digamma(x.add(1.)).subtract(x.pow(-1));
759 }
760
761 /**
762 * Computes the trigamma function of x.
763 * This function is derived by taking the derivative of the implementation
764 * of digamma.
765 *
766 * @param x Argument.
767 * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller
768 * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function">Trigamma</a>
769 * @see Gamma#digamma(double)
770 */
771 public static double trigamma(double x) {
772 if (Double.isNaN(x) || Double.isInfinite(x)) {
773 return x;
774 }
775
776 if (x > 0 && x <= S_LIMIT) {
777 return 1 / (x * x);
778 }
779
780 if (x >= C_LIMIT) {
781 double inv = 1 / (x * x);
782 // 1 1 1 1 1 1 5 691 7
783 // - + ---- + ---- - ----- + ----- - ----- + ------- - -------- + ------
784 // x 2 3 5 7 9 11 13 15
785 // 2 x 6 x 30 x 42 x 30 x 66 x 2730 x 6 x
786 return 1 / x + inv * 0.5 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv * (1.0 / 42 - inv * (1.0 / 30 + inv *
787 (5.0 / 66 - inv * (691. / 2730 + inv * 7. / 15))))));
788 }
789
790 return trigamma(x + 1) + 1 / (x * x);
791 }
792
793 /**
794 * Computes the trigamma function of x.
795 * This function is derived by taking the derivative of the implementation
796 * of digamma.
797 *
798 * @param x Argument.
799 * @param <T> Type of the field elements.
800 * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller
801 * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function">Trigamma</a>
802 * @see Gamma#digamma(double)
803 */
804 public static <T extends CalculusFieldElement<T>> T trigamma(T x) {
805 if (x.isNaN() || x.isInfinite()) {
806 return x;
807 }
808
809 if (x.getReal() > 0 && x.getReal() <= S_LIMIT) {
810 // use method 5 from Bernardo AS103
811 // accurate to O(x)
812 return x.square().reciprocal();
813 }
814
815 if (x.getReal() >= C_LIMIT) {
816 // use method 4 (accurate to O(1/x^8)
817 T inv = x.square().reciprocal();
818 T invCub = inv.multiply(x.reciprocal());
819 // 1 1 1 1 1 1 5 691 7
820 // - + ---- + ---- - ----- + ----- + ----- + ------- - -------- + ------
821 // x 2 3 5 7 9 11 13 15
822 // 2 x 6 x 30 x 42 x 30 x 66 x 2730 x 6 x
823 return x.pow(-1).add(
824 inv.multiply(0.5)).add(
825 invCub.multiply(
826 inv.multiply(
827 inv.multiply(
828 inv.multiply(
829 inv.multiply(
830 inv.multiply(inv.multiply(7. / 6)
831 .subtract(691. / 2730))
832 .add(5. / 66))
833 .subtract(1.0 / 30))
834 .add(1.0 / 42))
835 .subtract(1.0 / 30))
836 .add(1.0 / 6)));
837 }
838
839 return trigamma(x.add(1.)).add(x.square().reciprocal());
840 }
841
842 /**
843 * <p>
844 * Returns the Lanczos approximation used to compute the gamma function.
845 * The Lanczos approximation is related to the Gamma function by the
846 * following equation
847 * \[
848 * \Gamma(x) = \frac{\sqrt{2\pi}}{x} \times (x + g + \frac{1}{2}) ^ (x + \frac{1}{2})
849 * \times e^{-x - g - 0.5} \times \mathrm{lanczos}(x)
850 * \]
851 * where {@code g} is the Lanczos constant.
852 * </p>
853 *
854 * @param x Argument.
855 * @return The Lanczos approximation.
856 * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html">Lanczos Approximation</a>
857 * equations (1) through (5), and Paul Godfrey's
858 * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation
859 * of the convergent Lanczos complex Gamma approximation</a>
860 */
861 public static double lanczos(final double x) {
862 double sum = 0.0;
863 for (int i = LANCZOS.length - 1; i > 0; --i) {
864 sum += LANCZOS[i] / (x + i);
865 }
866 return sum + LANCZOS[0];
867 }
868
869 /**
870 * <p>
871 * Returns the Lanczos approximation used to compute the gamma function.
872 * The Lanczos approximation is related to the Gamma function by the
873 * following equation
874 * \[
875 * \Gamma(x) = \frac{\sqrt{2\pi}}{x} \times (x + g + \frac{1}{2}) ^ (x + \frac{1}{2})
876 * \times e^{-x - g - 0.5} \times \mathrm{lanczos}(x)
877 * \]
878 * where {@code g} is the Lanczos constant.
879 * </p>
880 *
881 * @param x Argument.
882 * @param <T> Type of the field elements.
883 * @return The Lanczos approximation.
884 * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html">Lanczos Approximation</a>
885 * equations (1) through (5), and Paul Godfrey's
886 * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation
887 * of the convergent Lanczos complex Gamma approximation</a>
888 */
889 public static <T extends CalculusFieldElement<T>> T lanczos(final T x) {
890 final Field<T> field = x.getField();
891 T sum = field.getZero();
892 for (int i = LANCZOS.length - 1; i > 0; --i) {
893 sum = sum.add(x.add(i).pow(-1.).multiply(LANCZOS[i]));
894 }
895 return sum.add(LANCZOS[0]);
896 }
897
898 /**
899 * Returns the value of 1 / Γ(1 + x) - 1 for -0.5 ≤ x ≤
900 * 1.5. This implementation is based on the double precision
901 * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
902 * {@code DGAM1}.
903 *
904 * @param x Argument.
905 * @return The value of {@code 1.0 / Gamma(1.0 + x) - 1.0}.
906 * @throws MathIllegalArgumentException if {@code x < -0.5}
907 * @throws MathIllegalArgumentException if {@code x > 1.5}
908 */
909 public static double invGamma1pm1(final double x) {
910
911 if (x < -0.5) {
912 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
913 x, -0.5);
914 }
915 if (x > 1.5) {
916 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
917 x, 1.5);
918 }
919
920 final double ret;
921 final double t = x <= 0.5 ? x : (x - 0.5) - 0.5;
922 if (t < 0.0) {
923 final double a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1;
924 double b = INV_GAMMA1P_M1_B8;
925 b = INV_GAMMA1P_M1_B7 + t * b;
926 b = INV_GAMMA1P_M1_B6 + t * b;
927 b = INV_GAMMA1P_M1_B5 + t * b;
928 b = INV_GAMMA1P_M1_B4 + t * b;
929 b = INV_GAMMA1P_M1_B3 + t * b;
930 b = INV_GAMMA1P_M1_B2 + t * b;
931 b = INV_GAMMA1P_M1_B1 + t * b;
932 b = 1.0 + t * b;
933
934 double c = INV_GAMMA1P_M1_C13 + t * (a / b);
935 c = INV_GAMMA1P_M1_C12 + t * c;
936 c = INV_GAMMA1P_M1_C11 + t * c;
937 c = INV_GAMMA1P_M1_C10 + t * c;
938 c = INV_GAMMA1P_M1_C9 + t * c;
939 c = INV_GAMMA1P_M1_C8 + t * c;
940 c = INV_GAMMA1P_M1_C7 + t * c;
941 c = INV_GAMMA1P_M1_C6 + t * c;
942 c = INV_GAMMA1P_M1_C5 + t * c;
943 c = INV_GAMMA1P_M1_C4 + t * c;
944 c = INV_GAMMA1P_M1_C3 + t * c;
945 c = INV_GAMMA1P_M1_C2 + t * c;
946 c = INV_GAMMA1P_M1_C1 + t * c;
947 c = INV_GAMMA1P_M1_C + t * c;
948 if (x > 0.5) {
949 ret = t * c / x;
950 } else {
951 ret = x * ((c + 0.5) + 0.5);
952 }
953 } else {
954 double p = INV_GAMMA1P_M1_P6;
955 p = INV_GAMMA1P_M1_P5 + t * p;
956 p = INV_GAMMA1P_M1_P4 + t * p;
957 p = INV_GAMMA1P_M1_P3 + t * p;
958 p = INV_GAMMA1P_M1_P2 + t * p;
959 p = INV_GAMMA1P_M1_P1 + t * p;
960 p = INV_GAMMA1P_M1_P0 + t * p;
961
962 double q = INV_GAMMA1P_M1_Q4;
963 q = INV_GAMMA1P_M1_Q3 + t * q;
964 q = INV_GAMMA1P_M1_Q2 + t * q;
965 q = INV_GAMMA1P_M1_Q1 + t * q;
966 q = 1.0 + t * q;
967
968 double c = INV_GAMMA1P_M1_C13 + (p / q) * t;
969 c = INV_GAMMA1P_M1_C12 + t * c;
970 c = INV_GAMMA1P_M1_C11 + t * c;
971 c = INV_GAMMA1P_M1_C10 + t * c;
972 c = INV_GAMMA1P_M1_C9 + t * c;
973 c = INV_GAMMA1P_M1_C8 + t * c;
974 c = INV_GAMMA1P_M1_C7 + t * c;
975 c = INV_GAMMA1P_M1_C6 + t * c;
976 c = INV_GAMMA1P_M1_C5 + t * c;
977 c = INV_GAMMA1P_M1_C4 + t * c;
978 c = INV_GAMMA1P_M1_C3 + t * c;
979 c = INV_GAMMA1P_M1_C2 + t * c;
980 c = INV_GAMMA1P_M1_C1 + t * c;
981 c = INV_GAMMA1P_M1_C0 + t * c;
982
983 if (x > 0.5) {
984 ret = (t / x) * ((c - 0.5) - 0.5);
985 } else {
986 ret = x * c;
987 }
988 }
989
990 return ret;
991 }
992
993 /**
994 * Returns the value of 1 / Γ(1 + x) - 1 for -0.5 ≤ x ≤
995 * 1.5. This implementation is based on the double precision
996 * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
997 * {@code DGAM1}.
998 *
999 * @param x Argument.
1000 * @param <T> Type of the field elements.
1001 * @return The value of {@code 1.0 / Gamma(1.0 + x) - 1.0}.
1002 * @throws MathIllegalArgumentException if {@code x < -0.5}
1003 * @throws MathIllegalArgumentException if {@code x > 1.5}
1004 */
1005 public static <T extends CalculusFieldElement<T>> T invGamma1pm1(final T x) {
1006 final T one = x.getField().getOne();
1007
1008 if (x.getReal() < -0.5) {
1009 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
1010 x, -0.5);
1011 }
1012 if (x.getReal() > 1.5) {
1013 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
1014 x, 1.5);
1015 }
1016
1017 final T ret;
1018 final T t = x.getReal() <= 0.5 ? x : x.subtract(1);
1019 if (t.getReal() < 0.0) {
1020 final T a = one.newInstance(INV_GAMMA1P_M1_A0).add(t.multiply(INV_GAMMA1P_M1_A1));
1021 T b = one.newInstance(INV_GAMMA1P_M1_B8);
1022 b = t.multiply(b).add(INV_GAMMA1P_M1_B7);
1023 b = t.multiply(b).add(INV_GAMMA1P_M1_B6);
1024 b = t.multiply(b).add(INV_GAMMA1P_M1_B5);
1025 b = t.multiply(b).add(INV_GAMMA1P_M1_B4);
1026 b = t.multiply(b).add(INV_GAMMA1P_M1_B3);
1027 b = t.multiply(b).add(INV_GAMMA1P_M1_B2);
1028 b = t.multiply(b).add(INV_GAMMA1P_M1_B1);
1029 b = t.multiply(b).add(1.);
1030
1031 T c = one.newInstance(INV_GAMMA1P_M1_C13).add(t.multiply(a.divide(b)));
1032 c = t.multiply(c).add(INV_GAMMA1P_M1_C12);
1033 c = t.multiply(c).add(INV_GAMMA1P_M1_C11);
1034 c = t.multiply(c).add(INV_GAMMA1P_M1_C10);
1035 c = t.multiply(c).add(INV_GAMMA1P_M1_C9);
1036 c = t.multiply(c).add(INV_GAMMA1P_M1_C8);
1037 c = t.multiply(c).add(INV_GAMMA1P_M1_C7);
1038 c = t.multiply(c).add(INV_GAMMA1P_M1_C6);
1039 c = t.multiply(c).add(INV_GAMMA1P_M1_C5);
1040 c = t.multiply(c).add(INV_GAMMA1P_M1_C4);
1041 c = t.multiply(c).add(INV_GAMMA1P_M1_C3);
1042 c = t.multiply(c).add(INV_GAMMA1P_M1_C2);
1043 c = t.multiply(c).add(INV_GAMMA1P_M1_C1);
1044 c = t.multiply(c).add(INV_GAMMA1P_M1_C);
1045 if (x.getReal() > 0.5) {
1046 ret = t.multiply(c).divide(x);
1047 }
1048 else {
1049 ret = x.multiply(c.add(1));
1050 }
1051 }
1052 else {
1053 T p = one.newInstance(INV_GAMMA1P_M1_P6);
1054 p = t.multiply(p).add(INV_GAMMA1P_M1_P5);
1055 p = t.multiply(p).add(INV_GAMMA1P_M1_P4);
1056 p = t.multiply(p).add(INV_GAMMA1P_M1_P3);
1057 p = t.multiply(p).add(INV_GAMMA1P_M1_P2);
1058 p = t.multiply(p).add(INV_GAMMA1P_M1_P1);
1059 p = t.multiply(p).add(INV_GAMMA1P_M1_P0);
1060
1061 T q = one.newInstance(INV_GAMMA1P_M1_Q4);
1062 q = t.multiply(q).add(INV_GAMMA1P_M1_Q3);
1063 q = t.multiply(q).add(INV_GAMMA1P_M1_Q2);
1064 q = t.multiply(q).add(INV_GAMMA1P_M1_Q1);
1065 q = t.multiply(q).add(1.);
1066
1067 T c = one.newInstance(INV_GAMMA1P_M1_C13).add(t.multiply(p.divide(q)));
1068 c = t.multiply(c).add(INV_GAMMA1P_M1_C12);
1069 c = t.multiply(c).add(INV_GAMMA1P_M1_C11);
1070 c = t.multiply(c).add(INV_GAMMA1P_M1_C10);
1071 c = t.multiply(c).add(INV_GAMMA1P_M1_C9);
1072 c = t.multiply(c).add(INV_GAMMA1P_M1_C8);
1073 c = t.multiply(c).add(INV_GAMMA1P_M1_C7);
1074 c = t.multiply(c).add(INV_GAMMA1P_M1_C6);
1075 c = t.multiply(c).add(INV_GAMMA1P_M1_C5);
1076 c = t.multiply(c).add(INV_GAMMA1P_M1_C4);
1077 c = t.multiply(c).add(INV_GAMMA1P_M1_C3);
1078 c = t.multiply(c).add(INV_GAMMA1P_M1_C2);
1079 c = t.multiply(c).add(INV_GAMMA1P_M1_C1);
1080 c = t.multiply(c).add(INV_GAMMA1P_M1_C0);
1081
1082 if (x.getReal() > 0.5) {
1083 ret = t.divide(x).multiply(c.subtract(1));
1084 }
1085 else {
1086 ret = x.multiply(c);
1087 }
1088 }
1089
1090 return ret;
1091 }
1092
1093 /**
1094 * Returns the value of log Γ(1 + x) for -0.5 ≤ x ≤ 1.5.
1095 * This implementation is based on the double precision implementation in
1096 * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}.
1097 *
1098 * @param x Argument.
1099 * @return The value of {@code log(Gamma(1 + x))}.
1100 * @throws MathIllegalArgumentException if {@code x < -0.5}.
1101 * @throws MathIllegalArgumentException if {@code x > 1.5}.
1102 */
1103 public static double logGamma1p(final double x)
1104 throws MathIllegalArgumentException {
1105
1106 if (x < -0.5) {
1107 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
1108 x, -0.5);
1109 }
1110 if (x > 1.5) {
1111 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
1112 x, 1.5);
1113 }
1114
1115 return -FastMath.log1p(invGamma1pm1(x));
1116 }
1117
1118 /**
1119 * Returns the value of log Γ(1 + x) for -0.5 ≤ x ≤ 1.5.
1120 * This implementation is based on the double precision implementation in
1121 * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}.
1122 *
1123 * @param x Argument.
1124 * @param <T> Type of the field elements.
1125 * @return The value of {@code log(Gamma(1 + x))}.
1126 * @throws MathIllegalArgumentException if {@code x < -0.5}.
1127 * @throws MathIllegalArgumentException if {@code x > 1.5}.
1128 */
1129 public static <T extends CalculusFieldElement<T>> T logGamma1p(final T x)
1130 throws MathIllegalArgumentException {
1131
1132 if (x.getReal() < -0.5) {
1133 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
1134 x, -0.5);
1135 }
1136 if (x.getReal() > 1.5) {
1137 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
1138 x, 1.5);
1139 }
1140
1141 return invGamma1pm1(x).log1p().negate();
1142 }
1143
1144
1145 /**
1146 * Returns the value of Γ(x). Based on the <em>NSWC Library of
1147 * Mathematics Subroutines</em> double precision implementation,
1148 * {@code DGAMMA}.
1149 *
1150 * @param x Argument.
1151 * @return the value of {@code Gamma(x)}.
1152 */
1153 public static double gamma(final double x) {
1154
1155 if ((x == FastMath.rint(x)) && (x <= 0.0)) {
1156 return Double.NaN;
1157 }
1158
1159 final double ret;
1160 final double absX = FastMath.abs(x);
1161 if (absX <= 20.0) {
1162 if (x >= 1.0) {
1163 /*
1164 * From the recurrence relation
1165 * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
1166 * then
1167 * Gamma(t) = 1 / [1 + invGamma1pm1(t - 1)],
1168 * where t = x - n. This means that t must satisfy
1169 * -0.5 <= t - 1 <= 1.5.
1170 */
1171 double prod = 1.0;
1172 double t = x;
1173 while (t > 2.5) {
1174 t -= 1.0;
1175 prod *= t;
1176 }
1177 ret = prod / (1.0 + invGamma1pm1(t - 1.0));
1178 } else {
1179 /*
1180 * From the recurrence relation
1181 * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
1182 * then
1183 * Gamma(x + n + 1) = 1 / [1 + invGamma1pm1(x + n)],
1184 * which requires -0.5 <= x + n <= 1.5.
1185 */
1186 double prod = x;
1187 double t = x;
1188 while (t < -0.5) {
1189 t += 1.0;
1190 prod *= t;
1191 }
1192 ret = 1.0 / (prod * (1.0 + invGamma1pm1(t)));
1193 }
1194 } else {
1195 final double y = absX + LANCZOS_G + 0.5;
1196 final double gammaAbs = SQRT_TWO_PI / absX *
1197 FastMath.pow(y, absX + 0.5) *
1198 FastMath.exp(-y) * lanczos(absX);
1199 if (x > 0.0) {
1200 ret = gammaAbs;
1201 } else {
1202 /*
1203 * From the reflection formula
1204 * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
1205 * and the recurrence relation
1206 * Gamma(1 - x) = -x * Gamma(-x),
1207 * it is found
1208 * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
1209 */
1210 ret = -FastMath.PI /
1211 (x * FastMath.sin(FastMath.PI * x) * gammaAbs);
1212 }
1213 }
1214 return ret;
1215 }
1216
1217 /**
1218 * Returns the value of Γ(x). Based on the <em>NSWC Library of
1219 * Mathematics Subroutines</em> double precision implementation,
1220 * {@code DGAMMA}.
1221 *
1222 * @param x Argument.
1223 * @param <T> Type of the field elements.
1224 * @return the value of {@code Gamma(x)}.
1225 */
1226 public static <T extends CalculusFieldElement<T>> T gamma(final T x) {
1227 final T one = x.getField().getOne();
1228
1229 if ((x.getReal() == x.rint().getReal()) && (x.getReal() <= 0.0)) {
1230 return one.multiply(Double.NaN);
1231 }
1232
1233 final T ret;
1234 final T absX = x.abs();
1235 if (absX.getReal() <= 20.0) {
1236 if (x.getReal() >= 1.0) {
1237 /*
1238 * From the recurrence relation
1239 * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
1240 * then
1241 * Gamma(t) = 1 / [1 + invGamma1pm1(t - 1)],
1242 * where t = x - n. This means that t must satisfy
1243 * -0.5 <= t - 1 <= 1.5.
1244 */
1245 T prod = one;
1246 T t = x;
1247 while (t.getReal() > 2.5) {
1248 t = t.subtract(1.0);
1249 prod = prod.multiply(t);
1250 }
1251 ret = prod.divide(invGamma1pm1(t.subtract(1.0)).add(1.0));
1252 }
1253 else {
1254 /*
1255 * From the recurrence relation
1256 * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
1257 * then
1258 * Gamma(x + n + 1) = 1 / [1 + invGamma1pm1(x + n)],
1259 * which requires -0.5 <= x + n <= 1.5.
1260 */
1261 T prod = x;
1262 T t = x;
1263 while (t.getReal() < -0.5) {
1264 t = t.add(1.0);
1265 prod = prod.multiply(t);
1266 }
1267 ret = prod.multiply(invGamma1pm1(t).add(1)).reciprocal();
1268 }
1269 }
1270 else {
1271 final T y = absX.add(LANCZOS_G + 0.5);
1272 final T gammaAbs = absX.reciprocal().multiply(SQRT_TWO_PI).multiply(y.pow(absX.add(0.5)))
1273 .multiply(y.negate().exp()).multiply(lanczos(absX));
1274 if (x.getReal() > 0.0) {
1275 ret = gammaAbs;
1276 }
1277 else {
1278 /*
1279 * From the reflection formula
1280 * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
1281 * and the recurrence relation
1282 * Gamma(1 - x) = -x * Gamma(-x),
1283 * it is found
1284 * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
1285 */
1286 ret = x.multiply(x.multiply(FastMath.PI).sin()).multiply(gammaAbs).reciprocal().multiply(-FastMath.PI);
1287 }
1288 }
1289 return ret;
1290 }
1291 }