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