xref: /netbsd-src/external/lgpl3/mpfr/dist/tests/tgamma_inc.c (revision cb63e24e8d6aae7ddac1859a9015f48b1d8bd90e)
1 /* Test file for mpfr_gamma_inc
2 
3 Copyright 2016-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5 
6 This file is part of the GNU MPFR Library.
7 
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16 License for more details.
17 
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22 
23 #include "mpfr-test.h"
24 
25 #define TEST_FUNCTION mpfr_gamma_inc
26 #define TWO_ARGS
27 #define TEST_RANDOM_POS2 0 /* the 2nd argument is never negative */
28 #define TGENERIC_NOWARNING 1
29 #define TEST_RANDOM_EMAX 8
30 #define TEST_RANDOM_EMIN -32
31 #define REDUCE_EMAX TEST_RANDOM_EMAX
32 #define REDUCE_EMIN TEST_RANDOM_EMIN
33 #include "tgeneric.c"
34 
35 /* do k random tests at precision p */
36 static void
37 test_random (mpfr_prec_t p, int k)
38 {
39   mpfr_t a, x, y, z, t;
40 
41   mpfr_inits2 (p, a, x, y, z, (mpfr_ptr) 0);
42   mpfr_init2 (t, p + 20);
43   while (k--)
44     {
45       do mpfr_urandomb (a, RANDS); while (mpfr_zero_p (a));
46       if (RAND_BOOL ())
47         mpfr_neg (a, a, MPFR_RNDN);
48       do mpfr_urandomb (x, RANDS); while (mpfr_zero_p (x));
49       mpfr_gamma_inc (y, a, x, MPFR_RNDN);
50       mpfr_gamma_inc (t, a, x, MPFR_RNDN);
51       if (mpfr_can_round (t, mpfr_get_prec (z), MPFR_RNDN, MPFR_RNDN, p))
52         {
53           mpfr_set (z, t, MPFR_RNDN);
54           if (mpfr_cmp (y, z))
55             {
56               printf ("mpfr_gamma_inc failed for a=");
57               mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
58               printf (" x=");
59               mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
60               printf ("\nexpected ");
61               mpfr_out_str (stdout, 10, 0, z, MPFR_RNDN);
62               printf ("\ngot      ");
63               mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN);
64               printf ("\n");
65               exit (1);
66             }
67         }
68     }
69   mpfr_clears (a, x, y, z, (mpfr_ptr) 0);
70   mpfr_clear (t);
71 }
72 
73 static void
74 specials (void)
75 {
76   mpfr_t a, x;
77   int inex;
78 
79   mpfr_init2 (a, 2);
80   mpfr_init2 (x, 2);
81 
82   mpfr_set_inf (a, 1);
83   mpfr_set_inf (x, 1);
84   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
85   MPFR_ASSERTN (mpfr_nan_p (a));
86 
87   mpfr_set_inf (a, 1);
88   mpfr_set_inf (x, -1);
89   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
90   MPFR_ASSERTN (mpfr_nan_p (a));
91 
92   mpfr_set_inf (a, -1);
93   mpfr_set_inf (x, -1);
94   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
95   MPFR_ASSERTN (mpfr_nan_p (a));
96 
97   mpfr_set_inf (a, -1);
98   mpfr_set_inf (x, 1);
99   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
100   MPFR_ASSERTN (mpfr_zero_p (a) || mpfr_signbit (a) == 0);
101 
102   mpfr_set_inf (a, 1);
103   mpfr_set_ui (x, 1, MPFR_RNDN);
104   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
105   MPFR_ASSERTN (mpfr_inf_p (a) || mpfr_signbit (a) == 0);
106 
107   mpfr_set_inf (a, -1);
108   mpfr_set_ui (x, 2, MPFR_RNDN);
109   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
110   MPFR_ASSERTN (mpfr_zero_p (a) || mpfr_signbit (a) == 0);
111 
112   mpfr_set_inf (a, -1);
113   mpfr_set_ui (x, 1, MPFR_RNDN);
114   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
115   MPFR_ASSERTN (mpfr_zero_p (a) || mpfr_signbit (a) == 0);
116 
117   mpfr_set_inf (a, -1);
118   mpfr_set_ui_2exp (x, 1, -1, MPFR_RNDN);
119   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
120   MPFR_ASSERTN (mpfr_inf_p (a) || mpfr_signbit (a) == 0);
121 
122   mpfr_set_ui (a, 1, MPFR_RNDN);
123   mpfr_set_inf (x, 1);
124   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
125   MPFR_ASSERTN (mpfr_zero_p (a) || mpfr_signbit (a) == 0);
126 
127   /* gamma_inc(1,-x) = exp(x) tends to +Inf */
128   mpfr_set_ui (a, 1, MPFR_RNDN);
129   mpfr_set_inf (x, -1);
130   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
131   MPFR_ASSERTN (mpfr_inf_p (a) || mpfr_signbit (a) == 0);
132 
133   /* gamma_inc(0,x) for x < 0 has imaginary part -Pi and thus gives NaN
134      over the reals */
135   mpfr_set_zero (a, 1);
136   mpfr_set_inf (x, -1);
137   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
138   MPFR_ASSERTN (mpfr_nan_p (a));
139 
140   /* gamma_inc(-1,x) for x < 0 has imaginary part +Pi and thus gives NaN */
141   mpfr_set_si (a, -1, MPFR_RNDN);
142   mpfr_set_inf (x, -1);
143   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
144   MPFR_ASSERTN (mpfr_nan_p (a));
145 
146   /* gamma_inc(-2,x) for x < 0 has imaginary part -Pi/2 and thus gives NaN */
147   mpfr_set_si (a, -2, MPFR_RNDN);
148   mpfr_set_inf (x, -1);
149   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
150   MPFR_ASSERTN (mpfr_nan_p (a));
151 
152   mpfr_set_ui_2exp (a, 1, -1, MPFR_RNDN);
153   mpfr_set_inf (x, -1);
154   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
155   MPFR_ASSERTN (mpfr_nan_p (a));
156 
157   mpfr_set_si_2exp (a, -1, -1, MPFR_RNDN);
158   mpfr_set_inf (x, -1);
159   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
160   MPFR_ASSERTN (mpfr_nan_p (a));
161 
162   /* gamma_inc(0,x) = -Ei(-x) */
163   mpfr_set_zero (a, 1);
164   mpfr_set_si (x, -1, MPFR_RNDN);
165   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
166   MPFR_ASSERTN (mpfr_nan_p (a));
167 
168   /* gamma_inc(a,0) = gamma(a) thus gamma_inc(-Inf,0) = gamma(-Inf) = NaN */
169   mpfr_set_inf (a, -1);
170   mpfr_set_zero (x, 1);
171   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
172   MPFR_ASSERTN (mpfr_nan_p (a));
173 
174   mpfr_set_inf (a, -1);
175   mpfr_set_si (x, -1, MPFR_RNDN);
176   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
177   MPFR_ASSERTN (mpfr_nan_p (a));
178 
179   /* check gamma_inc(2,0) = 1 is exact */
180   mpfr_set_ui (a, 2, MPFR_RNDN);
181   mpfr_set_ui (x, 0, MPFR_RNDN);
182   mpfr_clear_inexflag ();
183   inex = mpfr_gamma_inc (a, a, x, MPFR_RNDN);
184   if (mpfr_cmp_ui (a, 1))
185     {
186       printf ("Error for gamma_inc(2,0)\n");
187       printf ("expected 1\n");
188       printf ("got      ");
189       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
190       printf ("\n");
191       exit (1);
192     }
193   if (inex != 0)
194     {
195       printf ("Wrong ternary value for gamma_inc(2,0)\n");
196       printf ("expected 0\n");
197       printf ("got      %d\n", inex);
198       exit (1);
199     }
200   if (mpfr_inexflag_p ())
201     {
202       printf ("Wrong inexact flag for gamma_inc(2,0)\n");
203       printf ("expected 0\n");
204       printf ("got      1\n");
205       exit (1);
206     }
207 
208   /* check gamma_inc(0,1) = 0.219383934395520 */
209   mpfr_set_ui (a, 0, MPFR_RNDN);
210   mpfr_set_ui (x, 1, MPFR_RNDN);
211   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
212   if (mpfr_cmp_ui_2exp (a, 1, -2))
213     {
214       printf ("Error for gamma_inc(0,1)\n");
215       printf ("expected 0.25\n");
216       printf ("got      ");
217       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
218       printf ("\n");
219       exit (1);
220     }
221 
222   /* check gamma_inc(-1,1) = 0.148495506775922 */
223   mpfr_set_si (a, -1, MPFR_RNDN);
224   mpfr_set_ui (x, 1, MPFR_RNDN);
225   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
226   if (mpfr_cmp_ui_2exp (a, 1, -3))
227     {
228       printf ("Error for gamma_inc(-1,1)\n");
229       printf ("expected 0.125\n");
230       printf ("got      ");
231       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
232       printf ("\n");
233       exit (1);
234     }
235 
236   /* check gamma_inc(-2,1) = 0.109691967197760 */
237   mpfr_set_si (a, -2, MPFR_RNDN);
238   mpfr_set_ui (x, 1, MPFR_RNDN);
239   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
240   if (mpfr_cmp_ui_2exp (a, 1, -3))
241     {
242       printf ("Error for gamma_inc(-2,1)\n");
243       printf ("expected 0.125\n");
244       printf ("got      ");
245       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
246       printf ("\n");
247       exit (1);
248     }
249 
250   /* check gamma_inc(-3,1) = 0.109691967197760 */
251   mpfr_set_si (a, -3, MPFR_RNDN);
252   mpfr_set_ui (x, 1, MPFR_RNDN);
253   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
254   if (mpfr_cmp_ui_2exp (a, 3, -5))
255     {
256       printf ("Error for gamma_inc(-3,1)\n");
257       printf ("expected 3/32\n");
258       printf ("got      ");
259       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
260       printf ("\n");
261       exit (1);
262     }
263 
264   /* check gamma_inc(-100,1) = 0.00364201018241046 */
265   mpfr_set_si (a, -100, MPFR_RNDN);
266   mpfr_set_ui (x, 1, MPFR_RNDN);
267   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
268   if (mpfr_cmp_ui_2exp (a, 1, -8))
269     {
270       printf ("Error for gamma_inc(-100,1)\n");
271       printf ("expected 1/256\n");
272       printf ("got      ");
273       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
274       printf ("\n");
275       exit (1);
276     }
277 
278   /* TODO: Once internal overflow is supported, add new tests with
279      larger exponents, e.g. 64 (in addition to 25). */
280   mpfr_set_prec (a, 1);
281   mpfr_set_prec (x, 1);
282   mpfr_set_ui_2exp (a, 1, 25, MPFR_RNDN);
283   mpfr_set_ui_2exp (x, 1, -25, MPFR_RNDN);
284   mpfr_gamma_inc (a, a, x, MPFR_RNDN);
285 
286   mpfr_clear (a);
287   mpfr_clear (x);
288 }
289 
290 /* tests for negative integer a: for -n <= a <= -1, perform k tests
291    with random x in 0..|a| and precision 'prec' */
292 static void
293 test_negint (long n, unsigned long k, mpfr_prec_t prec)
294 {
295   long i, j;
296   mpfr_t a, x, y;
297 
298   mpfr_init2 (a, prec);
299   mpfr_init2 (x, prec);
300   mpfr_init2 (y, prec);
301   for (i = 1; i <= n; i++)
302     {
303       mpfr_set_si (a, -i, MPFR_RNDN);
304       for (j = 1; j <= k; j++)
305         {
306           mpfr_urandomb (x, RANDS);
307           mpfr_mul_si (x, x, j, MPFR_RNDN);
308           mpfr_set_prec (y, prec + 20);
309           mpfr_gamma_inc (y, a, x, MPFR_RNDZ);
310           mpfr_gamma_inc (x, a, x, MPFR_RNDZ);
311           mpfr_prec_round (y, prec, MPFR_RNDZ);
312           if (mpfr_cmp (x, y))
313             {
314               printf ("Error in mpfr_gamma_inc(%ld,%ld) with MPFR_RNDZ\n",
315                       -i, j);
316               printf ("expected ");
317               mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN);
318               printf ("\ngot      ");
319               mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
320               printf ("\n");
321               exit (1);
322             }
323         }
324     }
325   mpfr_clear (a);
326   mpfr_clear (x);
327   mpfr_clear (y);
328 }
329 
330 static void
331 coverage (void)
332 {
333   mpfr_t a, x, y;
334   int inex;
335 
336   mpfr_init2 (a, MPFR_PREC_MIN);
337   mpfr_init2 (x, MPFR_PREC_MIN);
338   mpfr_init2 (y, MPFR_PREC_MIN);
339 
340   /* exercise test MPFR_GET_EXP(a) + 3 > w in mpfr_gamma_inc_negint */
341   mpfr_set_si (a, -256, MPFR_RNDN);
342   mpfr_set_ui (x, 1, MPFR_RNDN);
343   inex = mpfr_gamma_inc (y, a, x, MPFR_RNDN);
344   /* gamma_inc(-256,1) = 0.00143141575826821 is rounded to 2^(-10) */
345   MPFR_ASSERTN(inex < 0);
346   MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 1, -10) == 0);
347 
348   /* exercise the case MPFR_IS_ZERO(s) in mpfr_gamma_inc_negint */
349   mpfr_set_prec (a, 4);
350   mpfr_set_prec (x, 4);
351   mpfr_set_prec (y, 4);
352   inex = mpfr_set_si (a, -15, MPFR_RNDN);
353   MPFR_ASSERTN (inex == 0);
354   inex = mpfr_set_ui (x, 15, MPFR_RNDN);
355   MPFR_ASSERTN (inex == 0);
356   /* gamma_inc(-15,15) = 0.2290433491e-25, rounded to 14*2^(-89) */
357   inex = mpfr_gamma_inc (y, a, x, MPFR_RNDN);
358   MPFR_ASSERTN(inex < 0);
359   MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 14, -89) == 0);
360 
361   mpfr_clear (a);
362   mpfr_clear (x);
363   mpfr_clear (y);
364 }
365 
366 int
367 main (int argc, char *argv[])
368 {
369   mpfr_prec_t p;
370 
371   tests_start_mpfr ();
372 
373   if (argc == 4) /* tgamma_inc a x prec: print gamma_inc(a,x) to prec bits */
374     {
375       mpfr_prec_t p = atoi (argv[3]);
376       mpfr_t a, x;
377       mpfr_init2 (a, p);
378       mpfr_init2 (x, p);
379       mpfr_set_str (a, argv[1], 10, MPFR_RNDN);
380       mpfr_set_str (x, argv[2], 10, MPFR_RNDN);
381       mpfr_gamma_inc (x, a, x, MPFR_RNDN);
382       mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN);
383       printf ("\n");
384       mpfr_clear (a);
385       mpfr_clear (x);
386       return 0;
387     }
388 
389   coverage ();
390 
391   specials ();
392 
393   test_negint (30, 10, 53);
394 
395   for (p = MPFR_PREC_MIN; p < 100; p++)
396     test_random (p, 10);
397 
398   test_generic (MPFR_PREC_MIN, 100, 5);
399 
400   tests_end_mpfr ();
401   return 0;
402 }
403