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 }