1/* Test file for mpfr_gamma_inc 2 3Copyright 2016-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 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 */ 36static void 37test_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 73static void 74specials (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' */ 292static void 293test_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 330static void 331coverage (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 366int 367main (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