1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 package org.hipparchus.special;
23
24 import org.hipparchus.CalculusFieldElement;
25 import org.hipparchus.Field;
26 import org.hipparchus.UnitTestUtils;
27 import org.hipparchus.exception.MathIllegalArgumentException;
28 import org.hipparchus.exception.MathIllegalStateException;
29 import org.hipparchus.util.Binary64;
30 import org.hipparchus.util.Binary64Field;
31 import org.hipparchus.util.FastMath;
32 import org.junit.jupiter.api.Test;
33
34 import static org.junit.jupiter.api.Assertions.assertEquals;
35 import static org.junit.jupiter.api.Assertions.assertThrows;
36 import static org.junit.jupiter.api.Assertions.assertTrue;
37
38
39
40 class GammaTest {
41
42 final Field<Binary64> field = Binary64Field.getInstance();
43 final Binary64 zero = field.getZero();
44 final Binary64 one = field.getOne();
45
46 private void testRegularizedGamma(double expected, double a, double x) {
47 double actualP = Gamma.regularizedGammaP(a, x);
48 double actualQ = Gamma.regularizedGammaQ(a, x);
49 UnitTestUtils.customAssertEquals(expected, actualP, 10e-15);
50 UnitTestUtils.customAssertEquals(actualP, 1.0 - actualQ, 10e-15);
51 }
52
53 private <T extends CalculusFieldElement<T>> void testRegularizedGammaField(double expected, T a, T x) {
54 T actualP = Gamma.regularizedGammaP(a, x);
55 T actualQ = Gamma.regularizedGammaQ(a, x);
56 UnitTestUtils.customAssertEquals(expected, actualP.getReal(), 10e-15);
57 UnitTestUtils.customAssertEquals(actualP.getReal(), actualQ.negate().add(1.0).getReal(), 10e-15);
58 }
59
60 private void testRegularizedGamma(double expected, double a, double x, double epsilon, int maxIterations) {
61 double actualP = Gamma.regularizedGammaP(a, x, epsilon, maxIterations);
62 double actualQ = Gamma.regularizedGammaQ(a, x, epsilon, maxIterations);
63 UnitTestUtils.customAssertEquals(expected, actualP, epsilon);
64 UnitTestUtils.customAssertEquals(actualP, 1.0 - actualQ, epsilon);
65 }
66
67 private <T extends CalculusFieldElement<T>> void testRegularizedGammaField(double expected, T a, T x, double epsilon, int maxIterations) {
68 T actualP = Gamma.regularizedGammaP(a, x, epsilon, maxIterations);
69 T actualQ = Gamma.regularizedGammaQ(a, x, epsilon, maxIterations);
70 UnitTestUtils.customAssertEquals(expected, actualP.getReal(), epsilon);
71 UnitTestUtils.customAssertEquals(actualP.getReal(), actualQ.negate().add(1.0).getReal(), epsilon);
72 }
73
74 private void testLogGamma(double expected, double x) {
75 double actual = Gamma.logGamma(x);
76 UnitTestUtils.customAssertEquals(expected, actual, 10e-15);
77 }
78
79 private <T extends CalculusFieldElement<T>> void testLogGammaField(T expected, T x) {
80 T actual = Gamma.logGamma(x);
81 UnitTestUtils.customAssertEquals(expected.getReal(), actual.getReal(), 10e-15);
82 }
83
84 @Test
85 void testRegularizedGammaNanPositive() {
86 testRegularizedGamma(Double.NaN, Double.NaN, 1.0);
87 }
88
89 @Test
90 void testRegularizedGammaNanPositiveField() {
91 testRegularizedGammaField(Double.NaN, new Binary64(Double.NaN), one);
92 }
93
94 @Test
95 void testRegularizedGammaPositiveNan() {
96 testRegularizedGamma(Double.NaN, 1.0, Double.NaN);
97 }
98
99 @Test
100 void testRegularizedGammaPositiveNanField() {
101 testRegularizedGammaField(Double.NaN, one, new Binary64(Double.NaN));
102 }
103
104 @Test
105 void testRegularizedGammaNegativePositive() {
106 testRegularizedGamma(Double.NaN, -1.5, 1.0);
107 }
108
109 @Test
110 void testRegularizedGammaNegativePositiveField() {
111 testRegularizedGammaField(Double.NaN, new Binary64(-1.5), one);
112 }
113
114 @Test
115 void testRegularizedGammaPositiveNegative() {
116 testRegularizedGamma(Double.NaN, 1.0, -1.0);
117 }
118
119 @Test
120 void testRegularizedGammaPositiveNegativeField() {
121 testRegularizedGammaField(Double.NaN, one, one.negate());
122 }
123
124 @Test
125 void testRegularizedGammaZeroPositive() {
126 testRegularizedGamma(Double.NaN, 0.0, 1.0);
127 }
128
129 @Test
130 void testRegularizedGammaZeroPositiveField() {
131 testRegularizedGammaField(Double.NaN, zero, one);
132 }
133
134 @Test
135 void testRegularizedGammaPositiveZero() {
136 testRegularizedGamma(0.0, 1.0, 0.0);
137 }
138
139 @Test
140 void testRegularizedGammaPositiveZeroField() {
141 testRegularizedGammaField(0.0, one, zero);
142 }
143
144 @Test
145 void testRegularizedGammaPositivePositive() {
146 testRegularizedGamma(0.632120558828558, 1.0, 1.0);
147 }
148
149 @Test
150 void testRegularizedGammaPositivePositiveField() {
151 testRegularizedGammaField(0.632120558828558, one, one);
152 }
153
154 @Test
155 void testRegularizedGammaPMaxNumberOfIterationsExceeded(){
156 assertThrows(MathIllegalStateException.class, () -> {
157 testRegularizedGamma(0.632120558828558, 1.0, 1.0, 1e-15, 1);
158 });
159 }
160
161 @Test
162 void testRegularizedGammaPMaxNumberOfIterationsExceededField(){
163 assertThrows(MathIllegalStateException.class, () -> {
164 testRegularizedGammaField(0.632120558828558, one, one, 1e-15, 1);
165 });
166 }
167
168 @Test
169 void testLogGammaNan() {
170 testLogGamma(Double.NaN, Double.NaN);
171 }
172
173 @Test
174 void testLogGammaNanField() {
175 testLogGammaField(new Binary64(Double.NaN), new Binary64(Double.NaN));
176 }
177
178 @Test
179 void testLogGammaNegative() {
180 testLogGamma(Double.NaN, -1.0);
181 }
182
183 @Test
184 void testLogGammaNegativeField() {
185 testLogGammaField(new Binary64(Double.NaN), one.negate());
186 }
187
188 @Test
189 void testLogGammaZero() {
190 testLogGamma(Double.NaN, 0.0);
191 }
192
193 @Test
194 void testLogGammaZeroField() {
195 testLogGammaField(new Binary64(Double.NaN), zero);
196 }
197
198 @Test
199 void testLogGammaPositive() {
200 testLogGamma(0.6931471805599457, 3.0);
201 }
202
203 @Test
204 void testLogGammaPositiveField() {
205 testLogGammaField(new Binary64(0.6931471805599457), new Binary64(3.0));
206 }
207
208 @Test
209 void testDigammaLargeArgs() {
210 double eps = 1e-8;
211 assertEquals(4.6001618527380874002, Gamma.digamma(100), eps);
212 assertEquals(3.9019896734278921970, Gamma.digamma(50), eps);
213 assertEquals(2.9705239922421490509, Gamma.digamma(20.), eps);
214 assertEquals(2.9958363947076465821, Gamma.digamma(20.5), eps);
215 assertEquals(2.2622143570941481605, Gamma.digamma(10.1), eps);
216 assertEquals(2.1168588189004379233, Gamma.digamma(8.8), eps);
217 assertEquals(1.8727843350984671394, Gamma.digamma(7), eps);
218 assertEquals(0.42278433509846713939, Gamma.digamma(2), eps);
219 assertEquals(-100.56088545786867450, Gamma.digamma(0.01), eps);
220 assertEquals(-4.0390398965921882955, Gamma.digamma(-0.8), eps);
221 assertEquals(4.2003210041401844726, Gamma.digamma(-6.3), eps);
222 }
223
224 @Test
225 void testDigammaLargeArgsField() {
226 double eps = 1e-8;
227 assertEquals(4.6001618527380874002, Gamma.digamma(new Binary64(100)).getReal(), eps);
228 assertEquals(3.9019896734278921970, Gamma.digamma(new Binary64(50)).getReal(), eps);
229 assertEquals(2.9705239922421490509, Gamma.digamma(new Binary64(20)).getReal(), eps);
230 assertEquals(2.9958363947076465821, Gamma.digamma(new Binary64(20.5)).getReal(), eps);
231 assertEquals(2.2622143570941481605, Gamma.digamma(new Binary64(10.1)).getReal(), eps);
232 assertEquals(2.1168588189004379233, Gamma.digamma(new Binary64(8.8)).getReal(), eps);
233 assertEquals(1.8727843350984671394, Gamma.digamma(new Binary64(7)).getReal(), eps);
234 assertEquals(0.42278433509846713939, Gamma.digamma(new Binary64(2)).getReal(), eps);
235 assertEquals(-100.56088545786867450, Gamma.digamma(new Binary64(0.01)).getReal(), eps);
236 assertEquals(-4.0390398965921882955, Gamma.digamma(new Binary64(-0.8)).getReal(), eps);
237 assertEquals(4.2003210041401844726, Gamma.digamma(new Binary64(-6.3)).getReal(), eps);
238 }
239
240 @Test
241 void testDigammaSmallArgs() {
242
243
244 double[] expected = {-10.423754940411076795, -100.56088545786867450, -1000.5755719318103005,
245 -10000.577051183514335, -100000.57719921568107, -1.0000005772140199687e6, -1.0000000577215500408e7,
246 -1.0000000057721564845e8, -1.0000000005772156633e9, -1.0000000000577215665e10, -1.0000000000057721566e11,
247 -1.0000000000005772157e12, -1.0000000000000577216e13, -1.0000000000000057722e14, -1.0000000000000005772e15, -1e+16,
248 -1e+17, -1e+18, -1e+19, -1e+20, -1e+21, -1e+22, -1e+23, -1e+24, -1e+25, -1e+26,
249 -1e+27, -1e+28, -1e+29, -1e+30};
250 for (double n = 1; n < 30; n++) {
251 checkRelativeError(String.format("Test %.0f: ", n), expected[(int) (n - 1)], Gamma.digamma(FastMath.pow(10.0, -n)), 1e-8);
252 }
253 }
254
255 @Test
256 void testDigammaSmallArgsField() {
257
258
259 double[] expected = {-10.423754940411076795, -100.56088545786867450, -1000.5755719318103005,
260 -10000.577051183514335, -100000.57719921568107, -1.0000005772140199687e6, -1.0000000577215500408e7,
261 -1.0000000057721564845e8, -1.0000000005772156633e9, -1.0000000000577215665e10, -1.0000000000057721566e11,
262 -1.0000000000005772157e12, -1.0000000000000577216e13, -1.0000000000000057722e14, -1.0000000000000005772e15, -1e+16,
263 -1e+17, -1e+18, -1e+19, -1e+20, -1e+21, -1e+22, -1e+23, -1e+24, -1e+25, -1e+26,
264 -1e+27, -1e+28, -1e+29, -1e+30};
265 for (double n = 1; n < 30; n++) {
266 checkRelativeError(String.format("Test %.0f: ", n), expected[(int) (n - 1)], Gamma.digamma(new Binary64(FastMath.pow(10.0, -n))).getReal(), 1e-8);
267 }
268 }
269
270 @Test
271 void testDigammaNonRealArgs() {
272 assertTrue(Double.isNaN(Gamma.digamma(Double.NaN)));
273 assertTrue(Double.isInfinite(Gamma.digamma(Double.POSITIVE_INFINITY)));
274 assertTrue(Double.isInfinite(Gamma.digamma(Double.NEGATIVE_INFINITY)));
275 }
276
277 @Test
278 void testDigammaNonRealArgsField() {
279 assertTrue(Gamma.digamma(new Binary64(Double.NaN)).isNaN());
280 assertTrue(Gamma.digamma(new Binary64(Double.POSITIVE_INFINITY)).isInfinite());
281 assertTrue(Gamma.digamma(new Binary64(Double.NEGATIVE_INFINITY)).isInfinite());
282 }
283
284 @Test
285 void testTrigamma() {
286 double eps = 1e-13;
287
288
289
290 double[] data = {
291 1e-6, 1.0000000000016449316627376e12,
292 1e-4, 1.0000000164469368793e8,
293 1e-3, 1.0000016425331958690e6,
294 1e-2, 10001.621213528313220,
295 1e-1, 101.43329915079275882,
296 1, 1.6449340668482264365,
297 2, 0.64493406684822643647,
298 3, 0.39493406684822643647,
299 4, 0.28382295573711532536,
300 5, 0.22132295573711532536,
301 10, 0.10516633568168574612,
302 20, 0.051270822935203119832,
303 50, 0.020201333226697125806,
304 100, 0.010050166663333571395
305 };
306 for (int i = data.length - 2; i >= 0; i -= 2) {
307 assertEquals(data[i + 1], Gamma.trigamma(data[i]), eps, String.format("trigamma %.0f", data[i]));
308 }
309 }
310
311 @Test
312 void testTrigammaField() {
313 double eps = 1e-13;
314
315
316
317 Binary64[] data = {
318 new Binary64(1e-6), new Binary64(1.00000000000164493166e12) ,
319 new Binary64(1e-4), new Binary64(1.0000000164469368793e8) ,
320 new Binary64(1e-3), new Binary64(1.0000016425331958690e6),
321 new Binary64(1e-2), new Binary64(10001.621213528313220),
322 new Binary64(1e-1), new Binary64(101.43329915079275882),
323 new Binary64(1), new Binary64(1.6449340668482264365),
324 new Binary64(2), new Binary64(0.64493406684822643647),
325 new Binary64(3), new Binary64(0.39493406684822643647),
326 new Binary64(4), new Binary64(0.28382295573711532536),
327 new Binary64(5), new Binary64(0.22132295573711532536),
328 new Binary64(10), new Binary64(0.10516633568168574612),
329 new Binary64(20), new Binary64(0.051270822935203119832),
330 new Binary64(50), new Binary64(0.020201333226697125806),
331 new Binary64(100), new Binary64(0.010050166663333571395)
332 };
333 for (int i = data.length - 2; i >= 0; i -= 2) {
334 assertEquals(data[i + 1].getReal(),
335 Gamma.trigamma(data[i]).getReal(), eps, String.format("trigamma %.6f", data[i].getReal()));
336 }
337 }
338
339 @Test
340 void testTrigammaSmallArg() {
341 double eps = 2;
342 assertEquals(1e16,
343 Gamma.trigamma(1e-8), eps, String.format("trigamma %.8f", 1e-8));
344 }
345
346 @Test
347 void testTrigammaSmallArgField() {
348 double eps = 2;
349 assertEquals(1e16,
350 Gamma.trigamma(new Binary64(1e-8)).getReal(), eps, String.format("trigamma %.8f", 1e-8));
351 }
352
353 @Test
354 void testTrigammaNonRealArgs() {
355 assertTrue(Double.isNaN(Gamma.trigamma(Double.NaN)));
356 assertTrue(Double.isInfinite(Gamma.trigamma(Double.POSITIVE_INFINITY)));
357 assertTrue(Double.isInfinite(Gamma.trigamma(Double.NEGATIVE_INFINITY)));
358 }
359
360 @Test
361 void testTrigammaNonRealArgsField() {
362 assertTrue(Gamma.trigamma(new Binary64(Double.NaN)).isNaN());
363 assertTrue(Gamma.trigamma(new Binary64(Double.POSITIVE_INFINITY)).isInfinite());
364 assertTrue(Gamma.trigamma(new Binary64(Double.NEGATIVE_INFINITY)).isInfinite());
365 }
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384 private static final double[][] LOG_GAMMA_REF = {
385 { 0.125 , 2.019418357553796 },
386 { 0.25 , 1.288022524698077 },
387 { 0.375 , .8630739822706475 },
388 { 0.5 , .5723649429247001 },
389 { 0.625 , .3608294954889402 },
390 { 0.75 , .2032809514312954 },
391 { 0.875 , .08585870722533433 },
392 { 0.890625 , .07353860936979656 },
393 { 0.90625 , .06169536624059108 },
394 { 0.921875 , .05031670080005688 },
395 { 0.9375 , 0.0393909017345823 },
396 { 0.953125 , .02890678734595923 },
397 { 0.96875 , .01885367233441289 },
398 { 0.984375 , .009221337197578781 },
399 { 1.0 , 0.0 },
400 { 1.015625 , - 0.00881970970573307 },
401 { 1.03125 , - .01724677500176807 },
402 { 1.046875 , - .02528981394675729 },
403 { 1.0625 , - .03295710029357782 },
404 { 1.078125 , - .04025658272400143 },
405 { 1.09375 , - .04719590272716985 },
406 { 1.109375 , - .05378241123619192 },
407 { 1.125 , - .06002318412603958 },
408 { 1.25 , - .09827183642181316 },
409 { 1.375 , - .1177552707410788 },
410 { 1.5 , - .1207822376352452 },
411 { 1.625 , - .1091741337567954 },
412 { 1.75 , - .08440112102048555 },
413 { 1.875 , - 0.0476726853991883 },
414 { 1.890625 , - .04229320615532515 },
415 { 1.90625 , - .03674470657266143 },
416 { 1.921875 , - .03102893865389552 },
417 { 1.9375 , - .02514761940298887 },
418 { 1.953125 , - .01910243184040138 },
419 { 1.96875 , - .01289502598016741 },
420 { 1.984375 , - .006527019770560387 },
421 { 2.0 , 0.0 },
422 { 2.015625 , .006684476830232185 },
423 { 2.03125 , .01352488366498562 },
424 { 2.046875 , .02051972208453692 },
425 { 2.0625 , .02766752152285702 },
426 { 2.078125 , 0.0349668385135861 },
427 { 2.09375 , .04241625596251728 },
428 { 2.109375 , .05001438244545164 },
429 { 2.125 , .05775985153034387 },
430 { 2.25 , .1248717148923966 },
431 { 2.375 , .2006984603774558 },
432 { 2.5 , .2846828704729192 },
433 { 2.625 , .3763336820249054 },
434 { 2.75 , .4752146669149371 },
435 { 2.875 , .5809359740231859 },
436 { 2.890625 , .5946142560817441 },
437 { 2.90625 , .6083932548009232 },
438 { 2.921875 , .6222723333588501 },
439 { 2.9375 , .6362508628423761 },
440 { 2.953125 , .6503282221022278 },
441 { 2.96875 , .6645037976116387 },
442 { 2.984375 , 0.678776983328359 },
443 { 3.0 , .6931471805599453 },
444 { 3.015625 , .7076137978322324 },
445 { 3.03125 , .7221762507608962 },
446 { 3.046875 , .7368339619260166 },
447 { 3.0625 , 0.751586360749556 },
448 { 3.078125 , .7664328833756681 },
449 { 3.09375 , .7813729725537568 },
450 { 3.109375 , .7964060775242092 },
451 { 3.125 , 0.811531653906724 },
452 { 3.25 , .9358019311087253 },
453 { 3.375 , 1.06569589786406 },
454 { 3.5 , 1.200973602347074 },
455 { 3.625 , 1.341414578068493 },
456 { 3.75 , 1.486815578593417 },
457 { 3.875 , 1.6369886482725 },
458 { 4.0 , 1.791759469228055 },
459 { 4.125 , 1.950965937095089 },
460 { 4.25 , 2.114456927450371 },
461 { 4.375 , 2.282091222188554 },
462 { 4.5 , 2.453736570842442 },
463 { 4.625 , 2.62926886637513 },
464 { 4.75 , 2.808571418575736 },
465 { 4.875 , 2.99153431107781 },
466 { 5.0 , 3.178053830347946 },
467 { 5.125 , 3.368031956881733 },
468 { 5.25 , 3.561375910386697 },
469 { 5.375 , 3.757997741998131 },
470 { 5.5 , 3.957813967618717 },
471 { 5.625 , 4.160745237339519 },
472 { 5.75 , 4.366716036622286 },
473 { 5.875 , 4.57565441552762 },
474 { 6.0 , 4.787491742782046 },
475 { 6.125 , 5.002162481906205 },
476 { 6.25 , 5.219603986990229 },
477 { 6.375 , 5.439756316011858 },
478 { 6.5 , 5.662562059857142 },
479 { 6.625 , 5.887966185430003 },
480 { 6.75 , 6.115915891431546 },
481 { 6.875 , 6.346360475557843 },
482 { 7.0 , 6.579251212010101 },
483 { 7.125 , 6.814541238336996 },
484 { 7.25 , 7.05218545073854 },
485 { 7.375 , 7.292140407056348 },
486 { 7.5 , 7.534364236758733 },
487 { 7.625 , 7.778816557302289 },
488 { 7.75 , 8.025458396315983 },
489 { 7.875 , 8.274252119110479 },
490 { 8.0 , 8.525161361065415 },
491 { 8.125 , 8.77815096449171 },
492 { 8.25 , 9.033186919605123 },
493 { 8.375 , 9.290236309282232 },
494 { 8.5 , 9.549267257300997 },
495 { 8.625 , 9.810248879795765 },
496 { 8.75 , 10.07315123968124 },
497 { 8.875 , 10.33794530382217 },
498 { 9.0 , 10.60460290274525 },
499 { 9.125 , 10.87309669270751 },
500 { 9.25 , 11.14340011995171 },
501 { 9.375 , 11.41548738699336 },
502 { 9.5 , 11.68933342079727 },
503 { 9.625 , 11.96491384271319 },
504 { 9.75 , 12.24220494005076 },
505 { 9.875 , 12.52118363918365 },
506 { 10.0 , 12.80182748008147 },
507 { 0.8 , .1520596783998376 },
508 { 100.0 , 359.1342053695754 },
509 { 1000.0 , 5905.220423209181 },
510 { 10000.0 , 82099.71749644238 },
511 { 100000.0 , 1051287.708973657 },
512 { 1000000.0 , 1.2815504569147612e+7 },
513 { 10000000.0 , 1.511809493694739e+8 },
514 { 1.e+8 , 1.7420680661038346e+9 },
515 { 1.e+9 , 1.972326582750371e+10 },
516 { 1.e+10 , 2.202585092888106e+11 },
517 };
518
519 @Test
520 void testLogGamma() {
521 final int ulps = 3;
522 for (int i = 0; i < LOG_GAMMA_REF.length; i++) {
523 final double[] data = LOG_GAMMA_REF[i];
524 final double x = data[0];
525 final double expected = data[1];
526 final double actual = Gamma.logGamma(x);
527 final double tol;
528 if (expected == 0.0) {
529 tol = 1E-15;
530 } else {
531 tol = ulps * FastMath.ulp(expected);
532 }
533 assertEquals(expected, actual, tol, Double.toString(x));
534 }
535 }
536
537 @Test
538 void testLogGammaField() {
539 final int ulps = 3;
540 for (int i = 0; i < LOG_GAMMA_REF.length; i++) {
541 final double[] data = LOG_GAMMA_REF[i];
542 final Binary64 x = new Binary64(data[0]);
543 final double expected = data[1];
544 final double actual = Gamma.logGamma(x).getReal();
545 final double tol;
546 if (expected == 0.0) {
547 tol = 1E-15;
548 } else {
549 tol = ulps * FastMath.ulp(expected);
550 }
551 assertEquals(expected, actual, tol, Double.toString(x.getReal()));
552 }
553 }
554
555 @Test
556 void testLogGammaPrecondition1() {
557 assertTrue(Double.isNaN(Gamma.logGamma(0.0)));
558 }
559
560 @Test
561 void testLogGammaPrecondition1Field() {
562 assertTrue(Gamma.logGamma(new Binary64(0.0)).isNaN());
563 }
564
565 @Test
566 void testLogGammaPrecondition2() {
567 assertTrue(Double.isNaN(Gamma.logGamma(-1.0)));
568 }
569
570 @Test
571 void testLogGammaPrecondition2Field() {
572 assertTrue(Gamma.logGamma(new Binary64(-1.0)).isNaN());
573 }
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596 private static final double[][] INV_GAMMA1P_M1_REF = {
597 { -0.5 , -.4358104164522437 },
598 { -0.375 , -.3029021533379859 },
599 { -0.25 , -0.183951060901737 },
600 { -0.125 , -.08227611018520711 },
601 { 0.0 , 0.0 },
602 { 0.125 , .06186116458306091 },
603 { 0.25 , .1032626513208373 },
604 { 0.375 , .1249687649039041 },
605 { 0.5 , .1283791670955126 },
606 { 0.625 , .1153565546592225 },
607 { 0.75 , 0.0880652521310173 },
608 { 0.875 , .04882730264547758 },
609 { 1.0 , 0.0 },
610 { 1.125 , -.05612340925950141 },
611 { 1.25 , -.1173898789433302 },
612 { 1.375 , -.1818408982517061 },
613 { 1.5 , -0.247747221936325 },
614 };
615
616 @Test
617 void testInvGamma1pm1() {
618
619 final int ulps = 3;
620 for (int i = 0; i < INV_GAMMA1P_M1_REF.length; i++) {
621 final double[] ref = INV_GAMMA1P_M1_REF[i];
622 final double x = ref[0];
623 final double expected = ref[1];
624 final double actual = Gamma.invGamma1pm1(x);
625 final double tol = ulps * FastMath.ulp(expected);
626 assertEquals(expected, actual, tol, Double.toString(x));
627 }
628 }
629
630 @Test
631 void testInvGamma1pm1Field() {
632
633 final int ulps = 3;
634 for (int i = 0; i < INV_GAMMA1P_M1_REF.length; i++) {
635 final double[] ref = INV_GAMMA1P_M1_REF[i];
636 final Binary64 x = new Binary64(ref[0]);
637 final double expected = ref[1];
638 final double actual = Gamma.invGamma1pm1(x).getReal();
639 final double tol = ulps * FastMath.ulp(expected);
640 assertEquals(expected, actual, tol, Double.toString(x.getReal()));
641 }
642 }
643
644 @Test
645 void testInvGamma1pm1Precondition1() {
646 assertThrows(MathIllegalArgumentException.class, () -> {
647
648 Gamma.invGamma1pm1(-0.51);
649 });
650 }
651
652 @Test
653 void testInvGamma1pm1Precondition2() {
654 assertThrows(MathIllegalArgumentException.class, () -> {
655
656 Gamma.invGamma1pm1(1.51);
657 });
658 }
659
660 @Test
661 void testInvGamma1pm1Precondition1Field() {
662 assertThrows(MathIllegalArgumentException.class, () -> {
663
664 Gamma.invGamma1pm1(new Binary64(-0.51));
665 });
666 }
667
668 @Test
669 void testInvGamma1pm1Precondition2Field() {
670 assertThrows(MathIllegalArgumentException.class, () -> {
671
672 Gamma.invGamma1pm1(new Binary64(1.51));
673 });
674 }
675
676 private static final double[][] LOG_GAMMA1P_REF = {
677 { - 0.5 , .5723649429247001 },
678 { - 0.375 , .3608294954889402 },
679 { - 0.25 , .2032809514312954 },
680 { - 0.125 , .08585870722533433 },
681 { 0.0 , 0.0 },
682 { 0.125 , - .06002318412603958 },
683 { 0.25 , - .09827183642181316 },
684 { 0.375 , - .1177552707410788 },
685 { 0.5 , - .1207822376352452 },
686 { 0.625 , - .1091741337567954 },
687 { 0.75 , - .08440112102048555 },
688 { 0.875 , - 0.0476726853991883 },
689 { 1.0 , 0.0 },
690 { 1.125 , .05775985153034387 },
691 { 1.25 , .1248717148923966 },
692 { 1.375 , .2006984603774558 },
693 { 1.5 , .2846828704729192 },
694 };
695
696 @Test
697 void testLogGamma1p() {
698
699 final int ulps = 3;
700 for (int i = 0; i < LOG_GAMMA1P_REF.length; i++) {
701 final double[] ref = LOG_GAMMA1P_REF[i];
702 final double x = ref[0];
703 final double expected = ref[1];
704 final double actual = Gamma.logGamma1p(x);
705 final double tol = ulps * FastMath.ulp(expected);
706 assertEquals(expected, actual, tol, Double.toString(x));
707 }
708 }
709
710 @Test
711 void testLogGamma1pField() {
712
713 final int ulps = 3;
714 for (int i = 0; i < LOG_GAMMA1P_REF.length; i++) {
715 final double[] ref = LOG_GAMMA1P_REF[i];
716 final Binary64 x = new Binary64(ref[0]);
717 final double expected = ref[1];
718 final double actual = Gamma.logGamma1p(x).getReal();
719 final double tol = ulps * FastMath.ulp(expected);
720 assertEquals(expected, actual, tol, Double.toString(x.getReal()));
721 }
722 }
723
724 @Test
725 void testLogGamma1pPrecondition1() {
726 assertThrows(MathIllegalArgumentException.class, () -> {
727
728 Gamma.logGamma1p(-0.51);
729 });
730 }
731
732 @Test
733 void testLogGamma1pPrecondition1Field() {
734 assertThrows(MathIllegalArgumentException.class, () -> {
735
736 Gamma.logGamma1p(new Binary64(-0.51));
737 });
738 }
739
740 @Test
741 void testLogGamma1pPrecondition2() {
742 assertThrows(MathIllegalArgumentException.class, () -> {
743
744 Gamma.logGamma1p(1.51);
745 });
746 }
747
748 @Test
749 void testLogGamma1pPrecondition2Field() {
750 assertThrows(MathIllegalArgumentException.class, () -> {
751
752 Gamma.logGamma1p(new Binary64(1.51));
753 });
754 }
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776 private static final double[][] GAMMA_REF = {
777 { - 19.875 , 4.920331854832504e-18 },
778 { - 19.75 , 3.879938752480031e-18 },
779 { - 19.625 , 4.323498423815027e-18 },
780 { - 19.5 , 5.811045977502237e-18 },
781 { - 19.375 , 9.14330910942125e-18 },
782 { - 19.25 , 1.735229114436739e-17 },
783 { - 19.125 , 4.653521565668223e-17 },
784 { - 18.875 , - 9.779159561479603e-17 },
785 { - 18.75 , - 7.662879036148061e-17 },
786 { - 18.625 , - 8.484865656736991e-17 },
787 { - 18.5 , - 1.133153965612936e-16 },
788 { - 18.375 , - 1.771516139950367e-16 },
789 { - 18.25 , - 3.340316045290721e-16 },
790 { - 18.125 , - 8.899859994340475e-16 },
791 { - 17.875 , 1.845816367229275e-15 },
792 { - 17.75 , 1.436789819277761e-15 },
793 { - 17.625 , 1.580306228567265e-15 },
794 { - 17.5 , 2.096334836383932e-15 },
795 { - 17.375 , 3.255160907158799e-15 },
796 { - 17.25 , 6.096076782655566e-15 },
797 { - 17.125 , 1.613099623974211e-14 },
798 { - 16.875 , - 3.29939675642233e-14 },
799 { - 16.75 , - 2.550301929218027e-14 },
800 { - 16.625 , - 2.785289727849803e-14 },
801 { - 16.5 , - 3.66858596367188e-14 },
802 { - 16.375 , - 5.655842076188414e-14 },
803 { - 16.25 , - 1.051573245008085e-13 },
804 { - 16.125 , - 2.762433106055837e-13 },
805 { - 15.875 , 5.567732026462681e-13 },
806 { - 15.75 , 4.271755731440195e-13 },
807 { - 15.625 , 4.630544172550298e-13 },
808 { - 15.5 , 6.053166840058604e-13 },
809 { - 15.375 , 9.261441399758529e-13 },
810 { - 15.25 , 1.708806523138138e-12 },
811 { - 15.125 , 4.454423383515037e-12 },
812 { - 14.875 , - 8.838774592009505e-12 },
813 { - 14.75 , - 6.728015277018307e-12 },
814 { - 14.625 , - 7.235225269609841e-12 },
815 { - 14.5 , - 9.382408602090835e-12 },
816 { - 14.375 , - 1.423946615212874e-11 },
817 { - 14.25 , - 2.605929947785661e-11 },
818 { - 14.125 , - 6.737315367566492e-11 },
819 { - 13.875 , 1.314767720561414e-10 },
820 { - 13.75 , 9.923822533602004e-11 },
821 { - 13.625 , 1.058151695680439e-10 },
822 { - 13.5 , 1.360449247303171e-10 },
823 { - 13.375 , 2.046923259368506e-10 },
824 { - 13.25 , 3.713450175594567e-10 },
825 { - 13.125 , 9.516457956687671e-10 },
826 { - 12.875 , - 1.8242402122789617e-9 },
827 { - 12.75 , - 1.3645255983702756e-9 },
828 { - 12.625 , - 1.4417316853645984e-9 },
829 { - 12.5 , - 1.836606483859281e-9 },
830 { - 12.375 , - 2.7377598594053765e-9 },
831 { - 12.25 , - 4.9203214826628017e-9 },
832 { - 12.125 , - 1.2490351068152569e-8 },
833 { - 11.875 , 2.3487092733091633e-8 },
834 { - 11.75 , 1.7397701379221012e-8 },
835 { - 11.625 , 1.8201862527728055e-8 },
836 { - 11.5 , 2.295758104824101e-8 },
837 { - 11.375 , 3.3879778260141535e-8 },
838 { - 11.25 , 6.027393816261931e-8 },
839 { - 11.125 , 1.5144550670134987e-7 },
840 { - 10.875 , - 2.7890922620546316e-7 },
841 { - 10.75 , - 2.044229912058469e-7 },
842 { - 10.625 , - 2.1159665188483867e-7 },
843 { - 10.5 , - 2.640121820547716e-7 },
844 { - 10.375 , - 3.8538247770911e-7 },
845 { - 10.25 , - 6.780818043294673e-7 },
846 { - 10.125 , - 1.6848312620525174e-6 },
847 { - 9.875 , 3.0331378349844124e-6 },
848 { - 9.75 , 2.1975471554628537e-6 },
849 { - 9.625 , 2.2482144262764103e-6 },
850 { - 9.5 , 2.772127911575102e-6 },
851 { - 9.375 , 3.998343206232017e-6 },
852 { - 9.25 , 6.95033849437704e-6 },
853 { - 9.125 , 1.7058916528281737e-5 },
854 { - 8.875 , - 2.9952236120471065e-5 },
855 { - 8.75 , - 2.1426084765762826e-5 },
856 { - 8.625 , - 2.163906385291045e-5 },
857 { - 8.5 , - 2.633521515996347e-5 },
858 { - 8.375 , - 3.748446755842515e-5 },
859 { - 8.25 , - 6.429063107298763e-5 },
860 { - 8.125 , - 1.5566261332057085e-4 },
861 { - 7.875 , 2.658260955691807e-4 },
862 { - 7.75 , 1.874782417004247e-4 },
863 { - 7.625 , 1.8663692573135265e-4 },
864 { - 7.5 , 2.238493288596895e-4 },
865 { - 7.375 , 3.1393241580181064e-4 },
866 { - 7.25 , 5.303977063521479e-4 },
867 { - 7.125 , .001264758733229638 },
868 { - 6.875 , - .002093380502607298 },
869 { - 6.75 , - .001452956373178292 },
870 { - 6.625 , - .001423106558701564 },
871 { - 6.5 , - .001678869966447671 },
872 { - 6.375 , - .002315251566538353 },
873 { - 6.25 , - .003845383371053072 },
874 { - 6.125 , - .009011405974261174 },
875 { - 5.875 , .01439199095542518 },
876 { - 5.75 , .009807455518953468 },
877 { - 5.625 , .009428080951397862 },
878 { - 5.5 , .01091265478190986 },
879 { - 5.375 , 0.014759728736682 },
880 { - 5.25 , 0.0240336460690817 },
881 { - 5.125 , .05519486159234969 },
882 { - 4.875 , - 0.0845529468631229 },
883 { - 4.75 , - .05639286923398244 },
884 { - 4.625 , - .05303295535161297 },
885 { - 4.5 , - .06001960130050425 },
886 { - 4.375 , - .07933354195966577 },
887 { - 4.25 , - .1261766418626789 },
888 { - 4.125 , - .2828736656607921 },
889 { - 3.875 , .4121956159577241 },
890 { - 3.75 , .2678661288614166 },
891 { - 3.625 , 0.24527741850121 },
892 { - 3.5 , .2700882058522691 },
893 { - 3.375 , .3470842460735378 },
894 { - 3.25 , .5362507279163854 },
895 { - 3.125 , 1.166853870850768 },
896 { - 2.875 , - 1.597258011836181 },
897 { - 2.75 , - 1.004497983230312 },
898 { - 2.625 , - .8891306420668862 },
899 { - 2.5 , - .9453087204829419 },
900 { - 2.375 , - 1.17140933049819 },
901 { - 2.25 , - 1.742814865728253 },
902 { - 2.125 , - 3.646418346408649 },
903 { - 1.875 , 4.59211678402902 },
904 { - 1.75 , 2.762369453883359 },
905 { - 1.625 , 2.333967935425576 },
906 { - 1.5 , 2.363271801207355 },
907 { - 1.375 , 2.782097159933201 },
908 { - 1.25 , 3.921333447888569 },
909 { - 1.125 , 7.748638986118379 },
910 { - 0.875 , - 8.610218970054413 },
911 { - 0.75 , - 4.834146544295877 },
912 { - 0.625 , - 3.792697895066561 },
913 { - 0.5 , - 3.544907701811032 },
914 { - 0.375 , - 3.825383594908152 },
915 { - 0.25 , - 4.901666809860711 },
916 { - 0.125 , - 8.717218859383175 },
917 { 0.125 , 7.533941598797612 },
918 { 0.25 , 3.625609908221908 },
919 { 0.375 , 2.370436184416601 },
920 { 0.5 , 1.772453850905516 },
921 { 0.625 , 1.434518848090557 },
922 { 0.75 , 1.225416702465178 },
923 { 0.875 , 1.089652357422897 },
924 { 1.0 , 1.0 },
925 { 1.125 , .9417426998497015 },
926 { 1.25 , 0.906402477055477 },
927 { 1.375 , .8889135691562253 },
928 { 1.5 , 0.886226925452758 },
929 { 1.625 , 0.896574280056598 },
930 { 1.75 , .9190625268488832 },
931 { 1.875 , .9534458127450348 },
932 { 2.0 , 1.0 },
933 { 2.125 , 1.059460537330914 },
934 { 2.25 , 1.133003096319346 },
935 { 2.375 , 1.22225615758981 },
936 { 2.5 , 1.329340388179137 },
937 { 2.625 , 1.456933205091972 },
938 { 2.75 , 1.608359421985546 },
939 { 2.875 , 1.78771089889694 },
940 { 3.0 , 2.0 },
941 { 3.125 , 2.251353641828193 },
942 { 3.25 , 2.549256966718529 },
943 { 3.375 , 2.902858374275799 },
944 { 3.5 , 3.323350970447843 },
945 { 3.625 , 3.824449663366426 },
946 { 3.75 , 4.422988410460251 },
947 { 3.875 , 5.139668834328703 },
948 { 4.0 , 6.0 },
949 { 4.125 , 7.035480130713102 },
950 { 4.25 , 8.28508514183522 },
951 { 4.375 , 9.797147013180819 },
952 { 4.5 , 11.63172839656745 },
953 { 4.625 , 13.86363002970329 },
954 { 4.75 , 16.58620653922594 },
955 { 4.875 , 19.91621673302373 },
956 { 5.0 , 24.0 },
957 { 5.125 , 29.02135553919155 },
958 { 5.25 , 35.21161185279968 },
959 { 5.375 , 42.86251818266609 },
960 { 5.5 , 52.34277778455352 },
961 { 5.625 , 64.11928888737773 },
962 { 5.75 , 78.78448106132322 },
963 { 5.875 , 97.09155657349066 },
964 { 6.0 , 120.0 },
965 { 6.125 , 148.7344471383567 },
966 { 6.25 , 184.8609622271983 },
967 { 6.375 , 230.3860352318302 },
968 { 6.5 , 287.8852778150443 },
969 { 6.625 , 360.6709999914997 },
970 { 6.75 , 453.0107661026085 },
971 { 6.875 , 570.4128948692577 },
972 { 7.0 , 720.0 },
973 { 7.125 , 910.9984887224346 },
974 { 7.25 , 1155.38101391999 },
975 { 7.375 , 1468.710974602918 },
976 { 7.5 , 1871.254305797788 },
977 { 7.625 , 2389.445374943686 },
978 { 7.75 , 3057.822671192607 },
979 { 7.875 , 3921.588652226146 },
980 { 8.0 , 5040.0 },
981 { 8.125 , 6490.864232147346 },
982 { 8.25 , 8376.512350919926 },
983 { 8.375 , 10831.74343769652 },
984 { 8.5 , 14034.40729348341 },
985 { 8.625 , 18219.5209839456 },
986 { 8.75 , 23698.12570174271 },
987 { 8.875 , 30882.5106362809 },
988 { 9.0 , 40320.0 },
989 { 9.125 , 52738.27188619719 },
990 { 9.25 , 69106.22689508938 },
991 { 9.375 , 90715.85129070834 },
992 { 9.5 , 119292.461994609 },
993 { 9.625 , 157143.3684865308 },
994 { 9.75 , 207358.5998902487 },
995 { 9.875 , 274082.281896993 },
996 { 10.0 , 362880.0 },
997 { 10.125 , 481236.7309615494 },
998 { 10.25 , 639232.5987795768 },
999 { 10.375 , 850461.1058503906 },
1000 { 10.5 , 1133278.388948786 },
1001 { 10.625 , 1512504.921682859 },
1002 { 10.75 , 2021746.348929925 },
1003 { 10.875 , 2706562.533732806 },
1004 { 11.0 , 3628800.0 },
1005 { 11.125 , 4872521.900985687 },
1006 { 11.25 , 6552134.137490662 },
1007 { 11.375 , 8823533.973197803 },
1008 { 11.5 , 1.1899423083962249e+7 },
1009 { 11.625 , 1.6070364792880382e+7 },
1010 { 11.75 , 2.1733773250996688e+7 },
1011 { 11.875 , 2.943386755434427e+7 },
1012 { 12.0 , 3.99168e+7 },
1013 { 12.125 , 5.420680614846578e+7 },
1014 { 12.25 , 7.371150904676994e+7 },
1015 { 12.375 , 1.0036769894512501e+8 },
1016 { 12.5 , 1.3684336546556586e+8 },
1017 { 12.625 , 1.8681799071723443e+8 },
1018 { 12.75 , 2.5537183569921107e+8 },
1019 { 12.875 , 3.4952717720783816e+8 },
1020 { 13.0 , 4.790016e+8 },
1021 { 13.125 , 6.572575245501475e+8 },
1022 { 13.25 , 9.029659858229319e+8 },
1023 { 13.375 , 1.2420502744459219e+9 },
1024 { 13.5 , 1.7105420683195732e+9 },
1025 { 13.625 , 2.3585771328050845e+9 },
1026 { 13.75 , 3.2559909051649416e+9 },
1027 { 13.875 , 4.500162406550916e+9 },
1028 { 14.0 , 6.2270208e+9 },
1029 { 14.125 , 8.626505009720685e+9 },
1030 { 14.25 , 1.196429931215385e+10 },
1031 { 14.375 , 1.66124224207142e+10 },
1032 { 14.5 , 2.309231792231424e+10 },
1033 { 14.625 , 3.213561343446927e+10 },
1034 { 14.75 , 4.476987494601794e+10 },
1035 { 14.875 , 6.243975339089396e+10 },
1036 { 15.0 , 8.71782912e+10 },
1037 { 15.125 , 1.218493832623047e+11 },
1038 { 15.25 , 1.704912651981923e+11 },
1039 { 15.375 , 2.388035722977667e+11 },
1040 { 15.5 , 3.348386098735565e+11 },
1041 { 15.625 , 4.699833464791132e+11 },
1042 { 15.75 , 6.603556554537646e+11 },
1043 { 15.875 , 9.287913316895475e+11 },
1044 { 16.0 , 1.307674368e+12 },
1045 { 16.125 , 1.842971921842358e+12 },
1046 { 16.25 , 2.599991794272433e+12 },
1047 { 16.375 , 3.671604924078163e+12 },
1048 { 16.5 , 5.189998453040126e+12 },
1049 { 16.625 , 7.343489788736144e+12 },
1050 { 16.75 , 1.040060157339679e+13 },
1051 { 16.875 , 1.474456239057157e+13 },
1052 { 17.0 , 2.0922789888e+13 },
1053 { 17.125 , 2.971792223970803e+13 },
1054 { 17.25 , 4.224986665692704e+13 },
1055 { 17.375 , 6.012253063177992e+13 },
1056 { 17.5 , 8.563497447516206e+13 },
1057 { 17.625 , 1.220855177377384e+14 },
1058 { 17.75 , 1.742100763543963e+14 },
1059 { 17.875 , 2.488144903408952e+14 },
1060 { 18.0 , 3.55687428096e+14 },
1061 { 18.125 , 5.08919418355e+14 },
1062 { 18.25 , 7.288101998319914e+14 },
1063 { 18.375 , 1.044628969727176e+15 },
1064 { 18.5 , 1.498612053315336e+15 },
1065 { 18.625 , 2.151757250127639e+15 },
1066 { 18.75 , 3.092228855290534e+15 },
1067 { 18.875 , 4.447559014843502e+15 },
1068 { 19.0 , 6.402373705728e+15 },
1069 { 19.125 , 9.224164457684374e+15 },
1070 { 19.25 , 1.330078614693384e+16 },
1071 { 19.375 , 1.919505731873686e+16 },
1072 { 19.5 , 2.772432298633372e+16 },
1073 { 19.625 , 4.007647878362728e+16 },
1074 { 19.75 , 5.797929103669752e+16 },
1075 { 19.875 , 8.39476764051711e+16 },
1076 { 20.0 , 1.21645100408832e+17 },
1077 { 20.5 , 5.406242982335075e+17 },
1078 { 21.0 , 2.43290200817664e+18 },
1079 { 21.5 , 1.10827981137869e+19 },
1080 { 22.0 , 5.109094217170944e+19 },
1081 { 22.5 , 2.382801594464184e+20 },
1082 { 23.0 , 1.124000727777608e+21 },
1083 { 23.5 , 5.361303587544415e+21 },
1084 { 24.0 , 2.585201673888498e+22 },
1085 { 24.5 , 1.259906343072938e+23 },
1086 { 25.0 , 6.204484017332395e+23 },
1087 { 25.5 , 3.086770540528697e+24 },
1088 { 26.0 , 1.551121004333099e+25 },
1089 { 26.5 , 7.871264878348176e+25 },
1090 { 27.0 , 4.032914611266056e+26 },
1091 { 27.5 , 2.085885192762267e+27 },
1092 { 28.0 , 1.088886945041835e+28 },
1093 { 28.5 , 5.736184280096234e+28 },
1094 { 29.0 , 3.048883446117139e+29 },
1095 { 29.5 , 1.634812519827427e+30 },
1096 { 30.0 , 8.841761993739702e+30 },
1097 { 30.5 , 4.822696933490909e+31 },
1098 { 31.0 , 2.65252859812191e+32 },
1099 { 31.5 , 1.470922564714727e+33 },
1100 { 32.0 , 8.222838654177922e+33 },
1101 { 32.5 , 4.633406078851391e+34 },
1102 { 33.0 , 2.631308369336935e+35 },
1103 { 33.5 , 1.505856975626702e+36 },
1104 { 34.0 , 8.683317618811885e+36 },
1105 { 34.5 , 5.044620868349451e+37 },
1106 { 35.0 , 2.952327990396041e+38 },
1107 { 35.5 , 1.740394199580561e+39 },
1108 { 36.0 , 1.033314796638614e+40 },
1109 { 36.5 , 6.178399408510991e+40 },
1110 { 37.0 , 3.719933267899013e+41 },
1111 { 37.5 , 2.255115784106512e+42 },
1112 { 38.0 , 1.376375309122634e+43 },
1113 { 38.5 , 8.456684190399419e+43 },
1114 { 39.0 , 5.230226174666011e+44 },
1115 { 39.5 , 3.255823413303776e+45 },
1116 { 40.0 , 2.039788208119745e+46 },
1117 { 40.5 , 1.286050248254992e+47 },
1118 { 41.0 , 8.159152832478976e+47 },
1119 { 41.5 , 5.208503505432716e+48 },
1120 { 42.0 , 3.345252661316381e+49 },
1121 { 42.5 , 2.161528954754577e+50 },
1122 { 43.0 , 1.40500611775288e+51 },
1123 { 43.5 , 9.186498057706953e+51 },
1124 { 44.0 , 6.041526306337383e+52 },
1125 { 44.5 , 3.996126655102524e+53 },
1126 { 45.0 , 2.658271574788449e+54 },
1127 { 45.5 , 1.778276361520623e+55 },
1128 { 46.0 , 1.196222208654802e+56 },
1129 { 46.5 , 8.091157444918836e+56 },
1130 { 47.0 , 5.502622159812088e+57 },
1131 { 47.5 , 3.762388211887259e+58 },
1132 { 48.0 , 2.586232415111682e+59 },
1133 { 48.5 , 1.787134400646448e+60 },
1134 { 49.0 , 1.241391559253607e+61 },
1135 { 49.5 , 8.667601843135273e+61 },
1136 { 50.0 , 6.082818640342675e+62 },
1137 { 50.5 , 4.290462912351959e+63 },
1138 { 51.0 , 3.041409320171338e+64 },
1139 { 51.5 , 2.16668377073774e+65 },
1140 { 52.0 , 1.551118753287382e+66 },
1141 { 52.5 , 1.115842141929936e+67 },
1142 { 53.0 , 8.065817517094388e+67 },
1143 { 53.5 , 5.858171245132164e+68 },
1144 { 54.0 , 4.274883284060025e+69 },
1145 { 54.5 , 3.134121616145708e+70 },
1146 { 55.0 , 2.308436973392413e+71 },
1147 { 55.5 , 1.70809628079941e+72 },
1148 { 56.0 , 1.269640335365828e+73 },
1149 { 56.5 , 9.479934358436728e+73 },
1150 { 57.0 , 7.109985878048635e+74 },
1151 { 57.5 , 5.356162912516752e+75 },
1152 { 58.0 , 4.052691950487721e+76 },
1153 { 58.5 , 3.079793674697132e+77 },
1154 { 59.0 , 2.350561331282878e+78 },
1155 { 59.5 , 1.801679299697822e+79 },
1156 { 60.0 , 1.386831185456898e+80 },
1157 { 60.5 , 1.071999183320204e+81 },
1158 { 61.0 , 8.320987112741391e+81 },
1159 { 61.5 , 6.485595059087236e+82 },
1160 { 62.0 , 5.075802138772249e+83 },
1161 { 62.5 , 3.98864096133865e+84 },
1162 { 63.0 , 3.146997326038794e+85 },
1163 { 63.5 , 2.492900600836656e+86 },
1164 { 64.0 , 1.98260831540444e+87 },
1165 { 64.5 , 1.582991881531277e+88 },
1166 { 65.0 , 1.268869321858841e+89 },
1167 { 65.5 , 1.021029763587673e+90 },
1168 { 66.0 , 8.247650592082472e+90 },
1169 { 66.5 , 6.687744951499262e+91 },
1170 { 67.0 , 5.443449390774431e+92 },
1171 { 67.5 , 4.447350392747009e+93 },
1172 { 68.0 , 3.647111091818868e+94 },
1173 { 68.5 , 3.001961515104231e+95 },
1174 { 69.0 , 2.48003554243683e+96 },
1175 { 69.5 , 2.056343637846398e+97 },
1176 { 70.0 , 1.711224524281413e+98 },
1177 { 70.5 , 1.429158828303247e+99 },
1178 { 71.0 , 1.19785716699699e+100 },
1179 { 71.5 , 1.00755697395379e+101 },
1180 { 72.0 , 8.50478588567862e+101 },
1181 { 72.5 , 7.20403236376959e+102 },
1182 { 73.0 , 6.12344583768861e+103 },
1183 { 73.5 , 5.22292346373295e+104 },
1184 { 74.0 , 4.47011546151268e+105 },
1185 { 74.5 , 3.83884874584372e+106 },
1186 { 75.0 , 3.30788544151939e+107 },
1187 { 75.5 , 2.85994231565357e+108 },
1188 { 76.0 , 2.48091408113954e+109 },
1189 { 76.5 , 2.15925644831845e+110 },
1190 { 77.0 , 1.88549470166605e+111 },
1191 { 77.5 , 1.65183118296361e+112 },
1192 { 78.0 , 1.45183092028286e+113 },
1193 { 78.5 , 1.2801691667968e+114 },
1194 { 79.0 , 1.13242811782063e+115 },
1195 { 79.5 , 1.00493279593549e+116 },
1196 { 80.0 , 8.94618213078298e+116 },
1197 { 80.5 , 7.98921572768712e+117 },
1198 { 81.0 , 7.15694570462638e+118 },
1199 { 81.5 , 6.43131866078814e+119 },
1200 { 82.0 , 5.79712602074737e+120 },
1201 { 82.5 , 5.24152470854233e+121 },
1202 { 83.0 , 4.75364333701284e+122 },
1203 { 83.5 , 4.32425788454742e+123 },
1204 { 84.0 , 3.94552396972066e+124 },
1205 { 84.5 , 3.6107553335971e+125 },
1206 { 85.0 , 3.31424013456535e+126 },
1207 { 85.5 , 3.05108825688955e+127 },
1208 { 86.0 , 2.81710411438055e+128 },
1209 { 86.5 , 2.60868045964056e+129 },
1210 { 87.0 , 2.42270953836727e+130 },
1211 { 87.5 , 2.25650859758909e+131 },
1212 { 88.0 , 2.10775729837953e+132 },
1213 { 88.5 , 1.97444502289045e+133 },
1214 { 89.0 , 1.85482642257398e+134 },
1215 { 89.5 , 1.74738384525805e+135 },
1216 { 90.0 , 1.65079551609085e+136 },
1217 { 90.5 , 1.56390854150595e+137 },
1218 { 91.0 , 1.48571596448176e+138 },
1219 { 91.5 , 1.41533723006289e+139 },
1220 { 92.0 , 1.3520015276784e+140 },
1221 { 92.5 , 1.29503356550754e+141 },
1222 { 93.0 , 1.24384140546413e+142 },
1223 { 93.5 , 1.19790604809448e+143 },
1224 { 94.0 , 1.15677250708164e+144 },
1225 { 94.5 , 1.12004215496834e+145 },
1226 { 95.0 , 1.08736615665674e+146 },
1227 { 95.5 , 1.05843983644508e+147 },
1228 { 96.0 , 1.03299784882391e+148 },
1229 { 96.5 , 1.01081004380505e+149 },
1230 { 97.0 , 9.9167793487095e+149 },
1231 { 97.5 , 9.75431692271873e+150 },
1232 { 98.0 , 9.61927596824821e+151 },
1233 { 98.5 , 9.51045899965076e+152 },
1234 { 99.0 , 9.42689044888325e+153 },
1235 { 99.5 , 9.367802114656e+154 },
1236 { 100.0 , 9.33262154439441e+155 },
1237 };
1238
1239 @Test
1240 void testGamma() {
1241
1242 for (int i = 0; i < GAMMA_REF.length; i++) {
1243 final double[] ref = GAMMA_REF[i];
1244 final double x = ref[0];
1245 final double expected = ref[1];
1246 final double actual = Gamma.gamma(x);
1247 final double absX = FastMath.abs(x);
1248 final int ulps;
1249 if (absX <= 8.0) {
1250 ulps = 3;
1251 } else if (absX <= 20.0) {
1252 ulps = 5;
1253 } else if (absX <= 30.0) {
1254 ulps = 50;
1255 } else if (absX <= 50.0) {
1256 ulps = 180;
1257 } else {
1258 ulps = 500;
1259 }
1260 final double tol = ulps * FastMath.ulp(expected);
1261 assertEquals(expected, actual, tol, Double.toString(x));
1262 }
1263 }
1264
1265 @Test
1266 void testGammaField() {
1267
1268 for (int i = 0; i < GAMMA_REF.length; i++) {
1269 final double[] ref = GAMMA_REF[i];
1270 final Binary64 x = new Binary64(ref[0]);
1271 final double expected = ref[1];
1272 final double actual = Gamma.gamma(x).getReal();
1273 final Binary64 absX = x.abs();
1274 final int ulps;
1275 if (absX.getReal() <= 8.0) {
1276 ulps = 3;
1277 } else if (absX.getReal() <= 20.0) {
1278 ulps = 5;
1279 } else if (absX.getReal() <= 30.0) {
1280 ulps = 50;
1281 } else if (absX.getReal() <= 50.0) {
1282 ulps = 180;
1283 } else {
1284 ulps = 500;
1285 }
1286 final double tol = ulps * FastMath.ulp(expected);
1287 assertEquals(expected, actual, tol, Double.toString(x.getReal()));
1288 }
1289 }
1290
1291 @Test
1292 void testGammaNegativeInteger() {
1293
1294 for (int i = -100; i <= 0; i++) {
1295 assertTrue(Double.isNaN(Gamma.gamma(i)), Integer.toString(i));
1296 }
1297 }
1298
1299 @Test
1300 void testGammaNegativeIntegerField() {
1301
1302 for (int i = -100; i <= 0; i++) {
1303 assertTrue(Gamma.gamma(new Binary64(i)).isNaN(), Integer.toString(i));
1304 }
1305 }
1306
1307 @Test
1308 void testGammaNegativeDouble() {
1309
1310
1311
1312 double previousGamma = Gamma.gamma(-18.5);
1313 for (double x = -19.5; x > -25; x -= 1.0) {
1314 double gamma = Gamma.gamma(x);
1315 assertEquals( (int) FastMath.signum(previousGamma),
1316 - (int) FastMath.signum(gamma));
1317
1318 previousGamma = gamma;
1319 }
1320 }
1321
1322 @Test
1323 void testGammaNegativeDoubleField() {
1324
1325
1326
1327 Binary64 previousGamma = Gamma.gamma(new Binary64(-18.5));
1328 for (Binary64 x = new Binary64(-19.5); x.getReal() > -25; x = x.subtract(1.0)) {
1329 Binary64 gamma = Gamma.gamma(x);
1330 assertEquals( (int) FastMath.signum(previousGamma.getReal()),
1331 - (int) FastMath.signum(gamma.getReal()));
1332
1333 previousGamma = gamma;
1334 }
1335 }
1336
1337 private void checkRelativeError(String msg, double expected, double actual,
1338 double tolerance) {
1339
1340 assertEquals(expected, actual, FastMath.abs(tolerance * actual), msg);
1341 }
1342 }