1/* Test file for mpfr_div_ui. 2 3Copyright 1999-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 <float.h> 24 25#include "mpfr-test.h" 26 27static void 28check (const char *ds, unsigned long u, mpfr_rnd_t rnd, const char *es) 29{ 30 mpfr_t x, y; 31 32 mpfr_init2 (x, 53); 33 mpfr_init2 (y, 53); 34 mpfr_set_str1 (x, ds); 35 mpfr_div_ui (y, x, u, rnd); 36 if (mpfr_cmp_str1 (y, es)) 37 { 38 printf ("mpfr_div_ui failed for x=%s, u=%lu, rnd=%s\n", ds, u, 39 mpfr_print_rnd_mode (rnd)); 40 printf ("expected result is %s, got", es); 41 mpfr_out_str(stdout, 10, 0, y, MPFR_RNDN); 42 exit (1); 43 } 44 mpfr_clear (x); 45 mpfr_clear (y); 46} 47 48static void 49special (void) 50{ 51 mpfr_t x, y; 52 unsigned xprec, yprec; 53 54 mpfr_init (x); 55 mpfr_init (y); 56 57 mpfr_set_prec (x, 32); 58 mpfr_set_prec (y, 32); 59 mpfr_set_ui (x, 1, MPFR_RNDN); 60 mpfr_div_ui (y, x, 3, MPFR_RNDN); 61 62 mpfr_set_prec (x, 100); 63 mpfr_set_prec (y, 100); 64 mpfr_urandomb (x, RANDS); 65 mpfr_div_ui (y, x, 123456, MPFR_RNDN); 66 mpfr_set_ui (x, 0, MPFR_RNDN); 67 mpfr_div_ui (y, x, 123456789, MPFR_RNDN); 68 if (mpfr_cmp_ui (y, 0)) 69 { 70 printf ("mpfr_div_ui gives non-zero for 0/ui\n"); 71 exit (1); 72 } 73 74 /* bug found by Norbert Mueller, 21 Aug 2001 */ 75 mpfr_set_prec (x, 110); 76 mpfr_set_prec (y, 60); 77 mpfr_set_str_binary (x, "0.110101110011111110011111001110011001110111000000111110001000111011000011E-44"); 78 mpfr_div_ui (y, x, 17, MPFR_RNDN); 79 mpfr_set_str_binary (x, "0.11001010100101100011101110000001100001010110101001010011011E-48"); 80 if (mpfr_cmp (x, y)) 81 { 82 printf ("Error in x/17 for x=1/16!\n"); 83 printf ("Expected "); 84 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN); 85 printf ("\nGot "); 86 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN); 87 printf ("\n"); 88 exit (1); 89 } 90 91 /* corner case */ 92 mpfr_set_prec (x, 2 * mp_bits_per_limb); 93 mpfr_set_prec (y, 2); 94 mpfr_set_ui (x, 4, MPFR_RNDN); 95 mpfr_nextabove (x); 96 mpfr_div_ui (y, x, 2, MPFR_RNDN); /* exactly in the middle */ 97 MPFR_ASSERTN(mpfr_cmp_ui (y, 2) == 0); 98 (mpfr_div_ui) (y, x, 2, MPFR_RNDN); /* exactly in the middle */ 99 MPFR_ASSERTN(mpfr_cmp_ui (y, 2) == 0); 100 101 mpfr_set_prec (x, 3 * mp_bits_per_limb); 102 mpfr_set_prec (y, 2); 103 mpfr_set_ui (x, 2, MPFR_RNDN); 104 mpfr_nextabove (x); 105 mpfr_div_ui (y, x, 2, MPFR_RNDN); 106 MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0); 107 (mpfr_div_ui) (y, x, 2, MPFR_RNDN); 108 MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0); 109 110 mpfr_set_prec (x, 3 * mp_bits_per_limb); 111 mpfr_set_prec (y, 2); 112 mpfr_set_si (x, -4, MPFR_RNDN); 113 mpfr_nextbelow (x); 114 mpfr_div_ui (y, x, 2, MPFR_RNDD); 115 MPFR_ASSERTN(mpfr_cmp_si (y, -3) == 0); 116 (mpfr_div_ui) (y, x, 2, MPFR_RNDD); 117 MPFR_ASSERTN(mpfr_cmp_si (y, -3) == 0); 118 119 for (xprec = 53; xprec <= 128; xprec++) 120 { 121 mpfr_set_prec (x, xprec); 122 mpfr_set_str_binary (x, "0.1100100100001111110011111000000011011100001100110111E2"); 123 for (yprec = 53; yprec <= 128; yprec++) 124 { 125 mpfr_set_prec (y, yprec); 126 mpfr_div_ui (y, x, 1, MPFR_RNDN); 127 if (mpfr_cmp(x,y)) 128 { 129 printf ("division by 1.0 fails for xprec=%u, yprec=%u\n", xprec, yprec); 130 printf ("expected "); mpfr_dump (x); 131 printf ("got "); mpfr_dump (y); 132 exit (1); 133 } 134 } 135 } 136 137 /* Bug reported by Mark Dickinson, 6 Nov 2007 */ 138 mpfr_set_si (x, 0, MPFR_RNDN); 139 mpfr_set_si (y, -1, MPFR_RNDN); 140 mpfr_div_ui (y, x, 4, MPFR_RNDN); 141 MPFR_ASSERTN(MPFR_IS_ZERO(y) && MPFR_IS_POS(y)); 142 (mpfr_div_ui) (y, x, 4, MPFR_RNDN); 143 MPFR_ASSERTN(MPFR_IS_ZERO(y) && MPFR_IS_POS(y)); 144 145 mpfr_clear (x); 146 mpfr_clear (y); 147} 148 149static void 150check_inexact (void) 151{ 152 mpfr_t x, y, z; 153 mpfr_prec_t px, py; 154 int inexact, cmp; 155 unsigned long int u; 156 int rnd; 157 158 mpfr_init (x); 159 mpfr_init (y); 160 mpfr_init (z); 161 162 for (px=2; px<300; px++) 163 { 164 mpfr_set_prec (x, px); 165 mpfr_urandomb (x, RANDS); 166 do 167 { 168 u = randlimb (); 169 } 170 while (u == 0); 171 for (py=2; py<300; py++) 172 { 173 mpfr_set_prec (y, py); 174 mpfr_set_prec (z, py + mp_bits_per_limb); 175 /* The following test fails with MPFR_RNDF ("Wrong ternary value") 176 when building with CFLAGS="-Wall -Werror -std=c90 -pedantic 177 -Wno-error=overlength-strings -Wno-error=format" so that 178 MPFR_LONG_WITHIN_LIMB is not defined (the implementation 179 is not the same in this case). But the ternary value is not 180 specified for MPFR_RNDF. Thus use RND_LOOP_NO_RNDF. */ 181 RND_LOOP_NO_RNDF (rnd) 182 { 183 inexact = mpfr_div_ui (y, x, u, (mpfr_rnd_t) rnd); 184 if (mpfr_mul_ui (z, y, u, (mpfr_rnd_t) rnd)) 185 { 186 printf ("z <- y * u should be exact for u=%lu\n", u); 187 printf ("y="); mpfr_dump (y); 188 printf ("z="); mpfr_dump (z); 189 exit (1); 190 } 191 cmp = mpfr_cmp (z, x); 192 if (((inexact == 0) && (cmp != 0)) || 193 ((inexact > 0) && (cmp <= 0)) || 194 ((inexact < 0) && (cmp >= 0))) 195 { 196 printf ("Wrong ternary value for u=%lu, rnd=%s\n", u, 197 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); 198 printf ("x="); mpfr_dump (x); 199 printf ("y="); mpfr_dump (y); 200 exit (1); 201 } 202 } 203 } 204 } 205 206 mpfr_clear (x); 207 mpfr_clear (y); 208 mpfr_clear (z); 209} 210 211/* Note that the preprocessor test and the if test in the function are 212 complementary. If the size of an unsigned long is a power of 2 and 213 this type has no padding bits, the test on ULONG_MAX and the if test 214 are equivalent. The preprocessor test can avoid compiler warnings 215 about the large shifts. */ 216#if GMP_NUMB_BITS == 64 && ULONG_MAX > 4294967295 217/* With r11140, on a 64-bit machine with GMP_CHECK_RANDOMIZE=1484406128: 218 Consistency error for i = 2577 219*/ 220static void 221test_20170105 (void) 222{ 223 mpfr_t x,z, t; 224 225 if (sizeof (unsigned long) * CHAR_BIT >= 64) 226 { 227 mpfr_init2 (x, 138); 228 mpfr_init2 (z, 128); 229 mpfr_init2 (t, 128); 230 mpfr_set_str_binary (x, "100110111001001000101111010010011101111110111111110001110100000001110111010100111010100011101010110000010100000011100100110101101011000000"); 231 /* up to exponents, x/y is exactly 232 367625553447399614694201910705139062483, which has 129 bits, 233 thus we are in the round-to-nearest-even case, and since the 234 penultimate bit of x/y is 1, we should round upwards */ 235 mpfr_set_str_binary (t, "10001010010010010000110110010110111111111100011011101010000000000110101000010001011110011011010000111010000000001100101101101010E-53"); 236 mpfr_div_ui (z, x, 36UL << 58, MPFR_RNDN); 237 MPFR_ASSERTN(mpfr_equal_p (z, t)); 238 239 mpfr_set_prec (x, 189); 240 mpfr_set_prec (z, 185); 241 mpfr_set_prec (t, 185); 242 mpfr_set_str_binary (x, "100001010000111100011110111010000011110000000110100010001010101011110001110000110111101000100100001101010011000111110100011111110110011011101000000000001010010010111011001100111111111101001"); 243 mpfr_set_str_binary (t, "10011000000100010100011111100100110101101110001011100101010101011010011010010110010000100111001010000101111011111111001011011010101111101100000000000000101111000100001110101001001001000E-60"); 244 mpfr_div_ui (z, x, 7UL << 61, MPFR_RNDN); 245 MPFR_ASSERTN(mpfr_equal_p (z, t)); 246 247 mpfr_clears (x, z, t, (mpfr_ptr) 0); 248 } 249} 250#endif 251 252static void 253bug20180126 (void) 254{ 255 mpfr_t w, x, y, z, t; 256 unsigned long u, v; 257 int i, k, m, n, p, pmax, q, r; 258 int inex; 259 260 /* Let m = n * q + r, with 0 <= r < v. 261 (2^m-1) / (2^n-1) = 2^r * (2^(n*q)-1) / (2^n-1) + (2^r-1) / (2^n-1) 262 = sum(i=0,q-1,2^(r+n*i)) + sum(i=1,inf,(2^r-1)*2^(-n*i)) 263 */ 264 n = 1; 265 for (u = 1; u != ULONG_MAX; u = (u << 1) + 1) 266 n++; 267 pmax = 6 * n; 268 mpfr_init2 (t, n); 269 for (m = n; m < 4 * n; m++) 270 { 271 q = m / n; 272 r = m % n; 273 mpfr_init2 (w, pmax + n + 1); 274 mpfr_set_zero (w, 1); 275 for (i = 0; i < q; i++) 276 { 277 inex = mpfr_set_ui_2exp (t, 1, r + n * i, MPFR_RNDN); 278 MPFR_ASSERTN (inex == 0); 279 inex = mpfr_add (w, w, t, MPFR_RNDN); 280 MPFR_ASSERTN (inex == 0); 281 } 282 v = (1UL << r) - 1; 283 for (i = 1; n * (q - 1 + i) <= MPFR_PREC (w); i++) 284 { 285 inex = mpfr_set_ui_2exp (t, v, - n * i, MPFR_RNDN); 286 MPFR_ASSERTN (inex == 0); 287 mpfr_add (w, w, t, MPFR_RNDN); 288 } 289 for (p = pmax; p >= MPFR_PREC_MIN; p--) 290 { 291 mpfr_inits2 (p, y, z, (mpfr_ptr) 0); 292 mpfr_set (z, w, MPFR_RNDN); /* the sticky bit is not 0 */ 293 mpfr_init2 (x, m); 294 inex = mpfr_set_ui_2exp (x, 1, m, MPFR_RNDN); 295 MPFR_ASSERTN (inex == 0); 296 inex = mpfr_sub_ui (x, x, 1, MPFR_RNDN); /* x = 2^m-1 */ 297 MPFR_ASSERTN (inex == 0); 298 for (k = 0; k < 2; k++) 299 { 300 if (k) 301 { 302 inex = mpfr_prec_round (x, 6 * n, MPFR_RNDN); 303 MPFR_ASSERTN (inex == 0); 304 } 305 inex = mpfr_div_ui (y, x, u, MPFR_RNDN); 306 if (! mpfr_equal_p (y, z)) 307 { 308 printf ("Error in bug20180126 for (2^%d-1)/(2^%d-1)" 309 " with px=%d py=%d\n", m, n, 310 (int) MPFR_PREC (x), p); 311 printf ("Expected "); 312 mpfr_dump (z); 313 printf ("Got "); 314 mpfr_dump (y); 315 exit (1); 316 } 317 } 318 mpfr_clears (x, y, z, (mpfr_ptr) 0); 319 } 320 mpfr_clear (w); 321 } 322 mpfr_clear (t); 323 324 /* This test expects that a limb fits in an unsigned long. 325 One failing case from function bug20180126() in tdiv.c, 326 for GMP_NUMB_BITS == 64. */ 327#if ULONG_MAX > 4294967295 328 if (GMP_NUMB_BITS == 64 && MPFR_LIMB_MAX <= ULONG_MAX) 329 { 330 mpfr_init2 (x, 133); 331 mpfr_init2 (y, 64); 332 mpfr_set_ui (x, 1, MPFR_RNDN); 333 mpfr_nextbelow (x); /* 1 - 2^(-133) = (2^133-1)/2^133 */ 334 u = MPFR_LIMB_MAX; /* 2^64 - 1 */ 335 inex = mpfr_div_ui (y, x, u, MPFR_RNDN); 336 /* 2^133*x/u = (2^133-1)/(2^64-1) = (2^64+1)*2^5 + 31/(2^64-1) 337 and should be rounded to 2^69+2^6, thus x/u should be rounded 338 to 2^(-133)*(2^69+2^6). */ 339 MPFR_ASSERTN (inex > 0); 340 mpfr_nextbelow (y); 341 MPFR_ASSERTN (mpfr_cmp_ui_2exp (y, 1, -64) == 0); 342 343 mpfr_set_prec (x, 49); 344 mpfr_set_str_binary (x, "0.1000000000000000000111111111111111111111100000000E0"); 345 /* x = 281476050452224/2^49 */ 346 /* let X = 2^256*x = q*u+r, then q has 192 bits, and 347 r = 8222597979955926678 > u/2 thus we should round to (q+1)/2^256 */ 348 mpfr_set_prec (y, 192); 349 /* The cast below avoid spurious warnings from GCC with a 32-bit ABI. */ 350 u = (mp_limb_t) 10865468317030705979U; 351 inex = mpfr_div_ui (y, x, u, MPFR_RNDN); 352 mpfr_init2 (z, 192); 353 mpfr_set_str_binary (z, "0.110110010100111111000100101011011110010101010010001101100110101111001010100011010111010011100001101000110100011101001010000001010000001001011100000100000110101111110100100101011000000110011111E-64"); 354 MPFR_ASSERTN (inex > 0); 355 MPFR_ASSERTN (mpfr_equal_p (y, z)); 356 mpfr_clear (x); 357 mpfr_clear (y); 358 mpfr_clear (z); 359 } 360#endif 361} 362 363/* check corner cases where the round bit is located in the upper bit of r */ 364static void 365corner_cases (int n) 366{ 367 mpfr_t x, y, t; 368 unsigned long u, v; 369 int i, xn; 370 371 if (MPFR_LIMB_MAX <= ULONG_MAX) 372 { 373 /* We need xn > yn + 1, thus we take xn = 3 and yn = 1. 374 Also take xn = 4 to 6 to cover more code. */ 375 for (xn = 3; xn < 6; xn++) 376 { 377 mpfr_init2 (x, xn * GMP_NUMB_BITS); 378 mpfr_init2 (y, GMP_NUMB_BITS); 379 mpfr_init2 (t, 2 * GMP_NUMB_BITS); 380 for (i = 0; i < n; i++) 381 { 382 do 383 u = randlimb (); 384 while (u == 0); 385 do 386 v = randlimb (); 387 while (v <= MPFR_LIMB_HIGHBIT); 388 mpfr_set_ui (t, v, MPFR_RNDN); 389 mpfr_sub_d (t, t, 0.5, MPFR_RNDN); 390 /* t = v-1/2 */ 391 mpfr_mul_ui (x, t, u, MPFR_RNDN); 392 393 /* when x = (v-1/2)*u, x/u should give v-1/2, which should 394 round to either v (if v is even) or v-1 (if v is odd) */ 395 mpfr_div_ui (y, x, u, MPFR_RNDN); 396 MPFR_ASSERTN(mpfr_cmp_ui (y, v - (v & 1)) == 0); 397 398 /* when x = (v-1/2)*u - epsilon, x/u should round to v-1 */ 399 mpfr_nextbelow (x); 400 mpfr_div_ui (y, x, u, MPFR_RNDN); 401 MPFR_ASSERTN(mpfr_cmp_ui (y, v - 1) == 0); 402 403 /* when x = (v-1/2)*u + epsilon, x/u should round to v */ 404 mpfr_nextabove (x); 405 mpfr_nextabove (x); 406 mpfr_div_ui (y, x, u, MPFR_RNDN); 407 MPFR_ASSERTN(mpfr_cmp_ui (y, v) == 0); 408 } 409 mpfr_clear (x); 410 mpfr_clear (y); 411 mpfr_clear (t); 412 } 413 } 414} 415 416static void 417midpoint_exact (void) 418{ 419 mpfr_t x, y1, y2; 420 unsigned long j; 421 int i, kx, ky, px, pxmin, py, pymin, r; 422 int inex1, inex2; 423 424 pymin = 1; 425 for (i = 3; i < 32; i += 2) 426 { 427 if ((i & (i-2)) == 1) 428 pymin++; 429 for (j = 1; j != 0; j++) 430 { 431 if (j == 31) 432 j = ULONG_MAX; 433 /* Test of (i*j) / j with various precisions. The target precisions 434 include: large, length(i), and length(i)-1; the latter case 435 corresponds to a midpoint. */ 436 mpfr_init2 (x, 5 + sizeof(long) * CHAR_BIT); 437 inex1 = mpfr_set_ui (x, j, MPFR_RNDN); 438 MPFR_ASSERTN (inex1 == 0); 439 inex1 = mpfr_mul_ui (x, x, i, MPFR_RNDN); 440 MPFR_ASSERTN (inex1 == 0); 441 /* x = (i*j) */ 442 pxmin = mpfr_min_prec (x); 443 if (pxmin < MPFR_PREC_MIN) 444 pxmin = MPFR_PREC_MIN; 445 for (kx = 0; kx < 8; kx++) 446 { 447 px = pxmin; 448 if (kx != 0) 449 px += randlimb () % (4 * GMP_NUMB_BITS); 450 inex1 = mpfr_prec_round (x, px, MPFR_RNDN); 451 MPFR_ASSERTN (inex1 == 0); 452 for (ky = 0; ky < 8; ky++) 453 { 454 py = pymin; 455 if (ky == 0) 456 py--; 457 else if (ky > 1) 458 py += randlimb () % (4 * GMP_NUMB_BITS); 459 if (py < MPFR_PREC_MIN) 460 break; 461 mpfr_inits2 (py, y1, y2, (mpfr_ptr) 0); 462 RND_LOOP_NO_RNDF (r) 463 { 464 inex1 = mpfr_set_ui (y1, i, (mpfr_rnd_t) r); 465 inex2 = mpfr_div_ui (y2, x, j, (mpfr_rnd_t) r); 466 if (! mpfr_equal_p (y1, y2) || 467 ! SAME_SIGN (inex1, inex2)) 468 { 469 printf ("Error in midpoint_exact for " 470 "i=%d j=%lu px=%d py=%d %s\n", i, j, px, py, 471 mpfr_print_rnd_mode ((mpfr_rnd_t) r)); 472 printf ("Expected "); 473 mpfr_dump (y1); 474 printf ("with inex = %d\n", inex1); 475 printf ("Got "); 476 mpfr_dump (y2); 477 printf ("with inex = %d\n", inex2); 478 exit (1); 479 } 480 } 481 mpfr_clears (y1, y2, (mpfr_ptr) 0); 482 } 483 } 484 mpfr_clear (x); 485 } 486 } 487} 488 489static void 490check_coverage (void) 491{ 492#ifdef MPFR_COV_CHECK 493 int i, j; 494 int err = 0; 495 496 if (MPFR_LIMB_MAX <= ULONG_MAX) 497 { 498 for (i = 0; i < numberof (__gmpfr_cov_div_ui_sb); i++) 499 for (j = 0; j < 2; j++) 500 if (!__gmpfr_cov_div_ui_sb[i][j]) 501 { 502 printf ("mpfr_div_ui not tested on case %d, sb=%d\n", i, j); 503 err = 1; 504 } 505 506 if (err) 507 exit (1); 508 } 509 else /* e.g. mips64 with the n32 ABI */ 510 printf ("Warning! Value coverage disabled (mp_limb_t > unsigned long).\n"); 511#endif 512} 513 514#define TEST_FUNCTION mpfr_div_ui 515#define ULONG_ARG2 516#define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS) 517#include "tgeneric.c" 518 519int 520main (int argc, char **argv) 521{ 522 mpfr_t x; 523 524 tests_start_mpfr (); 525 526 corner_cases (100); 527 528 bug20180126 (); 529 530 special (); 531 532 check_inexact (); 533 534 check("1.0", 3, MPFR_RNDN, "3.3333333333333331483e-1"); 535 check("1.0", 3, MPFR_RNDZ, "3.3333333333333331483e-1"); 536 check("1.0", 3, MPFR_RNDU, "3.3333333333333337034e-1"); 537 check("1.0", 3, MPFR_RNDD, "3.3333333333333331483e-1"); 538 check("1.0", 2116118, MPFR_RNDN, "4.7256343927890600483e-7"); 539 check("1.098612288668109782", 5, MPFR_RNDN, "0.21972245773362195087"); 540#if GMP_NUMB_BITS == 64 && ULONG_MAX > 4294967295 541 test_20170105 (); 542#endif 543 544 mpfr_init2 (x, 53); 545 mpfr_set_ui (x, 3, MPFR_RNDD); 546 mpfr_log (x, x, MPFR_RNDD); 547 mpfr_div_ui (x, x, 5, MPFR_RNDD); 548 if (mpfr_cmp_str1 (x, "0.21972245773362189536")) 549 { 550 printf ("Error in mpfr_div_ui for x=ln(3), u=5\n"); 551 exit (1); 552 } 553 mpfr_clear (x); 554 555 test_generic (MPFR_PREC_MIN, 200, 100); 556 midpoint_exact (); 557 558 check_coverage (); 559 tests_end_mpfr (); 560 return 0; 561} 562