1 /* Test file for mpfr_sub1sp. 2 3 Copyright 2003-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 static void check_underflow (mpfr_prec_t p); 28 static void check_corner (mpfr_prec_t p); 29 30 static void 31 bug20170109 (void) 32 { 33 mpfr_t a, b, c; 34 35 mpfr_init2 (a, 111); 36 mpfr_init2 (b, 111); 37 mpfr_init2 (c, 111); 38 mpfr_set_str_binary (b, "0.110010010000111111011010101000100010000101101000110000100011010011000100110001100110001010001011100000001101110E1"); 39 mpfr_set_str_binary (c, "0.111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111E-63"); 40 mpfr_sub (a, b, c, MPFR_RNDN); 41 mpfr_set_str_binary (b, "0.110010010000111111011010101000100010000101101000110000100011001111000100110001100110001010001011100000001101110E1"); 42 MPFR_ASSERTN(mpfr_equal_p (a, b)); 43 mpfr_clear (a); 44 mpfr_clear (b); 45 mpfr_clear (c); 46 } 47 48 /* check mpfr_sub1sp1 when: 49 (1) p = GMP_NUMB_BITS-1, d = GMP_NUMB_BITS and bp[0] = MPFR_LIMB_HIGHBIT 50 (2) p = 2*GMP_NUMB_BITS-1, d = 2*GMP_NUMB_BITS and b = 1000...000 51 (3) p = 3*GMP_NUMB_BITS-1, d = 3*GMP_NUMB_BITS and b = 1000...000 52 */ 53 static void 54 test20170208 (void) 55 { 56 mpfr_t a, b, c; 57 int inex; 58 59 mpfr_inits2 (GMP_NUMB_BITS - 1, a, b, c, (mpfr_ptr) 0); 60 61 /* test (1) */ 62 mpfr_set_ui_2exp (b, 1, GMP_NUMB_BITS, MPFR_RNDN); 63 mpfr_set_ui_2exp (c, 1, 0, MPFR_RNDN); 64 inex = mpfr_sub (a, b, c, MPFR_RNDN); 65 /* b-c = 2^GMP_NUMB_BITS-1 which has GMP_NUMB_BITS bits, thus we should 66 round to 2^GMP_NUMB_BITS (even rule) */ 67 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, GMP_NUMB_BITS) == 0); 68 MPFR_ASSERTN(inex > 0); 69 inex = mpfr_sub1sp (a, b, c, MPFR_RNDN); 70 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, GMP_NUMB_BITS) == 0); 71 MPFR_ASSERTN(inex > 0); 72 73 mpfr_set_ui_2exp (c, 2, 0, MPFR_RNDN); 74 mpfr_nextbelow (c); 75 /* now c = 2 - 2^(1-GMP_NUMB_BITS) */ 76 inex = mpfr_sub (a, b, c, MPFR_RNDN); 77 /* b-c = 2^GMP_NUMB_BITS-2+2^(1-GMP_NUMB_BITS), which should 78 round to 2^GMP_NUMB_BITS-2. We check by directly inspecting the bit 79 field of a, since mpfr_cmp_ui might not work if unsigned long is shorter 80 than mp_limb_t, and we don't want to use mpfr_add_ui or mpfr_sub_ui 81 to construct the expected result. */ 82 MPFR_ASSERTN(MPFR_MANT(a)[0] == (mp_limb_t) -2); 83 MPFR_ASSERTN(MPFR_EXP(a) == GMP_NUMB_BITS); 84 MPFR_ASSERTN(inex < 0); 85 inex = mpfr_sub1sp (a, b, c, MPFR_RNDN); 86 MPFR_ASSERTN(MPFR_MANT(a)[0] == (mp_limb_t) -2); 87 MPFR_ASSERTN(MPFR_EXP(a) == GMP_NUMB_BITS); 88 MPFR_ASSERTN(inex < 0); 89 90 /* test (2) */ 91 mpfr_set_prec (a, 2 * GMP_NUMB_BITS - 1); 92 mpfr_set_prec (b, 2 * GMP_NUMB_BITS - 1); 93 mpfr_set_prec (c, 2 * GMP_NUMB_BITS - 1); 94 mpfr_set_ui_2exp (b, 1, 2 * GMP_NUMB_BITS, MPFR_RNDN); 95 mpfr_set_ui_2exp (c, 1, 0, MPFR_RNDN); 96 inex = mpfr_sub (a, b, c, MPFR_RNDN); 97 /* b-c = 2^(2*GMP_NUMB_BITS)-1 which has 2*GMP_NUMB_BITS bits, thus we should 98 round to 2^(2*GMP_NUMB_BITS) (even rule) */ 99 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, 2 * GMP_NUMB_BITS) == 0); 100 MPFR_ASSERTN(inex > 0); 101 inex = mpfr_sub1sp (a, b, c, MPFR_RNDN); 102 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, 2 * GMP_NUMB_BITS) == 0); 103 MPFR_ASSERTN(inex > 0); 104 105 mpfr_set_ui_2exp (c, 2, 0, MPFR_RNDN); 106 mpfr_nextbelow (c); 107 /* now c = 2 - 2^(1-2*GMP_NUMB_BITS) */ 108 inex = mpfr_sub (a, b, c, MPFR_RNDN); 109 /* b-c = 2^(2*GMP_NUMB_BITS)-2+2^(1-2*GMP_NUMB_BITS), which should 110 round to 2^(2*GMP_NUMB_BITS)-2. We check by directly inspecting the bit 111 field of a, since mpfr_cmp_ui might not work if unsigned long is shorter 112 than mp_limb_t, and we don't want to use mpfr_add_ui or mpfr_sub_ui 113 to construct the expected result. */ 114 MPFR_ASSERTN(MPFR_MANT(a)[1] == (mp_limb_t) -1); 115 MPFR_ASSERTN(MPFR_MANT(a)[0] == (mp_limb_t) -2); 116 MPFR_ASSERTN(MPFR_EXP(a) == 2 * GMP_NUMB_BITS); 117 MPFR_ASSERTN(inex < 0); 118 inex = mpfr_sub1sp (a, b, c, MPFR_RNDN); 119 MPFR_ASSERTN(MPFR_MANT(a)[1] == (mp_limb_t) -1); 120 MPFR_ASSERTN(MPFR_MANT(a)[0] == (mp_limb_t) -2); 121 MPFR_ASSERTN(MPFR_EXP(a) == 2 * GMP_NUMB_BITS); 122 MPFR_ASSERTN(inex < 0); 123 124 /* test (3) */ 125 mpfr_set_prec (a, 3 * GMP_NUMB_BITS - 1); 126 mpfr_set_prec (b, 3 * GMP_NUMB_BITS - 1); 127 mpfr_set_prec (c, 3 * GMP_NUMB_BITS - 1); 128 mpfr_set_ui_2exp (b, 1, 3 * GMP_NUMB_BITS, MPFR_RNDN); 129 mpfr_set_ui_2exp (c, 1, 0, MPFR_RNDN); 130 inex = mpfr_sub (a, b, c, MPFR_RNDN); 131 /* b-c = 2^(3*GMP_NUMB_BITS)-1 which has 3*GMP_NUMB_BITS bits, thus we should 132 round to 3^(2*GMP_NUMB_BITS) (even rule) */ 133 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, 3 * GMP_NUMB_BITS) == 0); 134 MPFR_ASSERTN(inex > 0); 135 inex = mpfr_sub1sp (a, b, c, MPFR_RNDN); 136 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, 3 * GMP_NUMB_BITS) == 0); 137 MPFR_ASSERTN(inex > 0); 138 139 mpfr_set_ui_2exp (c, 2, 0, MPFR_RNDN); 140 mpfr_nextbelow (c); 141 /* now c = 2 - 2^(1-3*GMP_NUMB_BITS) */ 142 inex = mpfr_sub (a, b, c, MPFR_RNDN); 143 /* b-c = 2^(3*GMP_NUMB_BITS)-2+2^(1-3*GMP_NUMB_BITS), which should 144 round to 2^(3*GMP_NUMB_BITS)-2. We check by directly inspecting the bit 145 field of a, since mpfr_cmp_ui might not work if unsigned long is shorter 146 than mp_limb_t, and we don't want to use mpfr_add_ui or mpfr_sub_ui 147 to construct the expected result. */ 148 MPFR_ASSERTN(MPFR_MANT(a)[2] == (mp_limb_t) -1); 149 MPFR_ASSERTN(MPFR_MANT(a)[1] == (mp_limb_t) -1); 150 MPFR_ASSERTN(MPFR_MANT(a)[0] == (mp_limb_t) -2); 151 MPFR_ASSERTN(MPFR_EXP(a) == 3 * GMP_NUMB_BITS); 152 MPFR_ASSERTN(inex < 0); 153 inex = mpfr_sub1sp (a, b, c, MPFR_RNDN); 154 MPFR_ASSERTN(MPFR_MANT(a)[2] == (mp_limb_t) -1); 155 MPFR_ASSERTN(MPFR_MANT(a)[1] == (mp_limb_t) -1); 156 MPFR_ASSERTN(MPFR_MANT(a)[0] == (mp_limb_t) -2); 157 MPFR_ASSERTN(MPFR_EXP(a) == 3 * GMP_NUMB_BITS); 158 MPFR_ASSERTN(inex < 0); 159 160 mpfr_clears (a, b, c, (mpfr_ptr) 0); 161 } 162 163 static void 164 compare_sub_sub1sp (void) 165 { 166 mpfr_t a, b, c, a_ref; 167 mpfr_prec_t p; 168 unsigned long d; 169 int i, inex_ref, inex; 170 int r; 171 172 for (p = 1; p <= 3*GMP_NUMB_BITS; p++) 173 { 174 mpfr_inits2 (p, a, b, c, a_ref, (mpfr_ptr) 0); 175 for (d = 0; d <= p + 2; d++) 176 { 177 /* EXP(b) - EXP(c) = d */ 178 for (i = 0; i < 4; i++) 179 { 180 /* for i even, b is the smallest number, for b odd the largest */ 181 mpfr_set_ui_2exp (b, 1, d, MPFR_RNDN); 182 if (i & 1) 183 { 184 mpfr_mul_2ui (b, b, 1, MPFR_RNDN); 185 mpfr_nextbelow (b); 186 } 187 mpfr_set_ui_2exp (c, 1, 0, MPFR_RNDN); 188 if (i & 2) 189 { 190 mpfr_mul_2ui (c, c, 1, MPFR_RNDN); 191 mpfr_nextbelow (c); 192 } 193 RND_LOOP_NO_RNDF (r) 194 { 195 /* increase the precision of b to ensure sub1sp is not used */ 196 mpfr_prec_round (b, p + 1, MPFR_RNDN); 197 inex_ref = mpfr_sub (a_ref, b, c, (mpfr_rnd_t) r); 198 inex = mpfr_prec_round (b, p, MPFR_RNDN); 199 MPFR_ASSERTN(inex == 0); 200 inex = mpfr_sub1sp (a, b, c, (mpfr_rnd_t) r); 201 if (inex != inex_ref) 202 { 203 printf ("mpfr_sub and mpfr_sub1sp differ for r=%s\n", 204 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 205 printf ("b="); mpfr_dump (b); 206 printf ("c="); mpfr_dump (c); 207 printf ("expected inex=%d and ", inex_ref); 208 mpfr_dump (a_ref); 209 printf ("got inex=%d and ", inex); 210 mpfr_dump (a); 211 exit (1); 212 } 213 MPFR_ASSERTN(mpfr_equal_p (a, a_ref)); 214 } 215 } 216 } 217 mpfr_clears (a, b, c, a_ref, (mpfr_ptr) 0); 218 } 219 } 220 221 static void 222 bug20171213 (void) 223 { 224 mpfr_t a, b, c; 225 226 mpfr_init2 (a, 127); 227 mpfr_init2 (b, 127); 228 mpfr_init2 (c, 127); 229 mpfr_set_str_binary (b, "0.1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E1"); 230 mpfr_set_str_binary (c, "0.1000011010111101100101100110101111111001011001010000110000000000000000000000000000000000000000000000000000000000000000000000000E-74"); 231 mpfr_sub (a, b, c, MPFR_RNDN); 232 mpfr_set_str_binary (b, "0.1111111111111111111111111111111111111111111111111111111111111111111111111101111001010000100110100110010100000001101001101011110E0"); 233 MPFR_ASSERTN(mpfr_equal_p (a, b)); 234 mpfr_clear (a); 235 mpfr_clear (b); 236 mpfr_clear (c); 237 } 238 239 /* generic test for bug20171213: 240 b = 1.0 with precision p 241 c = 0.1xxx110...0E-e with precision p, with e >= 1, such that the part 1xxx1 has 242 exactly p+1-e bits, thus b-c = 0.111..01... is exact on p+1 bits. 243 Due to the round-to-even rule, b-c should be rounded to 0.111..0. 244 */ 245 static void 246 bug20171213_gen (mpfr_prec_t pmax) 247 { 248 mpfr_prec_t p; 249 mpfr_exp_t e; 250 mpfr_t a, b, c, d; 251 252 for (p = MPFR_PREC_MIN; p <= pmax; p++) 253 { 254 for (e = 1; e < p; e++) 255 { 256 mpfr_init2 (a, p); 257 mpfr_init2 (b, p); 258 mpfr_init2 (c, p); 259 mpfr_init2 (d, p); 260 mpfr_set_ui (b, 1, MPFR_RNDN); 261 mpfr_set_ui_2exp (c, 1, p + 1 - e, MPFR_RNDN); /* c = 2^(p + 1 - e) */ 262 mpfr_sub_ui (c, c, 1, MPFR_RNDN); /* c = 2^(p + 1 - e) - 1 */ 263 mpfr_div_2ui (c, c, p + 1, MPFR_RNDN); /* c = 2^(-e) - 2^(-p-1) */ 264 /* the exact difference is 1 - 2^(-e) + 2^(-p-1) */ 265 mpfr_sub (a, b, c, MPFR_RNDN); 266 /* check that a = 1 - 2^(-e) */ 267 mpfr_set_ui_2exp (d, 1, e, MPFR_RNDN); /* b = 2^e */ 268 mpfr_sub_ui (d, d, 1, MPFR_RNDN); /* b = 2^e - 1 */ 269 mpfr_div_2ui (d, d, e, MPFR_RNDN); /* b = 1 - 2^(-e) */ 270 if (! mpfr_equal_p (a, d)) 271 { 272 printf ("bug20171213_gen failed for p=%ld, e=%ld\n", 273 (long) p, (long) e); 274 printf ("b="); mpfr_dump (b); 275 printf ("c="); mpfr_dump (c); 276 printf ("got a="); mpfr_dump (a); 277 printf ("expected d="); mpfr_dump (d); 278 exit (1); 279 } 280 mpfr_clear (a); 281 mpfr_clear (b); 282 mpfr_clear (c); 283 mpfr_clear (d); 284 } 285 } 286 } 287 288 static void 289 coverage (void) 290 { 291 mpfr_t a, b, c, d, u; 292 int inex; 293 294 /* coverage test in mpfr_sub1sp: case d=1, limb > MPFR_LIMB_HIGHBIT, RNDF 295 and also RNDZ */ 296 mpfr_init2 (a, 3 * GMP_NUMB_BITS); 297 mpfr_init2 (b, 3 * GMP_NUMB_BITS); 298 mpfr_init2 (c, 3 * GMP_NUMB_BITS); 299 mpfr_init2 (d, 3 * GMP_NUMB_BITS); 300 mpfr_init2 (u, 3 * GMP_NUMB_BITS); 301 mpfr_set_ui (b, 1, MPFR_RNDN); 302 mpfr_nextbelow (b); /* b = 1 - 2^(-p) */ 303 mpfr_set_prec (c, GMP_NUMB_BITS); 304 mpfr_set_ui_2exp (c, 1, -1, MPFR_RNDN); 305 mpfr_nextbelow (c); 306 mpfr_nextbelow (c); /* c = 1/2 - 2*2^(-GMP_NUMB_BITS-1) */ 307 mpfr_prec_round (c, 3 * GMP_NUMB_BITS, MPFR_RNDN); 308 mpfr_nextbelow (c); /* c = 1/2 - 2*2^(-GMP_NUMB_BITS-1) - 2^(-p-1) */ 309 /* b-c = c */ 310 mpfr_sub (a, b, c, MPFR_RNDF); 311 mpfr_sub (d, b, c, MPFR_RNDD); 312 mpfr_sub (u, b, c, MPFR_RNDU); 313 /* check a = d or u */ 314 MPFR_ASSERTN(mpfr_equal_p (a, d) || mpfr_equal_p (a, u)); 315 316 /* coverage test in mpfr_sub1sp: case d=p, RNDN, sb = 0, significand of b 317 is even but b<>2^e, (case 1e) */ 318 mpfr_set_prec (a, 3 * GMP_NUMB_BITS); 319 mpfr_set_prec (b, 3 * GMP_NUMB_BITS); 320 mpfr_set_prec (c, 3 * GMP_NUMB_BITS); 321 mpfr_set_ui (b, 1, MPFR_RNDN); 322 mpfr_nextabove (b); 323 mpfr_nextabove (b); 324 mpfr_set_ui_2exp (c, 1, -3 * GMP_NUMB_BITS, MPFR_RNDN); 325 inex = mpfr_sub (a, b, c, MPFR_RNDN); 326 MPFR_ASSERTN(inex > 0); 327 MPFR_ASSERTN(mpfr_equal_p (a, b)); 328 329 mpfr_clear (a); 330 mpfr_clear (b); 331 mpfr_clear (c); 332 mpfr_clear (d); 333 mpfr_clear (u); 334 } 335 336 /* bug in mpfr_sub1sp1n, made generic */ 337 static void 338 bug20180217 (mpfr_prec_t pmax) 339 { 340 mpfr_t a, b, c; 341 int inex; 342 mpfr_prec_t p; 343 344 for (p = MPFR_PREC_MIN; p <= pmax; p++) 345 { 346 mpfr_init2 (a, p); 347 mpfr_init2 (b, p); 348 mpfr_init2 (c, p); 349 mpfr_set_ui (b, 1, MPFR_RNDN); /* b = 1 */ 350 mpfr_set_ui_2exp (c, 1, -p-1, MPFR_RNDN); /* c = 2^(-p-1) */ 351 /* a - b = 1 - 2^(-p-1) and should be rounded to 1 (case 2f of 352 mpfr_sub1sp) */ 353 inex = mpfr_sub (a, b, c, MPFR_RNDN); 354 MPFR_ASSERTN(inex > 0); 355 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0); 356 /* check also when a=b */ 357 mpfr_set_ui (a, 1, MPFR_RNDN); 358 inex = mpfr_sub (a, a, c, MPFR_RNDN); 359 MPFR_ASSERTN(inex > 0); 360 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0); 361 /* and when a=c */ 362 mpfr_set_ui (b, 1, MPFR_RNDN); /* b = 1 */ 363 mpfr_set_ui_2exp (a, 1, -p-1, MPFR_RNDN); 364 inex = mpfr_sub (a, b, a, MPFR_RNDN); 365 MPFR_ASSERTN(inex > 0); 366 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0); 367 mpfr_clear (a); 368 mpfr_clear (b); 369 mpfr_clear (c); 370 } 371 } 372 373 /* bug in revision 12985 with tlog and GMP_CHECK_RANDOMIZE=1534111552615050 374 (introduced in revision 12242, does not affect the 4.0 branch) */ 375 static void 376 bug20180813 (void) 377 { 378 mpfr_t a, b, c; 379 380 mpfr_init2 (a, 194); 381 mpfr_init2 (b, 194); 382 mpfr_init2 (c, 194); 383 mpfr_set_str_binary (b, "0.10000111101000100000010000100010110111011100110100000101100111000010101000110110010101011101101011110110001000111001000010110010111010010100011011010100001010001110000101000010101110100110001000E7"); 384 mpfr_set_str_binary (c, "0.10000000000000000100001111010001000000100001000101101110111001101000001011001110000101010001101100101010111011010111101100010001110010000101100101110100101000110110101000010100011100001010000101E24"); 385 mpfr_sub (a, b, c, MPFR_RNDN); 386 mpfr_set_str_binary (b, "-0.11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111E23"); 387 MPFR_ASSERTN(mpfr_equal_p (a, b)); 388 mpfr_clear (a); 389 mpfr_clear (b); 390 mpfr_clear (c); 391 } 392 393 /* bug in revision 13599 with tatan and GMP_CHECK_RANDOMIZE=1567609230659336: 394 the values are equal, but the ternary value differs between sub1 and sub1sp 395 (bug introduced with mpfr_sub1sp2n, does not affect the 4.0 branch) */ 396 static void 397 bug20190904 (void) 398 { 399 mpfr_t a, b, c; 400 int ret; 401 402 mpfr_init2 (a, 128); 403 mpfr_init2 (b, 128); 404 mpfr_init2 (c, 128); 405 mpfr_set_str_binary (b, "0.11001001000011111101101010100010001000010110100011000010001101001100010011000110011000101000101110000000110111000001110011010001E1"); 406 mpfr_set_str_binary (c, "0.10010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010010000000000000000000000E-102"); 407 ret = mpfr_sub (a, b, c, MPFR_RNDN); 408 mpfr_set_str_binary (b, "0.11001001000011111101101010100010001000010110100011000010001101001100010011000110011000101000101101111111101111000001110011010001E1"); 409 MPFR_ASSERTN(mpfr_equal_p (a, b)); 410 MPFR_ASSERTN(ret > 0); 411 mpfr_clear (a); 412 mpfr_clear (b); 413 mpfr_clear (c); 414 } 415 416 int 417 main (void) 418 { 419 mpfr_prec_t p; 420 421 tests_start_mpfr (); 422 423 bug20190904 (); 424 bug20180813 (); 425 bug20180217 (1024); 426 coverage (); 427 compare_sub_sub1sp (); 428 test20170208 (); 429 bug20170109 (); 430 bug20171213 (); 431 bug20171213_gen (256); 432 check_special (); 433 for (p = MPFR_PREC_MIN ; p < 200 ; p++) 434 { 435 check_underflow (p); 436 check_random (p); 437 check_corner (p); 438 } 439 440 tests_end_mpfr (); 441 return 0; 442 } 443 444 #define STD_ERROR \ 445 do \ 446 { \ 447 printf("ERROR: for %s and p=%lu and i=%d:\nY=", \ 448 mpfr_print_rnd_mode ((mpfr_rnd_t) r), (unsigned long) p, i); \ 449 mpfr_dump (y); \ 450 printf ("Z="); mpfr_dump (z); \ 451 printf ("Expected: "); mpfr_dump (x2); \ 452 printf ("Got : "); mpfr_dump (x); \ 453 abort(); \ 454 } \ 455 while (0) 456 457 #define STD_ERROR2 \ 458 do \ 459 { \ 460 printf("ERROR: for %s and p=%lu and i=%d:\nY=", \ 461 mpfr_print_rnd_mode ((mpfr_rnd_t) r), (unsigned long) p, i); \ 462 mpfr_dump (y); \ 463 printf ("Z="); mpfr_dump (z); \ 464 printf ("Expected: "); mpfr_dump (x2); \ 465 printf ("Got : "); mpfr_dump (x); \ 466 printf ("Wrong inexact flag. Expected %d. Got %d\n", \ 467 inexact1, inexact2); \ 468 exit(1); \ 469 } \ 470 while (0) 471 472 static void 473 check_random (mpfr_prec_t p) 474 { 475 mpfr_t x,y,z,x2; 476 int r; 477 int i, inexact1, inexact2; 478 479 mpfr_inits2 (p, x, y, z, x2, (mpfr_ptr) 0); 480 481 for (i = 0 ; i < 500 ; i++) 482 { 483 mpfr_urandomb (y, RANDS); 484 mpfr_urandomb (z, RANDS); 485 if (MPFR_IS_PURE_FP(y) && MPFR_IS_PURE_FP(z)) 486 for(r = 0 ; r < MPFR_RND_MAX ; r++) 487 { 488 if (r == MPFR_RNDF) 489 continue; /* mpfr_sub1 and mpfr_sub1sp could differ, 490 and inexact makes no sense */ 491 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 492 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 493 if (mpfr_cmp(x, x2)) 494 STD_ERROR; 495 if (inexact1 != inexact2) 496 STD_ERROR2; 497 } 498 } 499 500 mpfr_clears (x, y, z, x2, (mpfr_ptr) 0); 501 } 502 503 static void 504 check_special (void) 505 { 506 mpfr_t x,y,z,x2; 507 int r; 508 mpfr_prec_t p; 509 int i = -1, inexact1, inexact2; 510 mpfr_exp_t es; 511 512 mpfr_inits (x, y, z, x2, (mpfr_ptr) 0); 513 514 for (r = 0 ; r < MPFR_RND_MAX ; r++) 515 { 516 if (r == MPFR_RNDF) 517 continue; /* comparison between sub1 and sub1sp makes no sense here */ 518 519 p = 53; 520 mpfr_set_prec(x, 53); 521 mpfr_set_prec(x2, 53); 522 mpfr_set_prec(y, 53); 523 mpfr_set_prec(z, 53); 524 525 mpfr_set_str_binary (y, 526 "0.10110111101101110010010010011011000001101101011011001E31"); 527 528 mpfr_sub1sp (x, y, y, (mpfr_rnd_t) r); 529 if (mpfr_cmp_ui(x, 0)) 530 { 531 printf("Error for x-x with p=%lu. Expected 0. Got:", 532 (unsigned long) p); 533 mpfr_dump (x); 534 exit(1); 535 } 536 537 mpfr_set(z, y, (mpfr_rnd_t) r); 538 mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 539 if (mpfr_cmp_ui(x, 0)) 540 { 541 printf("Error for x-y with y=x and p=%lu. Expected 0. Got:", 542 (unsigned long) p); 543 mpfr_dump (x); 544 exit(1); 545 } 546 /* diff = 0 */ 547 mpfr_set_str_binary (y, 548 "0.10110111101101110010010010011011001001101101011011001E31"); 549 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 550 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 551 if (mpfr_cmp(x, x2)) 552 STD_ERROR; 553 if (inexact1 != inexact2) 554 STD_ERROR2; 555 556 /* Diff = 1 */ 557 mpfr_set_str_binary (y, 558 "0.10110111101101110010010010011011000001101101011011001E30"); 559 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 560 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 561 if (mpfr_cmp(x, x2)) 562 STD_ERROR; 563 if (inexact1 != inexact2) 564 STD_ERROR2; 565 566 /* Diff = 2 */ 567 mpfr_set_str_binary (y, 568 "0.10110111101101110010010010011011000101101101011011001E32"); 569 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 570 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 571 if (mpfr_cmp(x, x2)) 572 STD_ERROR; 573 if (inexact1 != inexact2) 574 STD_ERROR2; 575 576 /* Diff = 32 */ 577 mpfr_set_str_binary (y, 578 "0.10110111101101110010010010011011000001101101011011001E63"); 579 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 580 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 581 if (mpfr_cmp(x, x2)) 582 STD_ERROR; 583 if (inexact1 != inexact2) 584 STD_ERROR2; 585 586 /* Diff = 52 */ 587 mpfr_set_str_binary (y, 588 "0.10110111101101110010010010011011010001101101011011001E83"); 589 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 590 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 591 if (mpfr_cmp(x, x2)) 592 STD_ERROR; 593 if (inexact1 != inexact2) 594 STD_ERROR2; 595 596 /* Diff = 53 */ 597 mpfr_set_str_binary (y, 598 "0.10110111101101110010010010011111000001101101011011001E31"); 599 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 600 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 601 if (mpfr_cmp(x, x2)) 602 STD_ERROR; 603 if (inexact1 != inexact2) 604 STD_ERROR2; 605 606 /* Diff > 200 */ 607 mpfr_set_str_binary (y, 608 "0.10110111101101110010010010011011000001101101011011001E331"); 609 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 610 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 611 if (mpfr_cmp(x, x2)) 612 STD_ERROR; 613 if (inexact1 != inexact2) 614 STD_ERROR2; 615 616 mpfr_set_str_binary (y, 617 "0.10000000000000000000000000000000000000000000000000000E31"); 618 mpfr_set_str_binary (z, 619 "0.11111111111111111111111111111111111111111111111111111E30"); 620 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 621 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 622 if (mpfr_cmp(x, x2)) 623 STD_ERROR; 624 if (inexact1 != inexact2) 625 STD_ERROR2; 626 627 mpfr_set_str_binary (y, 628 "0.10000000000000000000000000000000000000000000000000000E31"); 629 mpfr_set_str_binary (z, 630 "0.11111111111111111111111111111111111111111111111111111E29"); 631 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 632 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 633 if (mpfr_cmp(x, x2)) 634 STD_ERROR; 635 if (inexact1 != inexact2) 636 STD_ERROR2; 637 638 mpfr_set_str_binary (y, 639 "0.10000000000000000000000000000000000000000000000000000E52"); 640 mpfr_set_str_binary (z, 641 "0.10000000000010000000000000000000000000000000000000000E00"); 642 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 643 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 644 if (mpfr_cmp(x, x2)) 645 STD_ERROR; 646 if (inexact1 != inexact2) 647 STD_ERROR2; 648 649 mpfr_set_str_binary (y, 650 "0.11100000000000000000000000000000000000000000000000000E53"); 651 mpfr_set_str_binary (z, 652 "0.10000000000000000000000000000000000000000000000000000E00"); 653 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 654 inexact2 = mpfr_sub1sp(z, y, z, (mpfr_rnd_t) r); 655 mpfr_set(x, z, (mpfr_rnd_t) r); 656 if (mpfr_cmp(x, x2)) 657 STD_ERROR; 658 if (inexact1 != inexact2) 659 STD_ERROR2; 660 661 mpfr_set_str_binary (y, 662 "0.10000000000000000000000000000000000000000000000000000E53"); 663 mpfr_set_str_binary (z, 664 "0.10100000000000000000000000000000000000000000000000000E00"); 665 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 666 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 667 if (mpfr_cmp(x, x2)) 668 STD_ERROR; 669 if (inexact1 != inexact2) 670 STD_ERROR2; 671 672 mpfr_set_str_binary (y, 673 "0.10000000000000000000000000000000000000000000000000000E54"); 674 mpfr_set_str_binary (z, 675 "0.10100000000000000000000000000000000000000000000000000E00"); 676 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 677 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 678 if (mpfr_cmp(x, x2)) 679 STD_ERROR; 680 if (inexact1 != inexact2) 681 STD_ERROR2; 682 683 p = 63; 684 mpfr_set_prec(x, p); 685 mpfr_set_prec(x2, p); 686 mpfr_set_prec(y, p); 687 mpfr_set_prec(z, p); 688 mpfr_set_str_binary (y, 689 "0.100000000000000000000000000000000000000000000000000000000000000E62"); 690 mpfr_set_str_binary (z, 691 "0.110000000000000000000000000000000000000000000000000000000000000E00"); 692 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 693 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 694 if (mpfr_cmp(x, x2)) 695 STD_ERROR; 696 if (inexact1 != inexact2) 697 STD_ERROR2; 698 699 p = 64; 700 mpfr_set_prec(x, 64); 701 mpfr_set_prec(x2, 64); 702 mpfr_set_prec(y, 64); 703 mpfr_set_prec(z, 64); 704 705 mpfr_set_str_binary (y, 706 "0.1100000000000000000000000000000000000000000000000000000000000000E31"); 707 mpfr_set_str_binary (z, 708 "0.1111111111111111111111111110000000000000000000000000011111111111E29"); 709 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 710 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 711 if (mpfr_cmp(x, x2)) 712 STD_ERROR; 713 if (inexact1 != inexact2) 714 STD_ERROR2; 715 716 mpfr_set_str_binary (y, 717 "0.1000000000000000000000000000000000000000000000000000000000000000E63"); 718 mpfr_set_str_binary (z, 719 "0.1011000000000000000000000000000000000000000000000000000000000000E00"); 720 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 721 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 722 if (mpfr_cmp(x, x2)) 723 STD_ERROR; 724 if (inexact1 != inexact2) 725 STD_ERROR2; 726 727 mpfr_set_str_binary (y, 728 "0.1000000000000000000000000000000000000000000000000000000000000000E63"); 729 mpfr_set_str_binary (z, 730 "0.1110000000000000000000000000000000000000000000000000000000000000E00"); 731 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 732 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 733 if (mpfr_cmp(x, x2)) 734 STD_ERROR; 735 if (inexact1 != inexact2) 736 STD_ERROR2; 737 738 mpfr_set_str_binary (y, 739 "0.10000000000000000000000000000000000000000000000000000000000000E63"); 740 mpfr_set_str_binary (z, 741 "0.10000000000000000000000000000000000000000000000000000000000000E00"); 742 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 743 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 744 if (mpfr_cmp(x, x2)) 745 STD_ERROR; 746 if (inexact1 != inexact2) 747 STD_ERROR2; 748 749 mpfr_set_str_binary (y, 750 "0.1000000000000000000000000000000000000000000000000000000000000000E64"); 751 mpfr_set_str_binary (z, 752 "0.1010000000000000000000000000000000000000000000000000000000000000E00"); 753 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 754 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 755 if (mpfr_cmp(x, x2)) 756 STD_ERROR; 757 if (inexact1 != inexact2) 758 STD_ERROR2; 759 760 MPFR_SET_NAN(x); 761 MPFR_SET_NAN(x2); 762 mpfr_set_str_binary (y, 763 "0.1000000000000000000000000000000000000000000000000000000000000000" 764 "E-1073741823"); 765 mpfr_set_str_binary (z, 766 "0.1100000000000000000000000000000000000000000000000000000000000000" 767 "E-1073741823"); 768 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 769 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 770 if (mpfr_cmp(x, x2)) 771 STD_ERROR; 772 if (inexact1 != inexact2) 773 STD_ERROR2; 774 775 p = 9; 776 mpfr_set_prec(x, p); 777 mpfr_set_prec(x2, p); 778 mpfr_set_prec(y, p); 779 mpfr_set_prec(z, p); 780 781 mpfr_set_str_binary (y, "0.100000000E1"); 782 mpfr_set_str_binary (z, "0.100000000E-8"); 783 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 784 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 785 if (mpfr_cmp(x, x2)) 786 STD_ERROR; 787 if (inexact1 != inexact2) 788 STD_ERROR2; 789 790 p = 34; 791 mpfr_set_prec(x, p); 792 mpfr_set_prec(x2, p); 793 mpfr_set_prec(y, p); 794 mpfr_set_prec(z, p); 795 796 mpfr_set_str_binary (y, "-0.1011110000111100010111011100110100E-18"); 797 mpfr_set_str_binary (z, "0.1000101010110011010101011110000000E-14"); 798 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 799 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 800 if (mpfr_cmp(x, x2)) 801 STD_ERROR; 802 if (inexact1 != inexact2) 803 STD_ERROR2; 804 805 p = 124; 806 mpfr_set_prec(x, p); 807 mpfr_set_prec(x2, p); 808 mpfr_set_prec(y, p); 809 mpfr_set_prec(z, p); 810 811 mpfr_set_str_binary (y, 812 "0.1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E1"); 813 mpfr_set_str_binary (z, 814 "0.1011111000100111000011001000011101010101101100101010101001000001110100001101110110001110111010000011101001100010111110001100E-31"); 815 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 816 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 817 if (mpfr_cmp(x, x2)) 818 STD_ERROR; 819 if (inexact1 != inexact2) 820 STD_ERROR2; 821 822 p = 288; 823 mpfr_set_prec(x, p); 824 mpfr_set_prec(x2, p); 825 mpfr_set_prec(y, p); 826 mpfr_set_prec(z, p); 827 828 mpfr_set_str_binary (y, 829 "0.111000110011000001000111101010111011110011101001101111111110000011100101000001001010110010101010011001010100000001110011110001010101101010001011101110100100001011110100110000101101100011010001001011011010101010000010001101001000110010010111111011110001111101001000101101001100101100101000E80"); 830 mpfr_set_str_binary (z, 831 "-0.100001111111101001011010001100110010100111001110000110011101001011010100001000000100111011010110110010000000000010101101011000010000110001110010100001100101011100100100001011000100011110000001010101000100011101001000010111100000111000111011001000100100011000100000010010111000000100100111E-258"); 832 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 833 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 834 if (mpfr_cmp(x, x2)) 835 STD_ERROR; 836 if (inexact1 != inexact2) 837 STD_ERROR2; 838 839 p = 85; 840 mpfr_set_prec(x, p); 841 mpfr_set_prec(x2, p); 842 mpfr_set_prec(y, p); 843 mpfr_set_prec(z, p); 844 845 mpfr_set_str_binary (y, 846 "0.1111101110100110110110100010101011101001100010100011110110110010010011101100101111100E-4"); 847 mpfr_set_str_binary (z, 848 "0.1111101110100110110110100010101001001000011000111000011101100101110100001110101010110E-4"); 849 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 850 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 851 if (mpfr_cmp(x, x2)) 852 STD_ERROR; 853 if (inexact1 != inexact2) 854 STD_ERROR2; 855 856 p = 64; 857 mpfr_set_prec(x, p); mpfr_set_prec(x2, p); 858 mpfr_set_prec(y, p); mpfr_set_prec(z, p); 859 860 mpfr_set_str_binary (y, 861 "0.11000000000000000000000000000000" 862 "00000000000000000000000000000000E1"); 863 mpfr_set_str_binary (z, 864 "0.10000000000000000000000000000000" 865 "00000000000000000000000000000001E0"); 866 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 867 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 868 if (mpfr_cmp(x, x2)) 869 STD_ERROR; 870 if (inexact1 != inexact2) 871 STD_ERROR2; 872 873 mpfr_set_str_binary (y, 874 "0.11000000000000000000000000000000" 875 "000000000000000000000000000001E1"); 876 mpfr_set_str_binary (z, 877 "0.10000000000000000000000000000000" 878 "00000000000000000000000000000001E0"); 879 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 880 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 881 if (mpfr_cmp(x, x2)) 882 STD_ERROR; 883 if (inexact1 != inexact2) 884 STD_ERROR2; 885 886 es = mpfr_get_emin (); 887 set_emin (-1024); 888 889 mpfr_set_str_binary (y, 890 "0.10000000000000000000000000000000" 891 "000000000000000000000000000000E-1023"); 892 mpfr_set_str_binary (z, 893 "0.10000000000000000000000000000000" 894 "00000000000000000000000000000001E-1023"); 895 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 896 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 897 if (mpfr_cmp(x, x2)) 898 STD_ERROR; 899 if (inexact1 != inexact2) 900 STD_ERROR2; 901 902 mpfr_set_str_binary (y, 903 "0.10000000000000000000000000000000" 904 "000000000000000000000000000000E-1023"); 905 mpfr_set_str_binary (z, 906 "0.1000000000000000000000000000000" 907 "000000000000000000000000000000E-1023"); 908 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r); 909 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r); 910 if (mpfr_cmp(x, x2)) 911 STD_ERROR; 912 if (inexact1 != inexact2) 913 STD_ERROR2; 914 915 set_emin (es); 916 } 917 918 mpfr_clears (x, y, z, x2, (mpfr_ptr) 0); 919 } 920 921 static void 922 check_underflow (mpfr_prec_t p) 923 { 924 mpfr_t x, y, z; 925 int inexact; 926 927 mpfr_inits2 (p, x, y, z, (mpfr_ptr) 0); 928 929 if (p >= 2) /* we need p >= 2 so that 3 is exact */ 930 { 931 mpfr_set_ui_2exp (y, 4, mpfr_get_emin () - 2, MPFR_RNDN); 932 mpfr_set_ui_2exp (z, 3, mpfr_get_emin () - 2, MPFR_RNDN); 933 inexact = mpfr_sub (x, y, z, MPFR_RNDN); 934 if (inexact >= 0 || (mpfr_cmp_ui (x, 0) != 0)) 935 { 936 printf ("4*2^(emin-2) - 3*2^(emin-2) with RNDN failed for p=%ld\n", 937 (long) p); 938 printf ("Expected inexact < 0 with x=0\n"); 939 printf ("Got inexact=%d with x=", inexact); 940 mpfr_dump (x); 941 exit (1); 942 } 943 } 944 945 if (p >= 3) /* we need p >= 3 so that 5 is exact */ 946 { 947 mpfr_set_ui_2exp (y, 5, mpfr_get_emin () - 2, MPFR_RNDN); 948 mpfr_set_ui_2exp (z, 4, mpfr_get_emin () - 2, MPFR_RNDN); 949 inexact = mpfr_sub (x, y, z, MPFR_RNDN); 950 if (inexact >= 0 || (mpfr_cmp_ui (x, 0) != 0)) 951 { 952 printf ("5*2^(emin-2) - 4*2^(emin-2) with RNDN failed for p=%ld\n", 953 (long) p); 954 printf ("Expected inexact < 0 with x=0\n"); 955 printf ("Got inexact=%d with x=", inexact); 956 mpfr_dump (x); 957 exit (1); 958 } 959 } 960 961 mpfr_clears (x, y, z, (mpfr_ptr) 0); 962 } 963 964 /* check corner cases of mpfr_sub1sp in case d = 1 and limb = MPFR_LIMB_HIGHBIT */ 965 static void 966 check_corner (mpfr_prec_t p) 967 { 968 mpfr_t x, y, z; 969 mpfr_exp_t e; 970 int inex, odd; 971 972 if (p < 4) /* ensures that the initial value of z is > 1 below */ 973 return; 974 975 mpfr_inits2 (p, x, y, z, (mpfr_ptr) 0); 976 mpfr_const_pi (y, MPFR_RNDN); 977 mpfr_set_ui (z, 2, MPFR_RNDN); 978 inex = mpfr_sub (z, y, z, MPFR_RNDN); /* z is near pi-2, thus y-z is near 2 */ 979 MPFR_ASSERTN(inex == 0); 980 for (e = 0; e < p; e++) 981 { 982 /* add 2^(-e) to z */ 983 mpfr_mul_2ui (z, z, e, MPFR_RNDN); 984 inex = mpfr_add_ui (z, z, 1, MPFR_RNDN); 985 MPFR_ASSERTN(inex == 0); 986 mpfr_div_2ui (z, z, e, MPFR_RNDN); 987 988 /* compute x = y - z which should be exact, near 2-2^(-e) */ 989 inex = mpfr_sub (x, y, z, MPFR_RNDN); 990 MPFR_ASSERTN(inex == 0); 991 MPFR_ASSERTN(mpfr_get_exp (x) == 1); 992 993 /* restore initial z */ 994 mpfr_mul_2ui (z, z, e, MPFR_RNDN); 995 inex = mpfr_sub_ui (z, z, 1, MPFR_RNDN); 996 MPFR_ASSERTN(inex == 0); 997 mpfr_div_2ui (z, z, e, MPFR_RNDN); 998 999 /* subtract 2^(-e) to z */ 1000 mpfr_mul_2ui (z, z, e, MPFR_RNDN); 1001 inex = mpfr_sub_ui (z, z, 1, MPFR_RNDN); 1002 MPFR_ASSERTN(inex == 0); 1003 mpfr_div_2ui (z, z, e, MPFR_RNDN); 1004 1005 /* ensure last significant bit of z is 0 so that y-z is exact */ 1006 odd = mpfr_min_prec (z) == p; 1007 if (odd) /* add one ulp to z */ 1008 mpfr_nextabove (z); 1009 1010 /* compute x = y - z which should be exact, near 2+2^(-e) */ 1011 inex = mpfr_sub (x, y, z, MPFR_RNDN); 1012 MPFR_ASSERTN(inex == 0); 1013 MPFR_ASSERTN(mpfr_get_exp (x) == 2); 1014 1015 /* restore initial z */ 1016 if (odd) 1017 mpfr_nextbelow (z); 1018 mpfr_mul_2ui (z, z, e, MPFR_RNDN); 1019 inex = mpfr_add_ui (z, z, 1, MPFR_RNDN); 1020 MPFR_ASSERTN(inex == 0); 1021 mpfr_div_2ui (z, z, e, MPFR_RNDN); 1022 } 1023 mpfr_clears (x, y, z, (mpfr_ptr) 0); 1024 } 1025