1/* Test file for mpfr_sub1sp. 2 3Copyright 2003-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include "mpfr-test.h" 24 25static void check_special (void); 26static void check_random (mpfr_prec_t p); 27static void check_underflow (mpfr_prec_t p); 28static void check_corner (mpfr_prec_t p); 29 30static void 31bug20170109 (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*/ 53static void 54test20170208 (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 163static void 164compare_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 221static void 222bug20171213 (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*/ 245static void 246bug20171213_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 288static void 289coverage (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 */ 337static void 338bug20180217 (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) */ 375static void 376bug20180813 (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) */ 396static void 397bug20190904 (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 416int 417main (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 472static void 473check_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 RND_LOOP (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 503static void 504check_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 RND_LOOP (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 921static void 922check_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 */ 965static void 966check_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