1 /* Test file for mpfr_div_ui.
2
3 Copyright 1999-2023 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 <float.h>
24
25 #include "mpfr-test.h"
26
27 static void
28 check (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
48 static void
49 special (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
149 static void
150 check_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 */
220 static void
221 test_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
252 static void
253 bug20180126 (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 */
364 static void
365 corner_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
416 static void
417 midpoint_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
489 static void
490 check_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
519 int
520 main (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 }