1 /* Test file for mpfr_can_round and mpfr_round_p. 2 3 Copyright 1999, 2001-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 #include "mpfr-test.h" 24 25 /* Simple cases where err - prec is large enough. 26 One can get failures as a consequence of r9883. */ 27 static void 28 test_simple (void) 29 { 30 int t[4] = { 2, 3, -2, -3 }; /* test powers of 2 and non-powers of 2 */ 31 int i; 32 int r1, r2; 33 34 for (i = 0; i < 4; i++) 35 RND_LOOP (r1) 36 RND_LOOP (r2) 37 { 38 mpfr_t b; 39 int p, err, prec, inex, c; 40 41 if (r2 == MPFR_RNDF) 42 continue; 43 p = 12 + (randlimb() % (2 * GMP_NUMB_BITS)); 44 err = p - 3; 45 prec = 4; 46 mpfr_init2 (b, p); 47 inex = mpfr_set_si (b, t[i], MPFR_RNDN); 48 MPFR_ASSERTN (inex == 0); 49 c = mpfr_can_round (b, err, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, prec); 50 /* If r1 == r2, we can round. 51 Note: For r1, RNDF is equivalent to RNDN. 52 TODO: complete this test for r1 != r2. */ 53 if ((r1 == MPFR_RNDF ? MPFR_RNDN : r1) == r2 && !c) 54 { 55 printf ("Error in test_simple for i=%d," 56 " err=%d r1=%s, r2=%s, p=%d\n", i, err, 57 mpfr_print_rnd_mode ((mpfr_rnd_t) r1), 58 mpfr_print_rnd_mode ((mpfr_rnd_t) r2), p); 59 printf ("b="); mpfr_dump (b); 60 exit (1); 61 } 62 mpfr_clear (b); 63 } 64 } 65 66 #define MAX_LIMB_SIZE 100 67 68 static void 69 check_round_p (void) 70 { 71 mp_limb_t buf[MAX_LIMB_SIZE]; 72 mp_size_t n, i; 73 mpfr_prec_t p; 74 mpfr_exp_t err; 75 int r1, r2; 76 77 for (n = 2 ; n <= MAX_LIMB_SIZE ; n++) 78 { 79 /* avoid mpn_random which leaks memory */ 80 for (i = 0; i < n; i++) 81 buf[i] = randlimb (); 82 /* force the number to be normalized */ 83 buf[n - 1] |= MPFR_LIMB_HIGHBIT; 84 p = randlimb() % ((n-1) * GMP_NUMB_BITS) + MPFR_PREC_MIN; 85 err = p + randlimb () % GMP_NUMB_BITS; 86 r1 = mpfr_round_p (buf, n, err, p); 87 r2 = mpfr_can_round_raw (buf, n, MPFR_SIGN_POS, err, 88 MPFR_RNDN, MPFR_RNDZ, p); 89 if (r1 != r2) 90 { 91 printf ("mpfr_round_p(%d) != mpfr_can_round(%d)!\n" 92 "bn = %ld, err0 = %ld, prec = %lu\nbp = ", 93 r1, r2, n, (long) err, (unsigned long) p); 94 #ifndef MPFR_USE_MINI_GMP 95 gmp_printf ("%NX\n", buf, n); 96 #endif 97 exit (1); 98 } 99 } 100 } 101 102 /* check x=2^i with precision px, error at most 1, and target precision prec */ 103 static void 104 test_pow2 (mpfr_exp_t i, mpfr_prec_t px, mpfr_rnd_t r1, mpfr_rnd_t r2, 105 mpfr_prec_t prec) 106 { 107 mpfr_t x; 108 int b, expected_b, b2; 109 110 mpfr_init2 (x, px); 111 mpfr_set_ui_2exp (x, 1, i, MPFR_RNDN); 112 /* for mpfr_can_round, r1=RNDF is equivalent to r1=RNDN (the sign of the 113 error is not known) */ 114 b = !!mpfr_can_round (x, i+1, r1, r2, prec); 115 /* Note: If mpfr_can_round succeeds for both 116 (r1,r2) = (MPFR_RNDD,MPFR_RNDN) and 117 (r1,r2) = (MPFR_RNDU,MPFR_RNDN), then it should succeed for 118 (r1,r2) = (MPFR_RNDN,MPFR_RNDN). So, the condition on prec below 119 for r1 = MPFR_RNDN should be the most restrictive between those 120 for r1 = any directed rounding mode. 121 For r1 like MPFR_RNDA, the unrounded, unknown number may be anyone 122 in [2^i-1,i]. As both 2^i-1 and 2^i fit on i bits, one cannot round 123 in any precision >= i bits, hence the condition prec < i; prec = i-1 124 will work here for r2 = MPFR_RNDN thanks to the even-rounding rule 125 (and also with rounding ties away from zero). */ 126 expected_b = 127 MPFR_IS_LIKE_RNDD (r1, MPFR_SIGN_POS) ? 128 (MPFR_IS_LIKE_RNDU (r2, MPFR_SIGN_POS) ? 0 : prec <= i) : 129 MPFR_IS_LIKE_RNDU (r1, MPFR_SIGN_POS) ? 130 (MPFR_IS_LIKE_RNDD (r2, MPFR_SIGN_POS) ? 0 : prec < i) : 131 (r2 != MPFR_RNDN ? 0 : prec < i); 132 if (b != expected_b) 133 { 134 printf ("Error for x=2^%d, px=%lu, err=%d, r1=%s, r2=%s, prec=%d\n", 135 (int) i, (unsigned long) px, (int) i + 1, 136 mpfr_print_rnd_mode (r1), mpfr_print_rnd_mode (r2), (int) prec); 137 printf ("Expected %d, got %d\n", expected_b, b); 138 exit (1); 139 } 140 141 if (r1 == MPFR_RNDN && r2 == MPFR_RNDZ) 142 { 143 /* Similar test to the one done in src/round_p.c 144 for MPFR_WANT_ASSERT >= 2. */ 145 b2 = !!mpfr_round_p (MPFR_MANT(x), MPFR_LIMB_SIZE(x), i+1, prec); 146 if (b2 != b) 147 { 148 printf ("Error for x=2^%d, px=%lu, err=%d, prec=%d\n", 149 (int) i, (unsigned long) px, (int) i + 1, (int) prec); 150 printf ("mpfr_can_round gave %d, mpfr_round_p gave %d\n", b, b2); 151 exit (1); 152 } 153 } 154 155 mpfr_clear (x); 156 } 157 158 static void 159 check_can_round (void) 160 { 161 mpfr_t x, xinf, xsup, yinf, ysup; 162 int precx, precy, err; 163 int rnd1, rnd2; 164 int i, u[3] = { 0, 1, 256 }; 165 int inex; 166 int expected, got; 167 168 mpfr_inits2 (4 * GMP_NUMB_BITS, x, xinf, xsup, yinf, ysup, (mpfr_ptr) 0); 169 170 for (precx = 3 * GMP_NUMB_BITS - 3; precx <= 3 * GMP_NUMB_BITS + 3; precx++) 171 { 172 mpfr_set_prec (x, precx); 173 for (precy = precx - 4; precy <= precx + 4; precy++) 174 { 175 mpfr_set_prec (yinf, precy); 176 mpfr_set_prec (ysup, precy); 177 178 for (i = 0; i <= 3; i++) 179 { 180 if (i <= 2) 181 { 182 /* Test x = 2^(precx - 1) + u */ 183 mpfr_set_ui_2exp (x, 1, precx - 1, MPFR_RNDN); 184 mpfr_add_ui (x, x, u[i], MPFR_RNDN); 185 } 186 else 187 { 188 /* Test x = 2^precx - 1 */ 189 mpfr_set_ui_2exp (x, 1, precx, MPFR_RNDN); 190 mpfr_sub_ui (x, x, 1, MPFR_RNDN); 191 } 192 MPFR_ASSERTN (mpfr_get_exp (x) == precx); 193 194 for (err = precy; err <= precy + 3; err++) 195 { 196 mpfr_set_ui_2exp (xinf, 1, precx - err, MPFR_RNDN); 197 inex = mpfr_add (xsup, x, xinf, MPFR_RNDN); 198 MPFR_ASSERTN (inex == 0); 199 inex = mpfr_sub (xinf, x, xinf, MPFR_RNDN); 200 MPFR_ASSERTN (inex == 0); 201 RND_LOOP (rnd1) 202 RND_LOOP (rnd2) 203 { 204 if (rnd2 == MPFR_RNDF) 205 continue; 206 mpfr_set (yinf, MPFR_IS_LIKE_RNDD (rnd1, 1) ? 207 x : xinf, (mpfr_rnd_t) rnd2); 208 mpfr_set (ysup, MPFR_IS_LIKE_RNDU (rnd1, 1) ? 209 x : xsup, (mpfr_rnd_t) rnd2); 210 expected = !! mpfr_equal_p (yinf, ysup); 211 got = !! mpfr_can_round (x, err, (mpfr_rnd_t) rnd1, 212 (mpfr_rnd_t) rnd2, precy); 213 if (got != expected) 214 { 215 printf ("Error in check_can_round on:\n" 216 "precx=%d, precy=%d, i=%d, err=%d, " 217 "rnd1=%s, rnd2=%s: got %d\n", 218 precx, precy, i, err, 219 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd1), 220 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd2), 221 got); 222 printf ("x="); mpfr_dump (x); 223 exit (1); 224 } 225 } 226 } 227 } 228 } 229 } 230 231 mpfr_clears (x, xinf, xsup, yinf, ysup, (mpfr_ptr) 0); 232 } 233 234 int 235 main (void) 236 { 237 mpfr_t x; 238 mpfr_prec_t i, j, k; 239 int r1, r2; 240 int n; 241 242 tests_start_mpfr (); 243 244 test_simple (); 245 246 /* checks that rounds to nearest sets the last 247 bit to zero in case of equal distance */ 248 mpfr_init2 (x, 59); 249 mpfr_set_str_binary (x, "-0.10010001010111000011110010111010111110000000111101100111111E663"); 250 if (mpfr_can_round (x, 54, MPFR_RNDZ, MPFR_RNDZ, 53)) 251 { 252 printf ("Error (1) in mpfr_can_round\n"); 253 exit (1); 254 } 255 256 mpfr_set_str_binary (x, "-Inf"); 257 if (mpfr_can_round (x, 2000, MPFR_RNDZ, MPFR_RNDZ, 2000)) 258 { 259 printf ("Error (2) in mpfr_can_round\n"); 260 exit (1); 261 } 262 263 mpfr_set_prec (x, 64); 264 mpfr_set_str_binary (x, "0.1011001000011110000110000110001111101011000010001110011000000000"); 265 if (mpfr_can_round (x, 65, MPFR_RNDN, MPFR_RNDN, 54)) 266 { 267 printf ("Error (3) in mpfr_can_round\n"); 268 exit (1); 269 } 270 271 mpfr_set_prec (x, 137); 272 mpfr_set_str_binary (x, "-0.10111001101001010110011000110100111010011101101010010100101100001110000100111111011101010110001010111100100101110111100001000010000000000E-97"); 273 if (! mpfr_can_round (x, 132, MPFR_RNDU, MPFR_RNDZ, 128)) 274 { 275 printf ("Error (4) in mpfr_can_round\n"); 276 exit (1); 277 } 278 279 /* in the following, we can round but cannot determine the inexact flag */ 280 mpfr_set_prec (x, 86); 281 mpfr_set_str_binary (x, "-0.11100100010011001111011010100111101010011000000000000000000000000000000000000000000000E-80"); 282 if (! mpfr_can_round (x, 81, MPFR_RNDU, MPFR_RNDZ, 44)) 283 { 284 printf ("Error (5) in mpfr_can_round\n"); 285 exit (1); 286 } 287 288 mpfr_set_prec (x, 62); 289 mpfr_set_str (x, "0.ff4ca619c76ba69", 16, MPFR_RNDZ); 290 for (i = 30; i < 99; i++) 291 for (j = 30; j < 99; j++) 292 RND_LOOP (r1) 293 RND_LOOP (r2) 294 if (r2 != MPFR_RNDF) 295 { 296 /* test for assertions */ 297 mpfr_can_round (x, i, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, j); 298 } 299 300 test_pow2 (32, 32, MPFR_RNDN, MPFR_RNDN, 32); 301 test_pow2 (174, 174, MPFR_RNDN, MPFR_RNDN, 174); 302 test_pow2 (174, 174, MPFR_RNDU, MPFR_RNDN, 174); 303 test_pow2 (176, 129, MPFR_RNDU, MPFR_RNDU, 174); 304 test_pow2 (176, 2, MPFR_RNDZ, MPFR_RNDZ, 174); 305 test_pow2 (176, 2, MPFR_RNDU, MPFR_RNDU, 176); 306 307 /* Tests for x = 2^i (E(x) = i+1) with error at most 1 = 2^0. */ 308 for (n = 0; n < 100; n++) 309 { 310 i = (randlimb() % 200) + 4; 311 for (j = i - 2; j < i + 2; j++) 312 RND_LOOP (r1) 313 RND_LOOP (r2) 314 if (r2 != MPFR_RNDF) 315 for (k = MPFR_PREC_MIN; k <= i + 2; k++) 316 test_pow2 (i, k, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, j); 317 } 318 319 mpfr_clear (x); 320 321 check_can_round (); 322 323 check_round_p (); 324 325 tests_end_mpfr (); 326 return 0; 327 } 328