1 /* Test file for mpfr_add1sp. 2 3 Copyright 2004-2020 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 static void check_special (void); 26 static void check_random (mpfr_prec_t p); 27 28 static int 29 mpfr_add_cf (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t r) 30 { 31 mpfr_clear_flags (); /* allows better checking */ 32 return mpfr_add (a, b, c, r); 33 } 34 35 static void 36 check_overflow (void) 37 { 38 mpfr_t x, y, z1, z2; 39 mpfr_exp_t emin, emax; 40 41 emin = mpfr_get_emin (); 42 emax = mpfr_get_emax (); 43 44 set_emin (-1021); 45 set_emax (1024); 46 47 mpfr_inits (x, y, z1, z2, (mpfr_ptr) 0); 48 49 mpfr_set_str1 (x, "8.00468257869324898448e+307"); 50 mpfr_set_str1 (y, "7.44784712422708645156e+307"); 51 mpfr_add1sp (z1, x, y, MPFR_RNDN); 52 mpfr_add1 (z2, x, y, MPFR_RNDN); 53 if (mpfr_cmp (z1, z2)) 54 { 55 printf ("Overflow bug in add1sp.\n"); 56 exit (1); 57 } 58 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0); 59 60 set_emin (emin); 61 set_emax (emax); 62 } 63 64 static void 65 bug20171217 (void) 66 { 67 mpfr_t a, b, c; 68 69 mpfr_init2 (a, 137); 70 mpfr_init2 (b, 137); 71 mpfr_init2 (c, 137); 72 mpfr_set_str_binary (b, "0.11111111111111111111111111111111111111111111111111111111111111111111000000000000000000000000000000000000000000000000000000000000000000000E-66"); 73 mpfr_set_str_binary (c, "0.11111111111111111111111111111111111111111111111111111111111111111000000000000000000000000000000000000000000000000000000000000000000110000E-2"); 74 mpfr_add_cf (a, b, c, MPFR_RNDN); 75 mpfr_set_str_binary (b, "0.10000000000000000000000000000000000000000000000000000000000000000100000000000000000000000000000000000000000000000000000000000000000001000E-1"); 76 MPFR_ASSERTN(mpfr_equal_p (a, b)); 77 mpfr_clear (a); 78 mpfr_clear (b); 79 mpfr_clear (c); 80 } 81 82 static void 83 bug20190903 (void) 84 { 85 mpfr_t a, b, c, d; 86 int inex; 87 mpfr_flags_t flags; 88 89 /* Bug in r13574, fixed in r13578. 90 Note: to reproduce the failure, GMP_NUMB_BITS == 64 is assumed. */ 91 mpfr_inits2 (128, a, b, c, d, (mpfr_ptr) 0); 92 mpfr_set_str_binary (b, "0.11111111111111111100000000000000000001001111111101111111110000101001111100111110110010011001111110000000101001001001110110101110E0"); 93 mpfr_set_str_binary (c, "0.10000001011101000010000111100111011110100000001001000010011011001000110100111111101100001101001101011101100100011000000101110111E-126"); 94 mpfr_add_cf (a, b, c, MPFR_RNDN); 95 mpfr_set_str_binary (b, "0.11111111111111111100000000000000000001001111111101111111110000101001111100111110110010011001111110000000101001001001110110110000E0"); 96 MPFR_ASSERTN (mpfr_equal_p (a, b)); 97 mpfr_clears (a, b, c, d, (mpfr_ptr) 0); 98 99 /* Bug in r13574, fixed in r13586. */ 100 /* Figure with GMP_NUMB_BITS = 4: 101 b = 1111 1000 102 c = 1000 0001 103 */ 104 mpfr_inits2 (2 * GMP_NUMB_BITS, a, b, c, d, (mpfr_ptr) 0); 105 mpfr_set_ui_2exp (d, 1, 3 * GMP_NUMB_BITS, MPFR_RNDN); 106 mpfr_set_ui_2exp (c, 1, 2 * GMP_NUMB_BITS - 1, MPFR_RNDN); 107 mpfr_sub (b, d, c, MPFR_RNDN); 108 mpfr_add_ui (c, c, 1, MPFR_RNDN); 109 inex = mpfr_add_cf (a, b, c, MPFR_RNDN); 110 flags = __gmpfr_flags; 111 MPFR_ASSERTN (mpfr_equal_p (a, d)); 112 MPFR_ASSERTN (inex < 0); 113 MPFR_ASSERTN (flags == MPFR_FLAGS_INEXACT); 114 inex = mpfr_add_cf (a, b, c, MPFR_RNDU); 115 flags = __gmpfr_flags; 116 mpfr_add_ui (d, d, 1, MPFR_RNDU); 117 MPFR_ASSERTN (mpfr_equal_p (a, d)); 118 MPFR_ASSERTN (inex > 0); 119 MPFR_ASSERTN (flags == MPFR_FLAGS_INEXACT); 120 mpfr_clears (a, b, c, d, (mpfr_ptr) 0); 121 } 122 123 /* Check corner case b = 1, c = 2^(-p) for MPFR_PREC_MIN <= p <= pmax. 124 With RNDN, result is 1, except for p=1, where it is 2. */ 125 static void 126 test_corner_1 (mpfr_prec_t pmax) 127 { 128 mpfr_prec_t p; 129 130 for (p = MPFR_PREC_MIN; p <= pmax; p++) 131 { 132 mpfr_t a, b, c; 133 int inex; 134 mpfr_init2 (a, p); 135 mpfr_init2 (b, p); 136 mpfr_init2 (c, p); 137 mpfr_set_ui (b, 1, MPFR_RNDN); 138 mpfr_set_ui_2exp (c, 1, -p, MPFR_RNDN); 139 inex = mpfr_add_cf (a, b, c, MPFR_RNDN); 140 if (p == 1) /* special case, since 2^(p-1) is odd */ 141 { 142 MPFR_ASSERTN(inex > 0); 143 MPFR_ASSERTN(mpfr_cmp_ui (a, 2) == 0); 144 } 145 else 146 { 147 MPFR_ASSERTN(inex < 0); 148 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0); 149 } 150 mpfr_clear (a); 151 mpfr_clear (b); 152 mpfr_clear (c); 153 } 154 } 155 156 static void 157 coverage (void) 158 { 159 mpfr_t a, b, c; 160 int inex; 161 mpfr_exp_t emax; 162 mpfr_prec_t p; 163 164 mpfr_init (a); 165 mpfr_init (b); 166 mpfr_init (c); 167 168 /* coverage test in mpfr_add1sp: case round away, where add_one_ulp 169 gives a carry, and the new exponent is below emax */ 170 for (p = MPFR_PREC_MIN; p <= 3 * GMP_NUMB_BITS; p++) 171 { 172 mpfr_set_prec (a, p); 173 mpfr_set_prec (b, p); 174 mpfr_set_prec (c, p); 175 mpfr_set_ui (b, 1, MPFR_RNDN); 176 mpfr_nextbelow (b); /* b = 1 - 2^(-p) (including for p=1) */ 177 mpfr_set_ui_2exp (c, 1, -p-1, MPFR_RNDN); 178 /* c = 2^(-p-1) thus b+c = 1 - 2^(-p-1) should be rounded to 1 */ 179 inex = mpfr_add_cf (a, b, c, MPFR_RNDU); 180 MPFR_ASSERTN(inex > 0); 181 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0); 182 } 183 184 /* coverage test in mpfr_add1sp2: case GMP_NUMB_BITS <= d < 2*GMP_NUMB_BITS 185 and a1 = 0 */ 186 mpfr_set_prec (a, GMP_NUMB_BITS + 2); 187 mpfr_set_prec (b, GMP_NUMB_BITS + 2); 188 mpfr_set_prec (c, GMP_NUMB_BITS + 2); 189 mpfr_set_ui (b, 1, MPFR_RNDN); 190 mpfr_nextbelow (b); /* b = 1 - 2^(-p) with p = GMP_NUMB_BITS+2 */ 191 mpfr_set_ui_2exp (c, 1, -GMP_NUMB_BITS-1, MPFR_RNDN); 192 mpfr_nextbelow (c); /* c = 2^(1-p) - 2^(1-2p) */ 193 /* a = 1 + 2^(-p) - 2^(1-2p) should be rounded to 1 with RNDN */ 194 inex = mpfr_add_cf (a, b, c, MPFR_RNDN); 195 MPFR_ASSERTN(inex < 0); 196 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0); 197 198 /* coverage test in mpfr_add1sp2: case round away, where add_one_ulp 199 gives a carry, and the new exponent is below emax */ 200 mpfr_set_prec (a, GMP_NUMB_BITS + 1); 201 mpfr_set_prec (b, GMP_NUMB_BITS + 1); 202 mpfr_set_prec (c, GMP_NUMB_BITS + 1); 203 mpfr_set_ui (b, 1, MPFR_RNDN); 204 mpfr_nextbelow (b); /* b = 1 - 2^(-p) */ 205 mpfr_set_ui_2exp (c, 1, -GMP_NUMB_BITS-2, MPFR_RNDN); 206 /* c = 2^(-p-1) */ 207 inex = mpfr_add_cf (a, b, c, MPFR_RNDU); 208 MPFR_ASSERTN(inex > 0); 209 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0); 210 211 /* coverage test in mpfr_add1sp3: case GMP_NUMB_BITS <= d < 2*GMP_NUMB_BITS 212 and a2 == 0 */ 213 mpfr_set_prec (a, 2 * GMP_NUMB_BITS + 2); 214 mpfr_set_prec (b, 2 * GMP_NUMB_BITS + 2); 215 mpfr_set_prec (c, 2 * GMP_NUMB_BITS + 2); 216 mpfr_set_ui (b, 1, MPFR_RNDN); 217 mpfr_nextbelow (b); /* b = 1 - 2^(-p) with p = 2*GMP_NUMB_BITS+2 */ 218 mpfr_set_ui_2exp (c, 1, -2*GMP_NUMB_BITS-1, MPFR_RNDN); 219 mpfr_nextbelow (c); /* c = 2^(1-p) - 2^(1-2p) */ 220 /* a = 1 + 2^(-p) - 2^(1-2p) should be rounded to 1 with RNDN */ 221 inex = mpfr_add_cf (a, b, c, MPFR_RNDN); 222 MPFR_ASSERTN(inex < 0); 223 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0); 224 225 /* coverage test in mpfr_add1sp3: case bx > emax */ 226 emax = mpfr_get_emax (); 227 mpfr_set_emax (1); 228 mpfr_set_prec (a, 2 * GMP_NUMB_BITS + 1); 229 mpfr_set_prec (b, 2 * GMP_NUMB_BITS + 1); 230 mpfr_set_prec (c, 2 * GMP_NUMB_BITS + 1); 231 mpfr_set_ui_2exp (b, 1, mpfr_get_emax () - 1, MPFR_RNDN); 232 mpfr_nextbelow (b); 233 mpfr_mul_2ui (b, b, 1, MPFR_RNDN); 234 /* now b is the largest number < +Inf */ 235 mpfr_div_2ui (c, b, GMP_NUMB_BITS - 1, MPFR_RNDN); 236 /* we are in the case d < GMP_NUMB_BITS of mpfr_add1sp3 */ 237 inex = mpfr_add_cf (a, b, b, MPFR_RNDU); 238 MPFR_ASSERTN(inex > 0); 239 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0); 240 mpfr_set_emax (emax); 241 242 /* coverage test in mpfr_add1sp3: case round away, where add_one_ulp gives 243 a carry, no overflow */ 244 mpfr_set_prec (a, 2 * GMP_NUMB_BITS + 1); 245 mpfr_set_prec (b, 2 * GMP_NUMB_BITS + 1); 246 mpfr_set_prec (c, 2 * GMP_NUMB_BITS + 1); 247 mpfr_set_ui (b, 1, MPFR_RNDN); 248 mpfr_nextbelow (b); /* b = 1 - 2^(-p) */ 249 mpfr_set_ui_2exp (c, 1, -2 * GMP_NUMB_BITS - 2, MPFR_RNDN); 250 /* c = 2^(-p-1) */ 251 inex = mpfr_add_cf (a, b, c, MPFR_RNDU); 252 MPFR_ASSERTN(inex > 0); 253 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0); 254 255 /* coverage test in mpfr_add1sp3: case round away, where add_one_ulp gives 256 a carry, with overflow */ 257 emax = mpfr_get_emax (); 258 mpfr_set_emax (1); 259 mpfr_set_prec (a, 2 * GMP_NUMB_BITS + 1); 260 mpfr_set_prec (b, 2 * GMP_NUMB_BITS + 1); 261 mpfr_set_prec (c, 2 * GMP_NUMB_BITS + 1); 262 mpfr_set_ui_2exp (b, 1, mpfr_get_emax () - 1, MPFR_RNDN); 263 mpfr_nextbelow (b); 264 mpfr_mul_2ui (b, b, 1, MPFR_RNDN); 265 /* now b is the largest number < +Inf */ 266 mpfr_set_ui_2exp (c, 1, mpfr_get_emin () - 1, MPFR_RNDN); 267 inex = mpfr_add_cf (a, b, c, MPFR_RNDU); 268 MPFR_ASSERTN(inex > 0); 269 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0); 270 mpfr_set_emax (emax); 271 272 mpfr_clear (a); 273 mpfr_clear (b); 274 mpfr_clear (c); 275 } 276 277 int 278 main (void) 279 { 280 mpfr_prec_t p; 281 int i; 282 283 tests_start_mpfr (); 284 285 coverage (); 286 test_corner_1 (1024); 287 bug20171217 (); 288 bug20190903 (); 289 check_special (); 290 for (p = MPFR_PREC_MIN; p < 200; p++) 291 check_random (p); 292 for (i = 0; i < 200; i++) 293 { 294 /* special precisions */ 295 check_random (GMP_NUMB_BITS); 296 check_random (2 * GMP_NUMB_BITS); 297 } 298 check_overflow (); 299 300 tests_end_mpfr (); 301 return 0; 302 } 303 304 #define STD_ERROR \ 305 do \ 306 { \ 307 printf("ERROR: for %s and p=%lu and i=%d:\nB=", \ 308 mpfr_print_rnd_mode ((mpfr_rnd_t) r), (unsigned long) p, i); \ 309 mpfr_dump (b); \ 310 printf ("C="); mpfr_dump (c); \ 311 printf ("add1 : "); mpfr_dump (a1); \ 312 printf ("add1sp: "); mpfr_dump (a2); \ 313 exit(1); \ 314 } \ 315 while (0) 316 317 #define STD_ERROR2 \ 318 do \ 319 { \ 320 printf("ERROR: Wrong inexact flag for %s and p=%lu and i=%d:\nB=", \ 321 mpfr_print_rnd_mode ((mpfr_rnd_t) r), (unsigned long) p, i); \ 322 mpfr_dump (b); \ 323 printf ("C="); mpfr_dump (c); \ 324 printf ("A="); mpfr_dump (a1); \ 325 printf ("Add1: %d. Add1sp: %d\n", inexact1, inexact2); \ 326 exit(1); \ 327 } \ 328 while (0) 329 330 #define SET_PREC(_p) \ 331 do \ 332 { \ 333 p = _p; \ 334 mpfr_set_prec(a1, _p); mpfr_set_prec(a2, _p); \ 335 mpfr_set_prec(b, _p); mpfr_set_prec(c, _p); \ 336 } \ 337 while (0) 338 339 static void 340 check_random (mpfr_prec_t p) 341 { 342 mpfr_t a1, a2, b, bs, c, cs; 343 int r; 344 int i, inexact1, inexact2; 345 346 mpfr_inits2 (p, a1, a2, b, c, (mpfr_ptr) 0); 347 348 for (i = 0 ; i < 500 ; i++) 349 { 350 mpfr_urandom (b, RANDS, MPFR_RNDA); 351 mpfr_urandom (c, RANDS, MPFR_RNDA); 352 if (MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c)) 353 { 354 if (randlimb () & 1) 355 mpfr_neg (b, b, MPFR_RNDN); 356 if (randlimb () & 1) 357 mpfr_neg (c, c, MPFR_RNDN); 358 if (MPFR_GET_EXP(b) < MPFR_GET_EXP(c)) 359 { 360 /* Exchange b and c, except the signs (actually, the sign 361 of cs doesn't matter). */ 362 MPFR_ALIAS (bs, c, MPFR_SIGN (b), MPFR_EXP (c)); 363 MPFR_ALIAS (cs, b, MPFR_SIGN (c), MPFR_EXP (b)); 364 } 365 else 366 { 367 MPFR_ALIAS (bs, b, MPFR_SIGN (b), MPFR_EXP (b)); 368 MPFR_ALIAS (cs, c, MPFR_SIGN (c), MPFR_EXP (c)); 369 } 370 for (r = 0 ; r < MPFR_RND_MAX ; r++) 371 { 372 mpfr_flags_t flags1, flags2; 373 374 if (r == MPFR_RNDF) /* inexact makes no sense, moreover 375 mpfr_add1 and mpfr_add1sp could 376 return different values */ 377 continue; 378 379 mpfr_clear_flags (); 380 inexact1 = mpfr_add1 (a1, bs, cs, (mpfr_rnd_t) r); 381 flags1 = __gmpfr_flags; 382 mpfr_clear_flags (); 383 inexact2 = mpfr_add1sp (a2, b, c, (mpfr_rnd_t) r); 384 flags2 = __gmpfr_flags; 385 if (! mpfr_equal_p (a1, a2)) 386 STD_ERROR; 387 if (inexact1 != inexact2) 388 STD_ERROR2; 389 MPFR_ASSERTN (flags1 == flags2); 390 } 391 } 392 } 393 394 mpfr_clears (a1, a2, b, c, (mpfr_ptr) 0); 395 } 396 397 static void 398 check_special (void) 399 { 400 mpfr_t a1, a2, b, c; 401 int r; 402 mpfr_prec_t p; 403 int i = -1, inexact1, inexact2; 404 405 mpfr_inits (a1, a2, b, c, (mpfr_ptr) 0); 406 407 for (r = 0 ; r < MPFR_RND_MAX ; r++) 408 { 409 if (r == MPFR_RNDF) 410 continue; /* inexact makes no sense, mpfr_add1 and mpfr_add1sp 411 could differ */ 412 SET_PREC(53); 413 mpfr_set_str1 (b, "1@100"); 414 mpfr_set_str1 (c, "1@1"); 415 inexact1 = mpfr_add1(a1, b, c, (mpfr_rnd_t) r); 416 inexact2 = mpfr_add1sp(a2, b, c, (mpfr_rnd_t) r); 417 if (mpfr_cmp(a1, a2)) 418 STD_ERROR; 419 if (inexact1 != inexact2) 420 STD_ERROR2; 421 mpfr_set_str_binary (b, "1E53"); 422 mpfr_set_str_binary (c, "1E0"); 423 inexact1 = mpfr_add1(a1, b, c, (mpfr_rnd_t) r); 424 inexact2 = mpfr_add1sp(a2, b, c, (mpfr_rnd_t) r); 425 if (mpfr_cmp(a1, a2)) 426 STD_ERROR; 427 if (inexact1 != inexact2) 428 STD_ERROR2; 429 } 430 431 mpfr_set_prec (c, 2); 432 mpfr_set_prec (a1, 2); 433 mpfr_set_prec (a2, 2); 434 mpfr_set_str_binary (c, "1.0e1"); 435 mpfr_set_str_binary (a2, "1.1e-1"); 436 mpfr_set_str_binary (a1, "0.11E2"); 437 mpfr_add1sp (a2, c, a2, MPFR_RNDN); 438 if (mpfr_cmp (a1, a2)) 439 { 440 printf ("Regression reuse test failed!\n"); 441 exit (1); 442 } 443 444 mpfr_set_prec (a1, 63); 445 mpfr_set_prec (b, 63); 446 mpfr_set_prec (c, 63); 447 mpfr_set_str_binary (b, "0.111111101010110111010100110101010110000101111011011011100011001E-3"); 448 mpfr_set_str_binary (c, "0.101111111101110000001100001000011000011011010001010011111100111E-4"); 449 mpfr_clear_inexflag (); 450 mpfr_add1sp (a1, b, c, MPFR_RNDN); 451 MPFR_ASSERTN (mpfr_inexflag_p ()); 452 453 mpfr_clears (a1, a2, b, c, (mpfr_ptr) 0); 454 } 455