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