View Javadoc
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.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         // values for negative powers of 10 from 1 to 30 as computed by webMathematica with 20 digits
243         // see functions.wolfram.com
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         // values for negative powers of 10 from 1 to 30 as computed by webMathematica with 20 digits
258         // see functions.wolfram.com
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         // computed using webMathematica.  For example, to compute trigamma($i) = Polygamma(1, $i), use
288         //
289         // http://functions.wolfram.com/webMathematica/Evaluated.jsp?name=PolyGamma2&plottype=0&vars={%221%22,%22$i%22}&digits=20
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         // computed using webMathematica.  For example, to compute trigamma($i) = Polygamma(1, $i), use
315         //
316         // http://functions.wolfram.com/webMathematica/Evaluated.jsp?name=PolyGamma2&plottype=0&vars={%221%22,%22$i%22}&digits=20
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      * Reference data for the {@link Gamma#logGamma(double)} function. This data
369      * was generated with the following <a
370      * href="http://maxima.sourceforge.net/">Maxima</a> script.
371      *
372      * <pre>
373      * kill(all);
374      *
375      * fpprec : 64;
376      * gamln(x) := log(gamma(x));
377      * x : append(makelist(bfloat(i / 8), i, 1, 80),
378      *     [0.8b0, 1b2, 1b3, 1b4, 1b5, 1b6, 1b7, 1b8, 1b9, 1b10]);
379      *
380      * for i : 1 while i <= length(x) do
381      *     print("{", float(x[i]), ",", float(gamln(x[i])), "},");
382      * </pre>
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      * <p>
577      * Reference values for the {@link Gamma#invGamma1pm1(double)} method.
578      * These values were generated with the following <a
579      * href="http://maxima.sourceforge.net/">Maxima</a> script
580      * </p>
581      *
582      * <pre>
583      * kill(all);
584      *
585      * fpprec : 64;
586      * gam1(x) := 1 / gamma(1 + x) - 1;
587      * x : makelist(bfloat(i / 8), i, -4, 12);
588      *
589      * for i : 1 while i <= length(x) do print("{",
590      *                                         float(x[i]),
591      *                                         ",",
592      *                                         float(gam1(x[i])),
593      *                                         "},");
594      * </pre>
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      * Reference data for the {@link Gamma#gamma(double)} function. This
758      * data was generated with the following <a
759      * href="http://maxima.sourceforge.net/">Maxima</a> script.
760      *
761      * <pre>
762      * kill(all);
763      *
764      * fpprec : 64;
765      *
766      * EPSILON : 10**(-fpprec + 1);
767      * isInteger(x) := abs(x - floor(x)) <= EPSILON * abs(x);
768      *
769      * x : makelist(bfloat(i / 8), i, -160, 160);
770      * x : append(x, makelist(bfloat(i / 2), i, 41, 200));
771      *
772      * for i : 1 while i <= length(x) do if not(isInteger(x[i])) then
773      *     print("{", float(x[i]), ",", float(gamma(x[i])), "},");
774      * </pre>
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         // check that the gamma function properly switches sign
1310         // see: https://en.wikipedia.org/wiki/Gamma_function
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         // check that the gamma function properly switches sign
1325         // see: https://en.wikipedia.org/wiki/Gamma_function
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 }