tadd.c revision 1.1.1.3.4.1
1 /* Test file for mpfr_add and mpfr_sub. 2 3 Copyright 1999-2018 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 http://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 #define N 30000 24 25 #include <float.h> 26 27 #include "mpfr-test.h" 28 29 /* If the precisions are the same, we want to test both mpfr_add1sp 30 and mpfr_add1. */ 31 32 /* FIXME: modify check() to test the ternary value and the flags. */ 33 34 static int usesp; 35 36 static int 37 test_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) 38 { 39 int res; 40 #ifdef CHECK_EXTERNAL 41 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c); 42 if (ok) 43 { 44 mpfr_print_raw (b); 45 printf (" "); 46 mpfr_print_raw (c); 47 } 48 #endif 49 if (usesp || MPFR_ARE_SINGULAR(b,c) || MPFR_SIGN(b) != MPFR_SIGN(c)) 50 res = mpfr_add (a, b, c, rnd_mode); 51 else 52 { 53 if (MPFR_GET_EXP(b) < MPFR_GET_EXP(c)) 54 res = mpfr_add1(a, c, b, rnd_mode); 55 else 56 res = mpfr_add1(a, b, c, rnd_mode); 57 } 58 #ifdef CHECK_EXTERNAL 59 if (ok) 60 { 61 printf (" "); 62 mpfr_print_raw (a); 63 printf ("\n"); 64 } 65 #endif 66 return res; 67 } 68 69 /* checks that xs+ys gives the expected result zs */ 70 static void 71 check (const char *xs, const char *ys, mpfr_rnd_t rnd_mode, 72 unsigned int px, unsigned int py, unsigned int pz, const char *zs) 73 { 74 mpfr_t xx,yy,zz; 75 76 mpfr_init2 (xx, px); 77 mpfr_init2 (yy, py); 78 mpfr_init2 (zz, pz); 79 80 mpfr_set_str1 (xx, xs); 81 mpfr_set_str1 (yy, ys); 82 test_add (zz, xx, yy, rnd_mode); 83 if (mpfr_cmp_str1 (zz, zs) ) 84 { 85 printf ("expected sum is %s, got ", zs); 86 mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN); 87 printf ("\nmpfr_add failed for x=%s y=%s with rnd_mode=%s\n", 88 xs, ys, mpfr_print_rnd_mode (rnd_mode)); 89 exit (1); 90 } 91 mpfr_clears (xx, yy, zz, (mpfr_ptr) 0); 92 } 93 94 static void 95 check2b (const char *xs, int px, 96 const char *ys, int py, 97 const char *rs, int pz, 98 mpfr_rnd_t rnd_mode) 99 { 100 mpfr_t xx, yy, zz; 101 102 mpfr_init2 (xx,px); 103 mpfr_init2 (yy,py); 104 mpfr_init2 (zz,pz); 105 mpfr_set_str_binary (xx, xs); 106 mpfr_set_str_binary (yy, ys); 107 test_add (zz, xx, yy, rnd_mode); 108 if (mpfr_cmp_str (zz, rs, 2, MPFR_RNDN)) 109 { 110 printf ("(2) x=%s,%d y=%s,%d pz=%d,rnd=%s\n", 111 xs, px, ys, py, pz, mpfr_print_rnd_mode (rnd_mode)); 112 printf ("got "); mpfr_dump (zz); 113 mpfr_set_str(zz, rs, 2, MPFR_RNDN); 114 printf ("instead of "); mpfr_dump (zz); 115 exit (1); 116 } 117 mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz); 118 } 119 120 static void 121 check64 (void) 122 { 123 mpfr_t x, t, u; 124 125 mpfr_init (x); 126 mpfr_init (t); 127 mpfr_init (u); 128 129 mpfr_set_prec (x, 29); 130 mpfr_set_str_binary (x, "1.1101001000101111011010010110e-3"); 131 mpfr_set_prec (t, 58); 132 mpfr_set_str_binary (t, "0.11100010011111001001100110010111110110011000000100101E-1"); 133 mpfr_set_prec (u, 29); 134 test_add (u, x, t, MPFR_RNDD); 135 mpfr_set_str_binary (t, "1.0101011100001000011100111110e-1"); 136 if (mpfr_cmp (u, t)) 137 { 138 printf ("mpfr_add(u, x, t) failed for prec(x)=29, prec(t)=58\n"); 139 printf ("expected "); mpfr_out_str (stdout, 2, 29, t, MPFR_RNDN); 140 puts (""); 141 printf ("got "); mpfr_out_str (stdout, 2, 29, u, MPFR_RNDN); 142 puts (""); 143 exit(1); 144 } 145 146 mpfr_set_prec (x, 4); 147 mpfr_set_str_binary (x, "-1.0E-2"); 148 mpfr_set_prec (t, 2); 149 mpfr_set_str_binary (t, "-1.1e-2"); 150 mpfr_set_prec (u, 2); 151 test_add (u, x, t, MPFR_RNDN); 152 if (MPFR_MANT(u)[0] << 2) 153 { 154 printf ("result not normalized for prec=2\n"); 155 mpfr_dump (u); 156 exit (1); 157 } 158 mpfr_set_str_binary (t, "-1.0e-1"); 159 if (mpfr_cmp (u, t)) 160 { 161 printf ("mpfr_add(u, x, t) failed for prec(x)=4, prec(t)=2\n"); 162 printf ("expected -1.0e-1\n"); 163 printf ("got "); mpfr_out_str (stdout, 2, 4, u, MPFR_RNDN); 164 puts (""); 165 exit (1); 166 } 167 168 mpfr_set_prec (x, 8); 169 mpfr_set_str_binary (x, "-0.10011010"); /* -77/128 */ 170 mpfr_set_prec (t, 4); 171 mpfr_set_str_binary (t, "-1.110e-5"); /* -7/128 */ 172 mpfr_set_prec (u, 4); 173 test_add (u, x, t, MPFR_RNDN); /* should give -5/8 */ 174 mpfr_set_str_binary (t, "-1.010e-1"); 175 if (mpfr_cmp (u, t)) { 176 printf ("mpfr_add(u, x, t) failed for prec(x)=8, prec(t)=4\n"); 177 printf ("expected -1.010e-1\n"); 178 printf ("got "); mpfr_out_str (stdout, 2, 4, u, MPFR_RNDN); 179 puts (""); 180 exit (1); 181 } 182 183 mpfr_set_prec (x, 112); mpfr_set_prec (t, 98); mpfr_set_prec (u, 54); 184 mpfr_set_str_binary (x, "-0.11111100100000000011000011100000101101010001000111E-401"); 185 mpfr_set_str_binary (t, "0.10110000100100000101101100011111111011101000111000101E-464"); 186 test_add (u, x, t, MPFR_RNDN); 187 if (mpfr_cmp (u, x)) 188 { 189 printf ("mpfr_add(u, x, t) failed for prec(x)=112, prec(t)=98\n"); 190 exit (1); 191 } 192 193 mpfr_set_prec (x, 92); mpfr_set_prec (t, 86); mpfr_set_prec (u, 53); 194 mpfr_set_str (x, "-5.03525136761487735093e-74", 10, MPFR_RNDN); 195 mpfr_set_str (t, "8.51539046314262304109e-91", 10, MPFR_RNDN); 196 test_add (u, x, t, MPFR_RNDN); 197 if (mpfr_cmp_str1 (u, "-5.0352513676148773509283672e-74") ) 198 { 199 printf ("mpfr_add(u, x, t) failed for prec(x)=92, prec(t)=86\n"); 200 exit (1); 201 } 202 203 mpfr_set_prec(x, 53); mpfr_set_prec(t, 76); mpfr_set_prec(u, 76); 204 mpfr_set_str_binary(x, "-0.10010010001001011011110000000000001010011011011110001E-32"); 205 mpfr_set_str_binary(t, "-0.1011000101110010000101111111011111010001110011110111100110101011110010011111"); 206 mpfr_sub(u, x, t, MPFR_RNDU); 207 mpfr_set_str_binary(t, "0.1011000101110010000101111111011100111111101010011011110110101011101000000100"); 208 if (mpfr_cmp(u,t)) 209 { 210 printf ("expect "); mpfr_dump (t); 211 printf ("mpfr_add failed for precisions 53-76\n"); 212 exit (1); 213 } 214 mpfr_set_prec(x, 53); mpfr_set_prec(t, 108); mpfr_set_prec(u, 108); 215 mpfr_set_str_binary(x, "-0.10010010001001011011110000000000001010011011011110001E-32"); 216 mpfr_set_str_binary(t, "-0.101100010111001000010111111101111101000111001111011110011010101111001001111000111011001110011000000000111111"); 217 mpfr_sub(u, x, t, MPFR_RNDU); 218 mpfr_set_str_binary(t, "0.101100010111001000010111111101110011111110101001101111011010101110100000001011000010101110011000000000111111"); 219 if (mpfr_cmp(u,t)) 220 { 221 printf ("expect "); mpfr_dump (t); 222 printf ("mpfr_add failed for precisions 53-108\n"); 223 exit (1); 224 } 225 mpfr_set_prec(x, 97); mpfr_set_prec(t, 97); mpfr_set_prec(u, 97); 226 mpfr_set_str_binary(x, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010000000000000000E-39"); 227 mpfr_set_ui(t, 1, MPFR_RNDN); 228 test_add (u, x, t, MPFR_RNDN); 229 mpfr_set_str_binary(x, "0.1000000000000000000000000000000000000000111110110000100000000101100011011110100000101111100010001E1"); 230 if (mpfr_cmp(u,x)) 231 { 232 printf ("mpfr_add failed for precision 97\n"); 233 exit (1); 234 } 235 mpfr_set_prec(x, 128); mpfr_set_prec(t, 128); mpfr_set_prec(u, 128); 236 mpfr_set_str_binary(x, "0.10101011111001001010111011001000101100111101000000111111111011010100001100011101010001010111111101111010100110111111100101100010E-4"); 237 mpfr_set(t, x, MPFR_RNDN); 238 mpfr_sub(u, x, t, MPFR_RNDN); 239 mpfr_set_prec(x, 96); mpfr_set_prec(t, 96); mpfr_set_prec(u, 96); 240 mpfr_set_str_binary(x, "0.111000000001110100111100110101101001001010010011010011100111100011010100011001010011011011000010E-4"); 241 mpfr_set(t, x, MPFR_RNDN); 242 mpfr_sub(u, x, t, MPFR_RNDN); 243 mpfr_set_prec(x, 85); mpfr_set_prec(t, 85); mpfr_set_prec(u, 85); 244 mpfr_set_str_binary(x, "0.1111101110100110110110100010101011101001100010100011110110110010010011101100101111100E-4"); 245 mpfr_set_str_binary(t, "0.1111101110100110110110100010101001001000011000111000011101100101110100001110101010110E-4"); 246 mpfr_sub(u, x, t, MPFR_RNDU); 247 mpfr_sub(x, x, t, MPFR_RNDU); 248 if (mpfr_cmp(x, u) != 0) 249 { 250 printf ("Error in mpfr_sub: u=x-t and x=x-t give different results\n"); 251 exit (1); 252 } 253 if (! MPFR_IS_NORMALIZED (u)) 254 { 255 printf ("Error in mpfr_sub: result is not msb-normalized (1)\n"); 256 exit (1); 257 } 258 mpfr_set_prec(x, 65); mpfr_set_prec(t, 65); mpfr_set_prec(u, 65); 259 mpfr_set_str_binary(x, "0.10011010101000110101010000000011001001001110001011101011111011101E623"); 260 mpfr_set_str_binary(t, "0.10011010101000110101010000000011001001001110001011101011111011100E623"); 261 mpfr_sub(u, x, t, MPFR_RNDU); 262 if (mpfr_cmp_ui_2exp(u, 1, 558)) 263 { /* 2^558 */ 264 printf ("Error (1) in mpfr_sub\n"); 265 exit (1); 266 } 267 268 mpfr_set_prec(x, 64); mpfr_set_prec(t, 64); mpfr_set_prec(u, 64); 269 mpfr_set_str_binary(x, "0.1000011110101111011110111111000011101011101111101101101100000100E-220"); 270 mpfr_set_str_binary(t, "0.1000011110101111011110111111000011101011101111101101010011111101E-220"); 271 test_add (u, x, t, MPFR_RNDU); 272 if ((MPFR_MANT(u)[0] & 1) != 1) 273 { 274 printf ("error in mpfr_add with rnd_mode=MPFR_RNDU\n"); 275 printf ("b= "); mpfr_dump (x); 276 printf ("c= "); mpfr_dump (t); 277 printf ("b+c="); mpfr_dump (u); 278 exit (1); 279 } 280 281 /* bug found by Norbert Mueller, 14 Sep 2000 */ 282 mpfr_set_prec(x, 56); mpfr_set_prec(t, 83); mpfr_set_prec(u, 10); 283 mpfr_set_str_binary(x, "0.10001001011011001111101100110100000101111010010111010111E-7"); 284 mpfr_set_str_binary(t, "0.10001001011011001111101100110100000101111010010111010111000000000111110110110000100E-7"); 285 mpfr_sub(u, x, t, MPFR_RNDU); 286 287 /* array bound write found by Norbert Mueller, 26 Sep 2000 */ 288 mpfr_set_prec(x, 109); mpfr_set_prec(t, 153); mpfr_set_prec(u, 95); 289 mpfr_set_str_binary(x,"0.1001010000101011101100111000110001111111111111111111111111111111111111111111111111111111111111100000000000000E33"); 290 mpfr_set_str_binary(t,"-0.100101000010101110110011100011000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011100101101000000100100001100110111E33"); 291 test_add (u, x, t, MPFR_RNDN); 292 293 /* array bound writes found by Norbert Mueller, 27 Sep 2000 */ 294 mpfr_set_prec(x, 106); mpfr_set_prec(t, 53); mpfr_set_prec(u, 23); 295 mpfr_set_str_binary(x, "-0.1000011110101111111001010001000100001011000000000000000000000000000000000000000000000000000000000000000000E-59"); 296 mpfr_set_str_binary(t, "-0.10000111101011111110010100010001101100011100110100000E-59"); 297 mpfr_sub(u, x, t, MPFR_RNDN); 298 mpfr_set_prec(x, 177); mpfr_set_prec(t, 217); mpfr_set_prec(u, 160); 299 mpfr_set_str_binary(x, "-0.111010001011010000111001001010010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E35"); 300 mpfr_set_str_binary(t, "0.1110100010110100001110010010100100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111011010011100001111001E35"); 301 test_add (u, x, t, MPFR_RNDN); 302 mpfr_set_prec(x, 214); mpfr_set_prec(t, 278); mpfr_set_prec(u, 207); 303 mpfr_set_str_binary(x, "0.1000100110100110101101101101000000010000100111000001001110001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E66"); 304 mpfr_set_str_binary(t, "-0.10001001101001101011011011010000000100001001110000010011100010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001111011111001001100011E66"); 305 test_add (u, x, t, MPFR_RNDN); 306 mpfr_set_prec(x, 32); mpfr_set_prec(t, 247); mpfr_set_prec(u, 223); 307 mpfr_set_str_binary(x, "0.10000000000000000000000000000000E1"); 308 mpfr_set_str_binary(t, "0.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100000110001110100000100011110000101110110011101110100110110111111011010111100100000000000000000000000000E0"); 309 mpfr_sub(u, x, t, MPFR_RNDN); 310 if (! MPFR_IS_NORMALIZED (u)) 311 { 312 printf ("Error in mpfr_sub: result is not msb-normalized (2)\n"); 313 exit (1); 314 } 315 316 /* bug found by Nathalie Revol, 21 March 2001 */ 317 mpfr_set_prec (x, 65); 318 mpfr_set_prec (t, 65); 319 mpfr_set_prec (u, 65); 320 mpfr_set_str_binary (x, "0.11100100101101001100111011111111110001101001000011101001001010010E-35"); 321 mpfr_set_str_binary (t, "0.10000000000000000000000000000000000001110010010110100110011110000E1"); 322 mpfr_sub (u, t, x, MPFR_RNDU); 323 if (! MPFR_IS_NORMALIZED (u)) 324 { 325 printf ("Error in mpfr_sub: result is not msb-normalized (3)\n"); 326 exit (1); 327 } 328 329 /* bug found by Fabrice Rouillier, 27 Mar 2001 */ 330 mpfr_set_prec (x, 107); 331 mpfr_set_prec (t, 107); 332 mpfr_set_prec (u, 107); 333 mpfr_set_str_binary (x, "0.10111001001111010010001000000010111111011011011101000001001000101000000000000000000000000000000000000000000E315"); 334 mpfr_set_str_binary (t, "0.10000000000000000000000000000000000101110100100101110110000001100101011111001000011101111100100100111011000E350"); 335 mpfr_sub (u, x, t, MPFR_RNDU); 336 if (! MPFR_IS_NORMALIZED (u)) 337 { 338 printf ("Error in mpfr_sub: result is not msb-normalized (4)\n"); 339 exit (1); 340 } 341 342 /* checks that NaN flag is correctly reset */ 343 mpfr_set_ui (t, 1, MPFR_RNDN); 344 mpfr_set_ui (u, 1, MPFR_RNDN); 345 mpfr_set_nan (x); 346 test_add (x, t, u, MPFR_RNDN); 347 if (mpfr_cmp_ui (x, 2)) 348 { 349 printf ("Error in mpfr_add: 1+1 gives "); 350 mpfr_out_str(stdout, 10, 0, x, MPFR_RNDN); 351 exit (1); 352 } 353 354 mpfr_clear(x); mpfr_clear(t); mpfr_clear(u); 355 } 356 357 /* check case when c does not overlap with a, but both b and c count 358 for rounding */ 359 static void 360 check_case_1b (void) 361 { 362 mpfr_t a, b, c; 363 unsigned int prec_a, prec_b, prec_c, dif; 364 365 mpfr_init (a); 366 mpfr_init (b); 367 mpfr_init (c); 368 369 { 370 prec_a = MPFR_PREC_MIN + (randlimb () % 63); 371 mpfr_set_prec (a, prec_a); 372 for (prec_b = prec_a + 2; prec_b <= 64; prec_b++) 373 { 374 dif = prec_b - prec_a; 375 mpfr_set_prec (b, prec_b); 376 /* b = 1 - 2^(-prec_a) + 2^(-prec_b) */ 377 mpfr_set_ui (b, 1, MPFR_RNDN); 378 mpfr_div_2exp (b, b, dif, MPFR_RNDN); 379 mpfr_sub_ui (b, b, 1, MPFR_RNDN); 380 mpfr_div_2exp (b, b, prec_a, MPFR_RNDN); 381 mpfr_add_ui (b, b, 1, MPFR_RNDN); 382 for (prec_c = dif; prec_c <= 64; prec_c++) 383 { 384 /* c = 2^(-prec_a) - 2^(-prec_b) */ 385 mpfr_set_prec (c, prec_c); 386 mpfr_set_si (c, -1, MPFR_RNDN); 387 mpfr_div_2exp (c, c, dif, MPFR_RNDN); 388 mpfr_add_ui (c, c, 1, MPFR_RNDN); 389 mpfr_div_2exp (c, c, prec_a, MPFR_RNDN); 390 test_add (a, b, c, MPFR_RNDN); 391 if (mpfr_cmp_ui (a, 1) != 0) 392 { 393 printf ("case (1b) failed for prec_a=%u, prec_b=%u," 394 " prec_c=%u\n", prec_a, prec_b, prec_c); 395 printf ("b="); mpfr_dump (b); 396 printf ("c="); mpfr_dump (c); 397 printf ("a="); mpfr_dump (a); 398 exit (1); 399 } 400 } 401 } 402 } 403 404 mpfr_clear (a); 405 mpfr_clear (b); 406 mpfr_clear (c); 407 } 408 409 /* check case when c overlaps with a */ 410 static void 411 check_case_2 (void) 412 { 413 mpfr_t a, b, c, d; 414 415 mpfr_init2 (a, 300); 416 mpfr_init2 (b, 800); 417 mpfr_init2 (c, 500); 418 mpfr_init2 (d, 800); 419 420 mpfr_set_str_binary(a, "1E110"); /* a = 2^110 */ 421 mpfr_set_str_binary(b, "1E900"); /* b = 2^900 */ 422 mpfr_set_str_binary(c, "1E500"); /* c = 2^500 */ 423 test_add (c, c, a, MPFR_RNDZ); /* c = 2^500 + 2^110 */ 424 mpfr_sub (d, b, c, MPFR_RNDZ); /* d = 2^900 - 2^500 - 2^110 */ 425 test_add (b, b, c, MPFR_RNDZ); /* b = 2^900 + 2^500 + 2^110 */ 426 test_add (a, b, d, MPFR_RNDZ); /* a = 2^901 */ 427 if (mpfr_cmp_ui_2exp (a, 1, 901)) 428 { 429 printf ("b + d fails for b=2^900+2^500+2^110, d=2^900-2^500-2^110\n"); 430 printf ("expected 1.0e901, got "); 431 mpfr_out_str (stdout, 2, 0, a, MPFR_RNDN); 432 printf ("\n"); 433 exit (1); 434 } 435 436 mpfr_clear (a); 437 mpfr_clear (b); 438 mpfr_clear (c); 439 mpfr_clear (d); 440 } 441 442 /* checks when source and destination are equal */ 443 static void 444 check_same (void) 445 { 446 mpfr_t x; 447 448 mpfr_init(x); mpfr_set_ui(x, 1, MPFR_RNDZ); 449 test_add (x, x, x, MPFR_RNDZ); 450 if (mpfr_cmp_ui (x, 2)) 451 { 452 printf ("Error when all 3 operands are equal\n"); 453 exit (1); 454 } 455 mpfr_clear(x); 456 } 457 458 #define check53(x, y, r, z) check(x, y, r, 53, 53, 53, z) 459 460 #define MAX_PREC 256 461 462 static void 463 check_inexact (void) 464 { 465 mpfr_t x, y, z, u; 466 mpfr_prec_t px, py, pu, pz; 467 int inexact, cmp; 468 mpfr_rnd_t rnd; 469 470 mpfr_init (x); 471 mpfr_init (y); 472 mpfr_init (z); 473 mpfr_init (u); 474 475 mpfr_set_prec (x, 2); 476 mpfr_set_str_binary (x, "0.1E-4"); 477 mpfr_set_prec (u, 33); 478 mpfr_set_str_binary (u, "0.101110100101101100000000111100000E-1"); 479 mpfr_set_prec (y, 31); 480 inexact = test_add (y, x, u, MPFR_RNDN); 481 482 if (inexact != 0) 483 { 484 printf ("Wrong ternary value (2): expected 0, got %d\n", inexact); 485 exit (1); 486 } 487 488 mpfr_set_prec (x, 2); 489 mpfr_set_str_binary (x, "0.1E-4"); 490 mpfr_set_prec (u, 33); 491 mpfr_set_str_binary (u, "0.101110100101101100000000111100000E-1"); 492 mpfr_set_prec (y, 28); 493 inexact = test_add (y, x, u, MPFR_RNDN); 494 495 if (inexact != 0) 496 { 497 printf ("Wrong ternary value (1): expected 0, got %d\n", inexact); 498 exit (1); 499 } 500 501 for (px = 2; px < MAX_PREC; px++) 502 { 503 mpfr_set_prec (x, px); 504 505 do 506 { 507 mpfr_urandomb (x, RANDS); 508 } 509 while (mpfr_cmp_ui (x, 0) == 0); 510 511 for (pu = 2; pu < MAX_PREC; pu++) 512 { 513 mpfr_set_prec (u, pu); 514 515 do 516 { 517 mpfr_urandomb (u, RANDS); 518 } 519 while (mpfr_cmp_ui (u, 0) == 0); 520 521 py = MPFR_PREC_MIN + (randlimb () % (MAX_PREC - 1)); 522 mpfr_set_prec (y, py); 523 pz = mpfr_cmpabs (x, u) >= 0 ? 524 MPFR_EXP(x) - MPFR_EXP(u) : 525 MPFR_EXP(u) - MPFR_EXP(x); 526 /* x + u is exactly representable with precision 527 abs(EXP(x)-EXP(u)) + max(prec(x), prec(u)) + 1 */ 528 pz = pz + MAX(MPFR_PREC(x), MPFR_PREC(u)) + 1; 529 mpfr_set_prec (z, pz); 530 531 rnd = RND_RAND_NO_RNDF (); 532 inexact = test_add (z, x, u, rnd); 533 if (inexact != 0) 534 { 535 printf ("z <- x + u should be exact\n"); 536 printf ("x="); mpfr_dump (x); 537 printf ("u="); mpfr_dump (u); 538 printf ("z="); mpfr_dump (z); 539 exit (1); 540 } 541 542 rnd = RND_RAND_NO_RNDF (); 543 inexact = test_add (y, x, u, rnd); 544 cmp = mpfr_cmp (y, z); 545 if ((inexact == 0 && cmp != 0) || 546 (inexact > 0 && cmp <= 0) || 547 (inexact < 0 && cmp >= 0)) 548 { 549 printf ("Wrong ternary value for rnd=%s\n", 550 mpfr_print_rnd_mode (rnd)); 551 printf ("expected %d, got %d\n", cmp, inexact); 552 printf ("x="); mpfr_dump (x); 553 printf ("u="); mpfr_dump (u); 554 printf ("y= "); mpfr_dump (y); 555 printf ("x+u="); mpfr_dump (z); 556 exit (1); 557 } 558 } 559 } 560 561 mpfr_clear (x); 562 mpfr_clear (y); 563 mpfr_clear (z); 564 mpfr_clear (u); 565 } 566 567 static void 568 check_nans (void) 569 { 570 mpfr_t s, x, y; 571 572 mpfr_init2 (x, 8L); 573 mpfr_init2 (y, 8L); 574 mpfr_init2 (s, 8L); 575 576 /* +inf + -inf == nan */ 577 mpfr_set_inf (x, 1); 578 mpfr_set_inf (y, -1); 579 test_add (s, x, y, MPFR_RNDN); 580 MPFR_ASSERTN (mpfr_nan_p (s)); 581 582 /* +inf + 1 == +inf */ 583 mpfr_set_inf (x, 1); 584 mpfr_set_ui (y, 1L, MPFR_RNDN); 585 test_add (s, x, y, MPFR_RNDN); 586 MPFR_ASSERTN (mpfr_inf_p (s)); 587 MPFR_ASSERTN (mpfr_sgn (s) > 0); 588 589 /* -inf + 1 == -inf */ 590 mpfr_set_inf (x, -1); 591 mpfr_set_ui (y, 1L, MPFR_RNDN); 592 test_add (s, x, y, MPFR_RNDN); 593 MPFR_ASSERTN (mpfr_inf_p (s)); 594 MPFR_ASSERTN (mpfr_sgn (s) < 0); 595 596 /* 1 + +inf == +inf */ 597 mpfr_set_ui (x, 1L, MPFR_RNDN); 598 mpfr_set_inf (y, 1); 599 test_add (s, x, y, MPFR_RNDN); 600 MPFR_ASSERTN (mpfr_inf_p (s)); 601 MPFR_ASSERTN (mpfr_sgn (s) > 0); 602 603 /* 1 + -inf == -inf */ 604 mpfr_set_ui (x, 1L, MPFR_RNDN); 605 mpfr_set_inf (y, -1); 606 test_add (s, x, y, MPFR_RNDN); 607 MPFR_ASSERTN (mpfr_inf_p (s)); 608 MPFR_ASSERTN (mpfr_sgn (s) < 0); 609 610 mpfr_clear (x); 611 mpfr_clear (y); 612 mpfr_clear (s); 613 } 614 615 static void 616 check_alloc (void) 617 { 618 mpfr_t a; 619 620 mpfr_init2 (a, 10000); 621 mpfr_set_prec (a, 53); 622 mpfr_set_ui (a, 15236, MPFR_RNDN); 623 test_add (a, a, a, MPFR_RNDN); 624 mpfr_mul (a, a, a, MPFR_RNDN); 625 mpfr_div (a, a, a, MPFR_RNDN); 626 mpfr_sub (a, a, a, MPFR_RNDN); 627 mpfr_clear (a); 628 } 629 630 static void 631 check_overflow (void) 632 { 633 mpfr_t a, b, c; 634 mpfr_prec_t prec_a, prec_b, prec_c; 635 int r, up; 636 637 mpfr_init (a); 638 mpfr_init (b); 639 mpfr_init (c); 640 641 RND_LOOP_NO_RNDF (r) 642 for (prec_a = 2; prec_a <= 128; prec_a += 2) 643 for (prec_b = 2; prec_b <= 128; prec_b += 2) 644 for (prec_c = 2; prec_c <= 128; prec_c += 2) 645 { 646 mpfr_set_prec (a, prec_a); 647 mpfr_set_prec (b, prec_b); 648 mpfr_set_prec (c, prec_c); 649 650 mpfr_setmax (b, mpfr_get_emax ()); 651 652 up = r == MPFR_RNDA || r == MPFR_RNDU || r == MPFR_RNDN; 653 654 /* set c with overlap with bits of b: will always overflow */ 655 mpfr_set_ui_2exp (c, 1, mpfr_get_emax () - prec_b / 2, MPFR_RNDN); 656 mpfr_nextbelow (c); 657 mpfr_clear_overflow (); 658 test_add (a, b, c, (mpfr_rnd_t) r); 659 if (!mpfr_overflow_p () || (up && !mpfr_inf_p (a))) 660 { 661 printf ("No overflow (1) in check_overflow for rnd=%s\n", 662 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 663 printf ("b="); mpfr_dump (b); 664 printf ("c="); mpfr_dump (c); 665 printf ("a="); mpfr_dump (a); 666 exit (1); 667 } 668 669 if (r == MPFR_RNDZ || r == MPFR_RNDD || prec_a >= prec_b + prec_c) 670 continue; 671 672 /* set c to 111...111 so that ufp(c) = 1/2 ulp(b): will only overflow 673 when prec_a < prec_b + prec_c, and rounding up or to nearest */ 674 mpfr_set_ui_2exp (c, 1, mpfr_get_emax () - prec_b, MPFR_RNDN); 675 mpfr_nextbelow (c); 676 mpfr_clear_overflow (); 677 test_add (a, b, c, (mpfr_rnd_t) r); 678 /* b + c is exactly representable iff prec_a >= prec_b + prec_c */ 679 if (!mpfr_overflow_p () || !mpfr_inf_p (a)) 680 { 681 printf ("No overflow (2) in check_overflow for rnd=%s\n", 682 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 683 printf ("b="); mpfr_dump (b); 684 printf ("c="); mpfr_dump (c); 685 printf ("a="); mpfr_dump (a); 686 exit (1); 687 } 688 } 689 690 mpfr_set_prec (b, 256); 691 mpfr_setmax (b, mpfr_get_emax ()); 692 mpfr_set_prec (c, 256); 693 mpfr_set_ui (c, 1, MPFR_RNDN); 694 mpfr_set_exp (c, mpfr_get_emax () - 512); 695 mpfr_set_prec (a, 256); 696 mpfr_clear_overflow (); 697 mpfr_add (a, b, c, MPFR_RNDU); 698 if (!mpfr_overflow_p ()) 699 { 700 printf ("No overflow (3) in check_overflow\n"); 701 printf ("b="); mpfr_dump (b); 702 printf ("c="); mpfr_dump (c); 703 printf ("a="); mpfr_dump (a); 704 exit (1); 705 } 706 707 mpfr_clear (a); 708 mpfr_clear (b); 709 mpfr_clear (c); 710 } 711 712 static void 713 check_1111 (void) 714 { 715 mpfr_t one; 716 long n; 717 718 mpfr_init2 (one, MPFR_PREC_MIN); 719 mpfr_set_ui (one, 1, MPFR_RNDN); 720 for (n = 0; n < N; n++) 721 { 722 mpfr_prec_t prec_a, prec_b, prec_c; 723 mpfr_exp_t tb=0, tc, diff; 724 mpfr_t a, b, c, s; 725 int m = 512; 726 int sb, sc; 727 int inex_a, inex_s; 728 mpfr_rnd_t rnd_mode; 729 730 prec_a = MPFR_PREC_MIN + (randlimb () % m); 731 prec_b = MPFR_PREC_MIN + (randlimb () % m); 732 /* we need prec_c > 1 so that % (prec_c - 1) is well defined below */ 733 do prec_c = MPFR_PREC_MIN + (randlimb () % m); while (prec_c == 1); 734 mpfr_init2 (a, prec_a); 735 mpfr_init2 (b, prec_b); 736 mpfr_init2 (c, prec_c); 737 /* we need prec_b - (sb != 2) > 0 below */ 738 do sb = randlimb () % 3; while (prec_b - (sb != 2) == 0); 739 if (sb != 0) 740 { 741 tb = 1 + (randlimb () % (prec_b - (sb != 2))); 742 mpfr_div_2ui (b, one, tb, MPFR_RNDN); 743 if (sb == 2) 744 mpfr_neg (b, b, MPFR_RNDN); 745 test_add (b, b, one, MPFR_RNDN); 746 } 747 else 748 mpfr_set (b, one, MPFR_RNDN); 749 tc = 1 + (randlimb () % (prec_c - 1)); 750 mpfr_div_2ui (c, one, tc, MPFR_RNDN); 751 sc = randlimb () % 2; 752 if (sc) 753 mpfr_neg (c, c, MPFR_RNDN); 754 test_add (c, c, one, MPFR_RNDN); 755 diff = (randlimb () % (2*m)) - m; 756 mpfr_mul_2si (c, c, diff, MPFR_RNDN); 757 rnd_mode = RND_RAND_NO_RNDF (); 758 inex_a = test_add (a, b, c, rnd_mode); 759 mpfr_init2 (s, MPFR_PREC_MIN + 2*m); 760 inex_s = test_add (s, b, c, MPFR_RNDN); /* exact */ 761 if (inex_s) 762 { 763 printf ("check_1111: result should have been exact.\n"); 764 exit (1); 765 } 766 inex_s = mpfr_prec_round (s, prec_a, rnd_mode); 767 if ((inex_a < 0 && inex_s >= 0) || 768 (inex_a == 0 && inex_s != 0) || 769 (inex_a > 0 && inex_s <= 0) || 770 !mpfr_equal_p (a, s)) 771 { 772 printf ("check_1111: results are different.\n"); 773 printf ("prec_a = %d, prec_b = %d, prec_c = %d\n", 774 (int) prec_a, (int) prec_b, (int) prec_c); 775 printf ("tb = %d, tc = %d, diff = %d, rnd = %s\n", 776 (int) tb, (int) tc, (int) diff, 777 mpfr_print_rnd_mode (rnd_mode)); 778 printf ("sb = %d, sc = %d\n", sb, sc); 779 printf ("a = "); mpfr_dump (a); 780 printf ("s = "); mpfr_dump (s); 781 printf ("inex_a = %d, inex_s = %d\n", inex_a, inex_s); 782 exit (1); 783 } 784 mpfr_clear (a); 785 mpfr_clear (b); 786 mpfr_clear (c); 787 mpfr_clear (s); 788 } 789 mpfr_clear (one); 790 } 791 792 static void 793 check_1minuseps (void) 794 { 795 static mpfr_prec_t prec_a[] = { 796 MPFR_PREC_MIN, 30, 31, 32, 33, 62, 63, 64, 65, 126, 127, 128, 129 797 }; 798 static int supp_b[] = { 799 0, 1, 2, 3, 4, 29, 30, 31, 32, 33, 34, 35, 61, 62, 63, 64, 65, 66, 67 800 }; 801 mpfr_t a, b, c; 802 unsigned int ia, ib, ic; 803 804 mpfr_init2 (c, MPFR_PREC_MIN); 805 806 for (ia = 0; ia < numberof (prec_a); ia++) 807 for (ib = 0; ib < numberof (supp_b); ib++) 808 { 809 mpfr_prec_t prec_b; 810 int rnd_mode; 811 812 prec_b = prec_a[ia] + supp_b[ib]; 813 814 mpfr_init2 (a, prec_a[ia]); 815 mpfr_init2 (b, prec_b); 816 817 mpfr_set_ui (c, 1, MPFR_RNDN); 818 mpfr_div_ui (b, c, prec_a[ia], MPFR_RNDN); 819 mpfr_sub (b, c, b, MPFR_RNDN); /* b = 1 - 2^(-prec_a) */ 820 821 for (ic = 0; ic < numberof (supp_b); ic++) 822 for (rnd_mode = 0; rnd_mode < MPFR_RND_MAX; rnd_mode++) 823 { 824 mpfr_t s; 825 int inex_a, inex_s; 826 827 if (rnd_mode == MPFR_RNDF) 828 continue; /* inex_a, inex_s make no sense */ 829 830 mpfr_set_ui (c, 1, MPFR_RNDN); 831 mpfr_div_ui (c, c, prec_a[ia] + supp_b[ic], MPFR_RNDN); 832 inex_a = test_add (a, b, c, (mpfr_rnd_t) rnd_mode); 833 mpfr_init2 (s, 256); 834 inex_s = test_add (s, b, c, MPFR_RNDN); /* exact */ 835 if (inex_s) 836 { 837 printf ("check_1minuseps: result should have been exact " 838 "(ia = %u, ib = %u, ic = %u)\n", ia, ib, ic); 839 exit (1); 840 } 841 inex_s = mpfr_prec_round (s, prec_a[ia], (mpfr_rnd_t) rnd_mode); 842 if ((inex_a < 0 && inex_s >= 0) || 843 (inex_a == 0 && inex_s != 0) || 844 (inex_a > 0 && inex_s <= 0) || 845 !mpfr_equal_p (a, s)) 846 { 847 printf ("check_1minuseps: results are different.\n"); 848 printf ("ia = %u, ib = %u, ic = %u\n", ia, ib, ic); 849 exit (1); 850 } 851 mpfr_clear (s); 852 } 853 854 mpfr_clear (a); 855 mpfr_clear (b); 856 } 857 858 mpfr_clear (c); 859 } 860 861 /* Test case bk == 0 in add1.c (b has entirely been read and 862 c hasn't been taken into account). */ 863 static void 864 coverage_bk_eq_0 (void) 865 { 866 mpfr_t a, b, c; 867 int inex; 868 869 mpfr_init2 (a, GMP_NUMB_BITS); 870 mpfr_init2 (b, 2 * GMP_NUMB_BITS); 871 mpfr_init2 (c, GMP_NUMB_BITS); 872 873 mpfr_set_ui_2exp (b, 1, 2 * GMP_NUMB_BITS, MPFR_RNDN); 874 mpfr_sub_ui (b, b, 1, MPFR_RNDN); 875 /* b = 111...111 (in base 2) where the 1's fit 2 whole limbs */ 876 877 mpfr_set_ui_2exp (c, 1, -1, MPFR_RNDN); /* c = 1/2 */ 878 879 inex = mpfr_add (a, b, c, MPFR_RNDU); 880 mpfr_set_ui_2exp (c, 1, 2 * GMP_NUMB_BITS, MPFR_RNDN); 881 if (! mpfr_equal_p (a, c)) 882 { 883 printf ("Error in coverage_bk_eq_0\n"); 884 printf ("Expected "); 885 mpfr_dump (c); 886 printf ("Got "); 887 mpfr_dump (a); 888 exit (1); 889 } 890 MPFR_ASSERTN (inex > 0); 891 892 mpfr_clear (a); 893 mpfr_clear (b); 894 mpfr_clear (c); 895 } 896 897 static void 898 tests (void) 899 { 900 check_alloc (); 901 check_nans (); 902 check_inexact (); 903 check_case_1b (); 904 check_case_2 (); 905 check64(); 906 coverage_bk_eq_0 (); 907 908 check("293607738.0", "1.9967571564050541e-5", MPFR_RNDU, 64, 53, 53, 909 "2.9360773800002003e8"); 910 check("880524.0", "-2.0769715792901673e-5", MPFR_RNDN, 64, 53, 53, 911 "8.8052399997923023e5"); 912 check("1196426492.0", "-1.4218093058435347e-3", MPFR_RNDN, 64, 53, 53, 913 "1.1964264919985781e9"); 914 check("982013018.0", "-8.941829477291838e-7", MPFR_RNDN, 64, 53, 53, 915 "9.8201301799999905e8"); 916 check("1092583421.0", "1.0880649218158844e9", MPFR_RNDN, 64, 53, 53, 917 "2.1806483428158846e9"); 918 check("1.8476886419022969e-6", "961494401.0", MPFR_RNDN, 53, 64, 53, 919 "9.6149440100000179e8"); 920 check("-2.3222118418069868e5", "1229318102.0", MPFR_RNDN, 53, 64, 53, 921 "1.2290858808158193e9"); 922 check("-3.0399171300395734e-6", "874924868.0", MPFR_RNDN, 53, 64, 53, 923 "8.749248679999969e8"); 924 check("9.064246624706179e1", "663787413.0", MPFR_RNDN, 53, 64, 53, 925 "6.6378750364246619e8"); 926 check("-1.0954322421551264e2", "281806592.0", MPFR_RNDD, 53, 64, 53, 927 "2.8180648245677572e8"); 928 check("5.9836930386056659e-8", "1016217213.0", MPFR_RNDN, 53, 64, 53, 929 "1.0162172130000001e9"); 930 check("-1.2772161928500301e-7", "1237734238.0", MPFR_RNDN, 53, 64, 53, 931 "1.2377342379999998e9"); 932 check("-4.567291988483277e8", "1262857194.0", MPFR_RNDN, 53, 64, 53, 933 "8.0612799515167236e8"); 934 check("4.7719471752925262e7", "196089880.0", MPFR_RNDN, 53, 53, 53, 935 "2.4380935175292528e8"); 936 check("4.7719471752925262e7", "196089880.0", MPFR_RNDN, 53, 64, 53, 937 "2.4380935175292528e8"); 938 check("-1.716113812768534e-140", "1271212614.0", MPFR_RNDZ, 53, 64, 53, 939 "1.2712126139999998e9"); 940 check("-1.2927455200185474e-50", "1675676122.0", MPFR_RNDD, 53, 64, 53, 941 "1.6756761219999998e9"); 942 943 check53("1.22191250737771397120e+20", "948002822.0", MPFR_RNDN, 944 "122191250738719408128.0"); 945 check53("9966027674114492.0", "1780341389094537.0", MPFR_RNDN, 946 "11746369063209028.0"); 947 check53("2.99280481918991653800e+272", "5.34637717585790933424e+271", 948 MPFR_RNDN, "3.5274425367757071711e272"); 949 check_same(); 950 check53("6.14384195492641560499e-02", "-6.14384195401037683237e-02", 951 MPFR_RNDU, "9.1603877261370314499e-12"); 952 check53("1.16809465359248765399e+196", "7.92883212101990665259e+196", 953 MPFR_RNDU, "9.0969267746123943065e196"); 954 check53("3.14553393112021279444e-67", "3.14553401015952024126e-67", MPFR_RNDU, 955 "6.2910679412797336946e-67"); 956 957 check53("5.43885304644369509058e+185","-1.87427265794105342763e-57",MPFR_RNDN, 958 "5.4388530464436950905e185"); 959 check53("5.43885304644369509058e+185","-1.87427265794105342763e-57",MPFR_RNDZ, 960 "5.4388530464436944867e185"); 961 check53("5.43885304644369509058e+185","-1.87427265794105342763e-57",MPFR_RNDU, 962 "5.4388530464436950905e185"); 963 check53("5.43885304644369509058e+185","-1.87427265794105342763e-57",MPFR_RNDD, 964 "5.4388530464436944867e185"); 965 966 check2b("1.001010101110011000000010100101110010111001010000001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e358",187, 967 "-1.11100111001101100010001111111110101101110001000000000000000000000000000000000000000000e160",87, 968 "1.001010101110011000000010100101110010111001010000000111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111e358",178, 969 MPFR_RNDD); 970 check2b("-1.111100100011100111010101010101001010100100000111001000000000000000000e481",70, 971 "1.1111000110100011110101111110110010010000000110101000000000000000e481",65, 972 "-1.001010111111101011010000001100011101100101000000000000000000e472",61, 973 MPFR_RNDD); 974 check2b("1.0100010111010000100101000000111110011100011001011010000000000000000000000000000000e516",83, 975 "-1.1001111000100001011100000001001100110011110010111111000000e541",59, 976 "-1.1001111000100001011011110111000001001011100000011110100000110001110011010011000000000000000000000000000000000000000000000000e541",125, 977 MPFR_RNDZ); 978 check2b("-1.0010111100000100110001011011010000000011000111101000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e261",155, 979 "-1.00111110100011e239",15, 980 "-1.00101111000001001100101010101110001100110001111010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e261",159, 981 MPFR_RNDD); 982 check2b("-1.110111000011111011000000001001111101101001010100111000000000000000000000000e880",76, 983 "-1.1010010e-634",8, 984 "-1.11011100001111101100000000100111110110100101010011100000000000000000000000e880",75, 985 MPFR_RNDZ); 986 check2b("1.00100100110110101001010010101111000001011100100101010000000000000000000000000000e-530",81, 987 "-1.101101111100000111000011001010110011001011101001110100000e-908",58, 988 "1.00100100110110101001010010101111000001011100100101010e-530",54, 989 MPFR_RNDN); 990 check2b("1.0101100010010111101000000001000010010010011000111011000000000000000000000000000000000000000000000000000000000000000000e374",119, 991 "1.11100101100101e358",15, 992 "1.01011000100110011000010110100100100100100110001110110000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e374",150, 993 MPFR_RNDZ); 994 check2b("-1.10011001000010100000010100100110110010011111101111000000000000000000000000000000000000000000000000000000000000000000e-172",117, 995 "1.111011100000101010110000100100110100100001001000011100000000e-173",61, 996 "-1.0100010000001001010110011011101001001011101011110001000000000000000e-173",68, 997 MPFR_RNDZ); 998 check2b("-1.011110000111101011100001100110100011100101000011011000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-189",175, 999 "1.1e631",2, 1000 "1.011111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111e631",115, 1001 MPFR_RNDZ); 1002 check2b("-1.101011101001011101100011001000001100010100001101011000000000000000000000000000000000000000000e-449",94, 1003 "-1.01111101010111000011000110011101000111001100110111100000000000000e-429",66, 1004 "-1.01111101010111000100110010000110100100101111111111101100010100001101011000000000000000000000000000000000000000e-429",111, 1005 MPFR_RNDU); 1006 check2b("-1.101011101001011101100011001000001100010100001101011000000000000000000000000000000000000000000e-449",94, 1007 "-1.01111101010111000011000110011101000111001100110111100000000000000e-429",66, 1008 "-1.01111101010111000100110010000110100100101111111111101100010100001101011000000000000000000000000000000000000000e-429",111, 1009 MPFR_RNDD); 1010 check2b("-1.1001000011101000110000111110010100100101110101111100000000000000000000000000000000000000000000000000000000e-72",107, 1011 "-1.001100011101100100010101101010101011010010111111010000000000000000000000000000e521",79, 1012 "-1.00110001110110010001010110101010101101001011111101000000000000000000000000000000000000000000000001e521",99, 1013 MPFR_RNDD); 1014 check2b("-1.01010001111000000101010100100100110101011011100001110000000000e498",63, 1015 "1.010000011010101111000100111100011100010101011110010100000000000e243",64, 1016 "-1.010100011110000001010101001001001101010110111000011100000000000e498",64, 1017 MPFR_RNDN); 1018 check2b("1.00101100010101000011010000011000111111011110010111000000000000000000000000000000000000000000000000000000000e178",108, 1019 "-1.10101101010101000110011011111001001101111111110000100000000e160",60, 1020 "1.00101100010100111100100011000011111001000010011101110010000000001111100000000000000000000000000000000000e178",105, 1021 MPFR_RNDN); 1022 check2b("1.00110011010100111110011010110100111101110101100100110000000000000000000000000000000000000000000000e559",99, 1023 "-1.011010110100111011100110100110011100000000111010011000000000000000e559",67, 1024 "-1.101111111101011111111111001001100100011100001001100000000000000000000000000000000000000000000e556",94, 1025 MPFR_RNDU); 1026 check2b("-1.100000111100101001100111011100011011000001101001111100000000000000000000000000e843",79, 1027 "-1.1101101010110000001001000100001100110011000110110111000000000000000000000000000000000000000000e414",95, 1028 "-1.1000001111001010011001110111000110110000011010100000e843",53, 1029 MPFR_RNDD); 1030 check2b("-1.110110010110100010100011000110111001010000010111110000000000e-415",61, 1031 "-1.0000100101100001111100110011111111110100011101101011000000000000000000e751",71, 1032 "-1.00001001011000011111001100111111111101000111011010110e751",54, 1033 MPFR_RNDN); 1034 check2b("-1.1011011011110001001101010101001000010100010110111101000000000000000000000e258",74, 1035 "-1.00011100010110110101001011000100100000100010101000010000000000000000000000000000000000000000000000e268",99, 1036 "-1.0001110011001001000011110001000111010110101011110010011011110100000000000000000000000000000000000000e268",101, 1037 MPFR_RNDD); 1038 check2b("-1.1011101010011101011000000100100110101101101110000001000000000e629",62, 1039 "1.111111100000011100100011100000011101100110111110111000000000000000000000000000000000000000000e525",94, 1040 "-1.101110101001110101100000010010011010110110111000000011111111111111111111111111111111111111111111111111101e629",106, 1041 MPFR_RNDD); 1042 check2b("1.111001000010001100010000001100000110001011110111011000000000000000000000000000000000000e152",88, 1043 "1.111110111001100100000100111111010111000100111111001000000000000000e152",67, 1044 "1.1110111111011110000010101001011011101010000110110100e153",53, 1045 MPFR_RNDN); 1046 check2b("1.000001100011110010110000110100001010101101111011110100e696",55, 1047 "-1.1011001111011100100001011110100101010101110111010101000000000000000000000000000000000000000000000000000000000000e730",113, 1048 "-1.1011001111011100100001011110100100010100010011100010e730",53, 1049 MPFR_RNDN); 1050 check2b("-1.11010111100001001111000001110101010010001111111001100000000000000000000000000000000000000000000000000000000000e530",111, 1051 "1.01110100010010000000010110111101011101000001111101100000000000000000000000000000000000000000000000e530",99, 1052 "-1.1000110011110011101010101101111101010011011111000000000000000e528",62, 1053 MPFR_RNDD); 1054 check2b("-1.0001100010010100111101101011101000100100010011100011000000000000000000000000000000000000000000000000000000000e733",110, 1055 "-1.001000000111110010100101010100110111001111011011001000000000000000000000000000000000000000000000000000000000e710",109, 1056 "-1.000110001001010011111000111110110001110110011000110110e733",55, 1057 MPFR_RNDN); 1058 check2b("-1.1101011110000100111100000111010101001000111111100110000000000000000000000e530",74, 1059 "1.01110100010010000000010110111101011101000001111101100000000000000000000000000000000000000000000000000000000000e530",111, 1060 "-1.10001100111100111010101011011111010100110111110000000000000000000000000000e528",75, 1061 MPFR_RNDU); 1062 check2b("1.00110011010100111110011010110100111101110101100100110000000000000000000000000000000000000000000000e559",99, 1063 "-1.011010110100111011100110100110011100000000111010011000000000000000e559",67, 1064 "-1.101111111101011111111111001001100100011100001001100000000000000000000000000000000000000000000e556",94, 1065 MPFR_RNDU); 1066 check2b("-1.100101111110110000000110111111011010011101101111100100000000000000e-624",67, 1067 "1.10111010101110100000010110101000000000010011100000100000000e-587",60, 1068 "1.1011101010111010000001011010011111110100011110001011111111001000000100101100010010000011100000000000000000000e-587",110, 1069 MPFR_RNDU); 1070 check2b("-1.10011001000010100000010100100110110010011111101111000000000000000000000000000000000000000000000000000000000000000000e-172",117, 1071 "1.111011100000101010110000100100110100100001001000011100000000e-173",61, 1072 "-1.0100010000001001010110011011101001001011101011110001000000000000000e-173",68, 1073 MPFR_RNDZ); 1074 check2b("1.1000111000110010101001010011010011101100010110001001000000000000000000000000000000000000000000000000e167",101, 1075 "1.0011110010000110000000101100100111000001110110110000000000000000000000000e167",74, 1076 "1.01100101010111000101001111111111010101110001100111001000000000000000000000000000000000000000000000000000e168",105, 1077 MPFR_RNDZ); 1078 check2b("1.100101111111110010100101110111100001110000100001010000000000000000000000000000000000000000000000e808",97, 1079 "-1.1110011001100000100000111111110000110010100111001011000000000000000000000000000000e807",83, 1080 "1.01001001100110001100011111000000000001011010010111010000000000000000000000000000000000000000000e807",96, 1081 MPFR_RNDN); 1082 check2b("1e128",128, 1083 "1e0",128, 1084 "100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001e0",256, 1085 MPFR_RNDN); 1086 1087 /* Checking double precision (53 bits) */ 1088 check53("-8.22183238641455905806e-19", "7.42227178769761587878e-19",MPFR_RNDD, 1089 "-7.9956059871694317927e-20"); 1090 check53("5.82106394662028628236e+234","-5.21514064202368477230e+89",MPFR_RNDD, 1091 "5.8210639466202855763e234"); 1092 check53("5.72931679569871602371e+122","-5.72886070363264321230e+122", 1093 MPFR_RNDN, "4.5609206607281141508e118"); 1094 check53("-5.09937369394650450820e+238", "2.70203299854862982387e+250", 1095 MPFR_RNDD, "2.7020329985435301323e250"); 1096 check53("-2.96695924472363684394e+27", "1.22842938251111500000e+16",MPFR_RNDD, 1097 "-2.96695924471135255027e27"); 1098 check53("1.74693641655743793422e-227", "-7.71776956366861843469e-229", 1099 MPFR_RNDN, "1.669758720920751867e-227"); 1100 /* x = -7883040437021647.0; for (i=0; i<468; i++) x = x / 2.0;*/ 1101 check53("-1.03432206392780011159e-125", "1.30127034799251347548e-133", 1102 MPFR_RNDN, 1103 "-1.0343220509150965661100887242027378881805094180354e-125"); 1104 check53("1.05824655795525779205e+71", "-1.06022698059744327881e+71",MPFR_RNDZ, 1105 "-1.9804226421854867632e68"); 1106 check53("-5.84204911040921732219e+240", "7.26658169050749590763e+240", 1107 MPFR_RNDD, "1.4245325800982785854e240"); 1108 check53("1.00944884131046636376e+221","2.33809162651471520268e+215",MPFR_RNDN, 1109 "1.0094511794020929787e221"); 1110 /*x = 7045852550057985.0; for (i=0; i<986; i++) x = x / 2.0;*/ 1111 check53("4.29232078932667367325e-278", 1112 "1.0773525047389793833221116707010783793203080117586e-281" 1113 , MPFR_RNDU, "4.2933981418314132787e-278"); 1114 check53("5.27584773801377058681e-80", "8.91207657803547196421e-91", MPFR_RNDN, 1115 "5.2758477381028917269e-80"); 1116 check53("2.99280481918991653800e+272", "5.34637717585790933424e+271", 1117 MPFR_RNDN, "3.5274425367757071711e272"); 1118 check53("4.67302514390488041733e-184", "2.18321376145645689945e-190", 1119 MPFR_RNDN, "4.6730273271186420541e-184"); 1120 check53("5.57294120336300389254e+71", "2.60596167942024924040e+65", MPFR_RNDZ, 1121 "5.5729438093246831053e71"); 1122 check53("6.6052588496951015469e24", "4938448004894539.0", MPFR_RNDU, 1123 "6.6052588546335505068e24"); 1124 check53("1.23056185051606761523e-190", "1.64589756643433857138e-181", 1125 MPFR_RNDU, "1.6458975676649006598e-181"); 1126 check53("2.93231171510175981584e-280", "3.26266919161341483877e-273", 1127 MPFR_RNDU, "3.2626694848445867288e-273"); 1128 check53("5.76707395945001907217e-58", "4.74752971449827687074e-51", MPFR_RNDD, 1129 "4.747530291205672325e-51"); 1130 check53("277363943109.0", "11.0", MPFR_RNDN, "277363943120.0"); 1131 check53("1.44791789689198883921e-140", "-1.90982880222349071284e-121", 1132 MPFR_RNDN, "-1.90982880222349071e-121"); 1133 1134 1135 /* tests for particular cases (Vincent Lefevre, 22 Aug 2001) */ 1136 check53("9007199254740992.0", "1.0", MPFR_RNDN, "9007199254740992.0"); 1137 check53("9007199254740994.0", "1.0", MPFR_RNDN, "9007199254740996.0"); 1138 check53("9007199254740992.0", "-1.0", MPFR_RNDN, "9007199254740991.0"); 1139 check53("9007199254740994.0", "-1.0", MPFR_RNDN, "9007199254740992.0"); 1140 check53("9007199254740996.0", "-1.0", MPFR_RNDN, "9007199254740996.0"); 1141 1142 check_overflow (); 1143 check_1111 (); 1144 check_1minuseps (); 1145 } 1146 1147 static void 1148 check_extreme (void) 1149 { 1150 mpfr_t u, v, w, x, y; 1151 int i, inex, r; 1152 1153 mpfr_inits2 (32, u, v, w, x, y, (mpfr_ptr) 0); 1154 mpfr_setmin (u, mpfr_get_emax ()); 1155 mpfr_setmax (v, mpfr_get_emin ()); 1156 mpfr_setmin (w, mpfr_get_emax () - 40); 1157 RND_LOOP (r) 1158 for (i = 0; i < 2; i++) 1159 { 1160 mpfr_add (x, u, v, (mpfr_rnd_t) r); 1161 mpfr_set_prec (y, 64); 1162 inex = mpfr_add (y, u, w, MPFR_RNDN); 1163 MPFR_ASSERTN (inex == 0); 1164 mpfr_prec_round (y, 32, (mpfr_rnd_t) r); 1165 /* for RNDF, the rounding directions are not necessarily consistent */ 1166 if (! mpfr_equal_p (x, y) && r != MPFR_RNDF) 1167 { 1168 printf ("Error in u + v (%s)\n", 1169 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 1170 printf ("u = "); 1171 mpfr_dump (u); 1172 printf ("v = "); 1173 mpfr_dump (v); 1174 printf ("Expected "); 1175 mpfr_dump (y); 1176 printf ("Got "); 1177 mpfr_dump (x); 1178 exit (1); 1179 } 1180 mpfr_neg (v, v, MPFR_RNDN); 1181 mpfr_neg (w, w, MPFR_RNDN); 1182 } 1183 mpfr_clears (u, v, w, x, y, (mpfr_ptr) 0); 1184 } 1185 1186 static void 1187 testall_rndf (mpfr_prec_t pmax) 1188 { 1189 mpfr_t a, b, c, d; 1190 mpfr_prec_t pa, pb, pc; 1191 mpfr_exp_t eb; 1192 1193 for (pa = MPFR_PREC_MIN; pa <= pmax; pa++) 1194 { 1195 mpfr_init2 (a, pa); 1196 mpfr_init2 (d, pa); 1197 for (pb = MPFR_PREC_MIN; pb <= pmax; pb++) 1198 { 1199 mpfr_init2 (b, pb); 1200 for (eb = 0; eb <= pmax + 3; eb ++) 1201 { 1202 mpfr_set_ui_2exp (b, 1, eb, MPFR_RNDN); 1203 while (mpfr_cmp_ui_2exp (b, 1, eb + 1) < 0) 1204 { 1205 for (pc = MPFR_PREC_MIN; pc <= pmax; pc++) 1206 { 1207 mpfr_init2 (c, pc); 1208 mpfr_set_ui (c, 1, MPFR_RNDN); 1209 while (mpfr_cmp_ui (c, 2) < 0) 1210 { 1211 mpfr_add (a, b, c, MPFR_RNDF); 1212 mpfr_add (d, b, c, MPFR_RNDD); 1213 if (!mpfr_equal_p (a, d)) 1214 { 1215 mpfr_add (d, b, c, MPFR_RNDU); 1216 if (!mpfr_equal_p (a, d)) 1217 { 1218 printf ("Error: mpfr_add(a,b,c,RNDF) does not " 1219 "match RNDD/RNDU\n"); 1220 printf ("b="); mpfr_dump (b); 1221 printf ("c="); mpfr_dump (c); 1222 printf ("a="); mpfr_dump (a); 1223 exit (1); 1224 } 1225 } 1226 mpfr_nextabove (c); 1227 } 1228 mpfr_clear (c); 1229 } 1230 mpfr_nextabove (b); 1231 } 1232 } 1233 mpfr_clear (b); 1234 } 1235 mpfr_clear (a); 1236 mpfr_clear (d); 1237 } 1238 } 1239 1240 static void 1241 test_rndf_exact (mp_size_t pmax) 1242 { 1243 mpfr_t a, b, c, d; 1244 mpfr_prec_t pa, pb, pc; 1245 mpfr_exp_t eb; 1246 1247 for (pa = MPFR_PREC_MIN; pa <= pmax; pa++) 1248 { 1249 /* only check pa mod GMP_NUMB_BITS = -2, -1, 0, 1, 2 */ 1250 if ((pa + 2) % GMP_NUMB_BITS > 4) 1251 continue; 1252 mpfr_init2 (a, pa); 1253 mpfr_init2 (d, pa); 1254 for (pb = MPFR_PREC_MIN; pb <= pmax; pb++) 1255 { 1256 if ((pb + 2) % GMP_NUMB_BITS > 4) 1257 continue; 1258 mpfr_init2 (b, pb); 1259 for (eb = 0; eb <= pmax + 3; eb ++) 1260 { 1261 mpfr_urandomb (b, RANDS); 1262 mpfr_mul_2exp (b, b, eb, MPFR_RNDN); 1263 for (pc = MPFR_PREC_MIN; pc <= pmax; pc++) 1264 { 1265 if ((pc + 2) % GMP_NUMB_BITS > 4) 1266 continue; 1267 mpfr_init2 (c, pc); 1268 mpfr_urandomb (c, RANDS); 1269 mpfr_add (a, b, c, MPFR_RNDF); 1270 mpfr_add (d, b, c, MPFR_RNDD); 1271 if (!mpfr_equal_p (a, d)) 1272 { 1273 mpfr_add (d, b, c, MPFR_RNDU); 1274 if (!mpfr_equal_p (a, d)) 1275 { 1276 printf ("Error: mpfr_add(a,b,c,RNDF) does not " 1277 "match RNDD/RNDU\n"); 1278 printf ("b="); mpfr_dump (b); 1279 printf ("c="); mpfr_dump (c); 1280 printf ("a="); mpfr_dump (a); 1281 exit (1); 1282 } 1283 } 1284 1285 /* now make the low bits from c match those from b */ 1286 mpfr_sub (c, b, d, MPFR_RNDN); 1287 mpfr_add (a, b, c, MPFR_RNDF); 1288 mpfr_add (d, b, c, MPFR_RNDD); 1289 if (!mpfr_equal_p (a, d)) 1290 { 1291 mpfr_add (d, b, c, MPFR_RNDU); 1292 if (!mpfr_equal_p (a, d)) 1293 { 1294 printf ("Error: mpfr_add(a,b,c,RNDF) does not " 1295 "match RNDD/RNDU\n"); 1296 printf ("b="); mpfr_dump (b); 1297 printf ("c="); mpfr_dump (c); 1298 printf ("a="); mpfr_dump (a); 1299 exit (1); 1300 } 1301 } 1302 1303 mpfr_clear (c); 1304 } 1305 } 1306 mpfr_clear (b); 1307 } 1308 mpfr_clear (a); 1309 mpfr_clear (d); 1310 } 1311 } 1312 1313 #define TEST_FUNCTION test_add 1314 #define TWO_ARGS 1315 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS) 1316 #include "tgeneric.c" 1317 1318 int 1319 main (int argc, char *argv[]) 1320 { 1321 tests_start_mpfr (); 1322 1323 usesp = 0; 1324 tests (); 1325 1326 #ifndef CHECK_EXTERNAL /* no need to check twice */ 1327 usesp = 1; 1328 tests (); 1329 #endif 1330 1331 test_rndf_exact (200); 1332 testall_rndf (7); 1333 check_extreme (); 1334 1335 test_generic (MPFR_PREC_MIN, 1000, 100); 1336 1337 tests_end_mpfr (); 1338 return 0; 1339 } 1340
Indexes created Thu Mar 26 00:23:24 UTC 2026