1 /* Test file for mpfr_can_round and mpfr_round_p.
2
3 Copyright 1999, 2001-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 "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
test_simple(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;
40 int s1, s2;
41 int expected, got;
42
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 got = mpfr_can_round (b, err, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, prec);
50 s1 = r1;
51 s2 = r2;
52 if (s1 == MPFR_RNDD)
53 s1 = (t[i] > 0) ? MPFR_RNDZ : MPFR_RNDA;
54 if (s1 == MPFR_RNDU)
55 s1 = (t[i] < 0) ? MPFR_RNDZ : MPFR_RNDA;
56 if (s1 == MPFR_RNDF)
57 s1 = MPFR_RNDN; /* For s1, RNDF is equivalent to RNDN. */
58 if (s2 == MPFR_RNDD)
59 s2 = (t[i] > 0) ? MPFR_RNDZ : MPFR_RNDA;
60 if (s2 == MPFR_RNDU)
61 s2 = (t[i] < 0) ? MPFR_RNDZ : MPFR_RNDA;
62 /* If s1 == s2, we can round.
63 s1 s2 can round
64 xxx xxx yes
65 RNDZ RNDA no
66 RNDZ RNDN yes
67 RNDA RNDZ no
68 RNDA RNDN yes
69 RNDN RNDZ no
70 RNDN RNDA no
71 xxx RNDF yes
72 */
73 expected = 1;
74 if ((s1 == MPFR_RNDZ && s2 == MPFR_RNDA) ||
75 (s1 == MPFR_RNDA && s2 == MPFR_RNDZ) ||
76 (s1 == MPFR_RNDN && s2 == MPFR_RNDZ) ||
77 (s1 == MPFR_RNDN && s2 == MPFR_RNDA))
78 expected = 0;
79 if (!!got != !!expected)
80 {
81 printf ("Error in test_simple for i=%d,"
82 " err=%d r1=%s, r2=%s, p=%d, prec=%d\n", i, err,
83 mpfr_print_rnd_mode ((mpfr_rnd_t) r1),
84 mpfr_print_rnd_mode ((mpfr_rnd_t) r2), p, prec);
85 printf ("b="); mpfr_dump (b);
86 printf ("expected %d, got %d\n", expected, got);
87 exit (1);
88 }
89 mpfr_clear (b);
90 }
91 }
92
93 #define MAX_LIMB_SIZE 100
94
95 static void
check_round_p(void)96 check_round_p (void)
97 {
98 mp_limb_t buf[MAX_LIMB_SIZE];
99 mp_size_t n, i;
100 mpfr_prec_t p;
101 mpfr_exp_t err;
102 int r1, r2;
103
104 for (n = 2 ; n <= MAX_LIMB_SIZE ; n++)
105 {
106 /* avoid mpn_random which leaks memory */
107 for (i = 0; i < n; i++)
108 buf[i] = randlimb ();
109 /* force the number to be normalized */
110 buf[n - 1] |= MPFR_LIMB_HIGHBIT;
111 p = randlimb() % ((n-1) * GMP_NUMB_BITS) + MPFR_PREC_MIN;
112 err = p + randlimb () % GMP_NUMB_BITS;
113 r1 = mpfr_round_p (buf, n, err, p);
114 r2 = mpfr_can_round_raw (buf, n, MPFR_SIGN_POS, err,
115 MPFR_RNDN, MPFR_RNDZ, p);
116 if (r1 != r2)
117 {
118 printf ("mpfr_round_p(%d) != mpfr_can_round(%d,RNDZ)!\n"
119 "bn = %ld, err0 = %ld, prec = %lu\n",
120 r1, r2, (long) n, (long) err, (unsigned long) p);
121 n_trace ("b", buf, n);
122 exit (1);
123 }
124 /* Same with RNDF: with rnd1=RNDN, rnd2=RNDF is converted to RNDN. */
125 r1 = mpfr_can_round_raw (buf, n, MPFR_SIGN_POS, err,
126 MPFR_RNDN, MPFR_RNDN, p);
127 r2 = mpfr_can_round_raw (buf, n, MPFR_SIGN_POS, err,
128 MPFR_RNDN, MPFR_RNDF, p);
129 if (r1 != r2)
130 {
131 printf ("mpfr_can_round(%d,RNDN) != mpfr_can_round(%d,RNDF)!\n"
132 "bn = %ld, err0 = %ld, prec = %lu\n",
133 r1, r2, (long) n, (long) err, (unsigned long) p);
134 n_trace ("b", buf, n);
135 exit (1);
136 }
137 /* PZ: disabled those tests for now, since when {buf, n} is exactly
138 representable in the target precision p, then mpfr_can_round_raw(RNDA)
139 should give 0, and mpfr_can_round_raw(MPFR_RNDF) should give 1 if the
140 error is small enough. */
141 #if 0
142 /* Same with RNDF: with rnd1=RNDZ, rnd2=RNDF is converted to RNDA. */
143 r1 = mpfr_can_round_raw (buf, n, MPFR_SIGN_POS, err,
144 MPFR_RNDZ, MPFR_RNDA, p);
145 r2 = mpfr_can_round_raw (buf, n, MPFR_SIGN_POS, err,
146 MPFR_RNDZ, MPFR_RNDF, p);
147 if (r1 != r2)
148 {
149 printf ("mpfr_can_round(%d,RNDA) != mpfr_can_round(%d,RNDF)!\n"
150 "bn = %ld, err0 = %ld, prec = %lu\n",
151 r1, r2, (long) n, (long) err, (unsigned long) p);
152 n_trace ("b", buf, n);
153 exit (1);
154 }
155 /* Same with RNDF: with rnd1=RNDA, rnd2=RNDF is converted to RNDZ. */
156 r1 = mpfr_can_round_raw (buf, n, MPFR_SIGN_POS, err,
157 MPFR_RNDA, MPFR_RNDZ, p);
158 r2 = mpfr_can_round_raw (buf, n, MPFR_SIGN_POS, err,
159 MPFR_RNDA, MPFR_RNDF, p);
160 if (r1 != r2)
161 {
162 printf ("mpfr_can_round(%d,RNDZ) != mpfr_can_round(%d,RNDF)!\n"
163 "bn = %ld, err0 = %ld, prec = %lu\n",
164 r1, r2, (long) n, (long) err, (unsigned long) p);
165 n_trace ("b", buf, n);
166 exit (1);
167 }
168 #endif
169 }
170 }
171
172 /* check x=2^i with precision px, error at most 1, and target precision prec */
173 static void
test_pow2(mpfr_exp_t i,mpfr_prec_t px,mpfr_rnd_t r1,mpfr_rnd_t r2,mpfr_prec_t prec)174 test_pow2 (mpfr_exp_t i, mpfr_prec_t px, mpfr_rnd_t r1, mpfr_rnd_t r2,
175 mpfr_prec_t prec)
176 {
177 mpfr_t x;
178 int b, expected_b, b2;
179
180 mpfr_init2 (x, px);
181 mpfr_set_ui_2exp (x, 1, i, MPFR_RNDN);
182 /* for mpfr_can_round, r1=RNDF is equivalent to r1=RNDN (the sign of the
183 error is not known) */
184 b = !!mpfr_can_round (x, i+1, r1, r2, prec);
185 /* Note: If mpfr_can_round succeeds for both
186 (r1,r2) = (MPFR_RNDD,MPFR_RNDN) and
187 (r1,r2) = (MPFR_RNDU,MPFR_RNDN), then it should succeed for
188 (r1,r2) = (MPFR_RNDN,MPFR_RNDN). So, the condition on prec below
189 for r1 = MPFR_RNDN should be the most restrictive between those
190 for r1 = any directed rounding mode.
191 For r1 like MPFR_RNDA, the unrounded, unknown number may be anyone
192 in [2^i-1,i]. As both 2^i-1 and 2^i fit on i bits, one cannot round
193 in any precision >= i bits, hence the condition prec < i; prec = i-1
194 will work here for r2 = MPFR_RNDN thanks to the even-rounding rule
195 (and also with rounding ties away from zero). */
196 expected_b =
197 MPFR_IS_LIKE_RNDD (r1, MPFR_SIGN_POS) ?
198 (MPFR_IS_LIKE_RNDU (r2, MPFR_SIGN_POS) ? 0 : prec <= i) :
199 MPFR_IS_LIKE_RNDU (r1, MPFR_SIGN_POS) ?
200 (MPFR_IS_LIKE_RNDD (r2, MPFR_SIGN_POS) ? 0 : prec < i) :
201 (r2 != MPFR_RNDN ? 0 : prec < i);
202 if (b != expected_b)
203 {
204 printf ("Error for x=2^%d, px=%lu, err=%d, r1=%s, r2=%s, prec=%d\n",
205 (int) i, (unsigned long) px, (int) i + 1,
206 mpfr_print_rnd_mode (r1), mpfr_print_rnd_mode (r2), (int) prec);
207 printf ("Expected %d, got %d\n", expected_b, b);
208 exit (1);
209 }
210
211 if (r1 == MPFR_RNDN && r2 == MPFR_RNDZ)
212 {
213 /* Similar test to the one done in src/round_p.c
214 for MPFR_WANT_ASSERT >= 2. */
215 b2 = !!mpfr_round_p (MPFR_MANT(x), MPFR_LIMB_SIZE(x), i+1, prec);
216 if (b2 != b)
217 {
218 printf ("Error for x=2^%d, px=%lu, err=%d, prec=%d\n",
219 (int) i, (unsigned long) px, (int) i + 1, (int) prec);
220 printf ("mpfr_can_round gave %d, mpfr_round_p gave %d\n", b, b2);
221 exit (1);
222 }
223 }
224
225 mpfr_clear (x);
226 }
227
228 static void
check_can_round(void)229 check_can_round (void)
230 {
231 mpfr_t x, xinf, xsup, yinf, ysup;
232 int precx, precy, err;
233 int rnd1, rnd2;
234 int i, u[3] = { 0, 1, 256 };
235 int inex;
236 int expected, got;
237 int maxerr;
238
239 mpfr_inits2 (4 * GMP_NUMB_BITS, x, xinf, xsup, yinf, ysup, (mpfr_ptr) 0);
240
241 for (precx = 3 * GMP_NUMB_BITS - 3; precx <= 3 * GMP_NUMB_BITS + 3; precx++)
242 {
243 mpfr_set_prec (x, precx);
244 for (precy = precx - 4; precy <= precx + 4; precy++)
245 {
246 mpfr_set_prec (yinf, precy);
247 mpfr_set_prec (ysup, precy);
248
249 for (i = 0; i <= 3; i++)
250 {
251 if (i <= 2)
252 {
253 /* Test x = 2^(precx - 1) + u */
254 mpfr_set_ui_2exp (x, 1, precx - 1, MPFR_RNDN);
255 mpfr_add_ui (x, x, u[i], MPFR_RNDN);
256 }
257 else
258 {
259 /* Test x = 2^precx - 1 */
260 mpfr_set_ui_2exp (x, 1, precx, MPFR_RNDN);
261 mpfr_sub_ui (x, x, 1, MPFR_RNDN);
262 }
263 MPFR_ASSERTN (mpfr_get_exp (x) == precx);
264
265 maxerr = precy + 3;
266 if (4 * GMP_NUMB_BITS < maxerr)
267 maxerr = 4 * GMP_NUMB_BITS;
268 for (err = precy; err <= maxerr; err++)
269 {
270 mpfr_set_ui_2exp (xinf, 1, precx - err, MPFR_RNDN);
271 inex = mpfr_add (xsup, x, xinf, MPFR_RNDN);
272 /* Since EXP(x) = precx, and xinf = 2^(precx-err),
273 x + xinf is exactly representable on 4 * GMP_NUMB_BITS
274 nbits as long as err <= 4 * GMP_NUMB_BITS */
275 MPFR_ASSERTN (inex == 0);
276 inex = mpfr_sub (xinf, x, xinf, MPFR_RNDN);
277 MPFR_ASSERTN (inex == 0);
278 RND_LOOP (rnd1)
279 RND_LOOP (rnd2)
280 {
281 /* TODO: Test r2 == MPFR_RNDF. The following "continue"
282 was added while this case had not been specified
283 yet, but this is no longer the case. */
284 if (rnd2 == MPFR_RNDF)
285 continue;
286 mpfr_set (yinf, MPFR_IS_LIKE_RNDD (rnd1, 1) ?
287 x : xinf, (mpfr_rnd_t) rnd2);
288 mpfr_set (ysup, MPFR_IS_LIKE_RNDU (rnd1, 1) ?
289 x : xsup, (mpfr_rnd_t) rnd2);
290 expected = !! mpfr_equal_p (yinf, ysup);
291 got = !! mpfr_can_round (x, err, (mpfr_rnd_t) rnd1,
292 (mpfr_rnd_t) rnd2, precy);
293 if (got != expected)
294 {
295 printf ("Error in check_can_round on:\n"
296 "precx=%d, precy=%d, i=%d, err=%d, "
297 "rnd1=%s, rnd2=%s: expected %d, got %d\n",
298 precx, precy, i, err,
299 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd1),
300 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd2),
301 expected, got);
302 printf ("x="); mpfr_dump (x);
303 printf ("yinf="); mpfr_dump (yinf);
304 printf ("ysup="); mpfr_dump (ysup);
305 exit (1);
306 }
307 }
308 }
309 }
310 }
311 }
312
313 mpfr_clears (x, xinf, xsup, yinf, ysup, (mpfr_ptr) 0);
314 }
315
316 /* test of RNDNA (nearest with ties to away) */
317 static void
test_rndna(void)318 test_rndna (void)
319 {
320 mpfr_t x;
321 int inex;
322
323 mpfr_init2 (x, 10);
324 mpfr_set_str_binary (x, "1111111101"); /* 1021 */
325 inex = mpfr_prec_round (x, 9, MPFR_RNDNA);
326 MPFR_ASSERTN(inex > 0);
327 MPFR_ASSERTN(mpfr_cmp_ui (x, 1022) == 0);
328 mpfr_set_prec (x, 10);
329 mpfr_set_str_binary (x, "1111111101"); /* 1021 */
330 inex = mpfr_prec_round (x, 9, MPFR_RNDN);
331 MPFR_ASSERTN(inex < 0);
332 MPFR_ASSERTN(mpfr_cmp_ui (x, 1020) == 0);
333 mpfr_set_prec (x, 10);
334 mpfr_set_str_binary (x, "1111111011"); /* 1019 */
335 inex = mpfr_prec_round (x, 9, MPFR_RNDNA);
336 MPFR_ASSERTN(inex > 0);
337 MPFR_ASSERTN(mpfr_cmp_ui (x, 1020) == 0);
338 mpfr_set_prec (x, 10);
339 mpfr_set_str_binary (x, "1111111011"); /* 1019 */
340 inex = mpfr_prec_round (x, 9, MPFR_RNDN);
341 MPFR_ASSERTN(inex > 0);
342 MPFR_ASSERTN(mpfr_cmp_ui (x, 1020) == 0);
343 mpfr_clear (x);
344 }
345
346 int
main(void)347 main (void)
348 {
349 mpfr_t x;
350 mpfr_prec_t i, j, k;
351 int r1, r2;
352 int n;
353
354 tests_start_mpfr ();
355
356 test_rndna ();
357 test_simple ();
358
359 /* checks that rounds to nearest sets the last
360 bit to zero in case of equal distance */
361 mpfr_init2 (x, 59);
362 mpfr_set_str_binary (x, "-0.10010001010111000011110010111010111110000000111101100111111E663");
363 if (mpfr_can_round (x, 54, MPFR_RNDZ, MPFR_RNDZ, 53))
364 {
365 printf ("Error (1) in mpfr_can_round\n");
366 exit (1);
367 }
368
369 mpfr_set_str_binary (x, "-Inf");
370 if (mpfr_can_round (x, 2000, MPFR_RNDZ, MPFR_RNDZ, 2000))
371 {
372 printf ("Error (2) in mpfr_can_round\n");
373 exit (1);
374 }
375
376 mpfr_set_prec (x, 64);
377 mpfr_set_str_binary (x, "0.1011001000011110000110000110001111101011000010001110011000000000");
378 if (mpfr_can_round (x, 65, MPFR_RNDN, MPFR_RNDN, 54))
379 {
380 printf ("Error (3) in mpfr_can_round\n");
381 exit (1);
382 }
383
384 mpfr_set_prec (x, 137);
385 mpfr_set_str_binary (x, "-0.10111001101001010110011000110100111010011101101010010100101100001110000100111111011101010110001010111100100101110111100001000010000000000E-97");
386 if (! mpfr_can_round (x, 132, MPFR_RNDU, MPFR_RNDZ, 128))
387 {
388 printf ("Error (4) in mpfr_can_round\n");
389 exit (1);
390 }
391
392 /* in the following, we can round but cannot determine the inexact flag */
393 mpfr_set_prec (x, 86);
394 mpfr_set_str_binary (x, "-0.11100100010011001111011010100111101010011000000000000000000000000000000000000000000000E-80");
395 if (! mpfr_can_round (x, 81, MPFR_RNDU, MPFR_RNDZ, 44))
396 {
397 printf ("Error (5) in mpfr_can_round\n");
398 exit (1);
399 }
400
401 mpfr_set_prec (x, 62);
402 mpfr_set_str (x, "0.ff4ca619c76ba69", 16, MPFR_RNDZ);
403 for (i = 30; i < 99; i++)
404 for (j = 30; j < 99; j++)
405 RND_LOOP (r1)
406 RND_LOOP (r2)
407 {
408 /* test for assertions */
409 mpfr_can_round (x, i, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, j);
410 }
411
412 test_pow2 (32, 32, MPFR_RNDN, MPFR_RNDN, 32);
413 test_pow2 (174, 174, MPFR_RNDN, MPFR_RNDN, 174);
414 test_pow2 (174, 174, MPFR_RNDU, MPFR_RNDN, 174);
415 test_pow2 (176, 129, MPFR_RNDU, MPFR_RNDU, 174);
416 test_pow2 (176, 2, MPFR_RNDZ, MPFR_RNDZ, 174);
417 test_pow2 (176, 2, MPFR_RNDU, MPFR_RNDU, 176);
418
419 /* Tests for x = 2^i (E(x) = i+1) with error at most 1 = 2^0. */
420 for (n = 0; n < 100; n++)
421 {
422 /* TODO: Test r2 == MPFR_RNDF (add its support in test_pow2). The
423 exclusion below was added while this case had not been specified
424 yet, but this is no longer the case. */
425 i = (randlimb() % 200) + 4;
426 for (j = i - 2; j < i + 2; j++)
427 RND_LOOP (r1)
428 RND_LOOP (r2)
429 if (r2 != MPFR_RNDF)
430 for (k = MPFR_PREC_MIN; k <= i + 2; k++)
431 test_pow2 (i, k, (mpfr_rnd_t) r1, (mpfr_rnd_t) r2, j);
432 }
433
434 mpfr_clear (x);
435
436 check_can_round ();
437
438 check_round_p ();
439
440 tests_end_mpfr ();
441 return 0;
442 }
443