1 /* Test file for mpfr_sub.
2
3 Copyright 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 static int
test_sub(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)26 test_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
27 {
28 #ifdef CHECK_EXTERNAL
29 int res;
30 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c);
31
32 if (ok)
33 {
34 mpfr_print_raw (b);
35 printf (" ");
36 mpfr_print_raw (c);
37 }
38 res = mpfr_sub (a, b, c, rnd_mode);
39 if (ok)
40 {
41 printf (" ");
42 mpfr_print_raw (a);
43 printf ("\n");
44 }
45 return res;
46 #else /* reuse test */
47 int inex;
48
49 inex = mpfr_sub (a, b, c, rnd_mode);
50
51 if (a != b && a != c && ! MPFR_IS_NAN (a))
52 {
53 mpfr_t t;
54 int reuse_b, reuse_c, inex_r;
55
56 reuse_b = MPFR_PREC (a) == MPFR_PREC (b);
57 reuse_c = MPFR_PREC (a) == MPFR_PREC (c);
58
59 if (reuse_b || reuse_c)
60 mpfr_init2 (t, MPFR_PREC (a));
61
62 if (reuse_b)
63 {
64 mpfr_set (t, b, MPFR_RNDN);
65 inex_r = mpfr_sub (t, t, c, rnd_mode);
66 if (!(mpfr_equal_p (t, a) && SAME_SIGN (inex_r, inex)))
67 {
68 printf ("reuse of b error in b - c in %s for\n",
69 mpfr_print_rnd_mode (rnd_mode));
70 printf ("b = ");
71 mpfr_dump (b);
72 printf ("c = ");
73 mpfr_dump (c);
74 printf ("Expected "); mpfr_dump (a);
75 printf (" with inex = %d\n", inex);
76 printf ("Got "); mpfr_dump (t);
77 printf (" with inex = %d\n", inex_r);
78 exit (1);
79 }
80 }
81
82 if (reuse_c)
83 {
84 mpfr_set (t, c, MPFR_RNDN);
85 inex_r = mpfr_sub (t, b, t, rnd_mode);
86 if (!(mpfr_equal_p (t, a) && SAME_SIGN (inex_r, inex)))
87 {
88 printf ("reuse of c error in b - c in %s for\n",
89 mpfr_print_rnd_mode (rnd_mode));
90 printf ("b = ");
91 mpfr_dump (b);
92 printf ("c = ");
93 mpfr_dump (c);
94 printf ("Expected "); mpfr_dump (a);
95 printf (" with inex = %d\n", inex);
96 printf ("Got "); mpfr_dump (t);
97 printf (" with inex = %d\n", inex_r);
98 exit (1);
99 }
100 }
101
102 if (reuse_b || reuse_c)
103 mpfr_clear (t);
104 }
105
106 return inex;
107 #endif
108 }
109
110 static void
check_diverse(void)111 check_diverse (void)
112 {
113 mpfr_t x, y, z;
114 int inexact;
115
116 mpfr_init (x);
117 mpfr_init (y);
118 mpfr_init (z);
119
120 /* check corner case cancel=0, but add_exp=1 */
121 mpfr_set_prec (x, 2);
122 mpfr_set_prec (y, 4);
123 mpfr_set_prec (z, 2);
124 mpfr_setmax (y, __gmpfr_emax);
125 mpfr_set_str_binary (z, "0.1E-10"); /* tiny */
126 test_sub (x, y, z, MPFR_RNDN); /* should round to 2^emax, i.e. overflow */
127 if (!mpfr_inf_p (x) || mpfr_sgn (x) < 0)
128 {
129 printf ("Error in mpfr_sub(a,b,c,RNDN) for b=maxfloat, prec(a)<prec(b), c tiny\n");
130 exit (1);
131 }
132
133 /* other coverage test */
134 mpfr_set_prec (x, 2);
135 mpfr_set_prec (y, 2);
136 mpfr_set_prec (z, 2);
137 mpfr_set_ui (y, 1, MPFR_RNDN);
138 mpfr_set_si (z, -2, MPFR_RNDN);
139 test_sub (x, y, z, MPFR_RNDD);
140 if (mpfr_cmp_ui (x, 3))
141 {
142 printf ("Error in mpfr_sub(1,-2,RNDD)\n");
143 exit (1);
144 }
145
146 /* yet another coverage test */
147 mpfr_set_prec (x, 2);
148 mpfr_set_prec (y, 3);
149 mpfr_set_prec (z, 1);
150 mpfr_set_ui_2exp (y, 1, mpfr_get_emax (), MPFR_RNDZ);
151 /* y = (1 - 2^(-3))*2^emax */
152 mpfr_set_ui_2exp (z, 1, mpfr_get_emax () - 4, MPFR_RNDZ);
153 /* z = 2^(emax - 4) */
154 /* y - z = (1 - 2^(-3) - 2^(-4))*2^emax > (1-2^(-2))*2^emax */
155 inexact = mpfr_sub (x, y, z, MPFR_RNDU);
156 MPFR_ASSERTN(inexact > 0);
157 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);
158
159 mpfr_set_prec (x, 288);
160 mpfr_set_prec (y, 288);
161 mpfr_set_prec (z, 288);
162 mpfr_set_str_binary (y, "0.111000110011000001000111101010111011110011101001101111111110000011100101000001001010110010101010011001010100000001110011110001010101101010001011101110100100001011110100110000101101100011010001001011011010101010000010001101001000110010010111111011110001111101001000101101001100101100101000E80");
163 mpfr_set_str_binary (z, "0.100001111111101001011010001100110010100111001110000110011101001011010100001000000100111011010110110010000000000010101101011000010000110001110010100001100101011100100100001011000100011110000001010101000100011101001000010111100000111000111011001000100100011000100000010010111000000100100111E-258");
164 inexact = test_sub (x, y, z, MPFR_RNDN);
165 if (inexact <= 0)
166 {
167 printf ("Wrong inexact flag for prec=288\n");
168 exit (1);
169 }
170
171 mpfr_set_prec (x, 32);
172 mpfr_set_prec (y, 63);
173 mpfr_set_prec (z, 63);
174 mpfr_set_str_binary (x, "0.101101111011011100100100100111E31");
175 mpfr_set_str_binary (y, "0.111110010010100100110101101010001001100101110001000101110111111E-1");
176 test_sub (z, x, y, MPFR_RNDN);
177 mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
178 if (mpfr_cmp (z, y))
179 {
180 printf ("Error in mpfr_sub (5)\n");
181 printf ("expected "); mpfr_dump (y);
182 printf ("got "); mpfr_dump (z);
183 exit (1);
184 }
185
186 mpfr_set_prec (y, 63);
187 mpfr_set_prec (z, 63);
188 mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
189 mpfr_sub_ui (z, y, 1541116494, MPFR_RNDN);
190 mpfr_set_str_binary (y, "-0.11111001001010010011010110101E-1");
191 if (mpfr_cmp (z, y))
192 {
193 printf ("Error in mpfr_sub (7)\n");
194 printf ("expected "); mpfr_dump (y);
195 printf ("got "); mpfr_dump (z);
196 exit (1);
197 }
198
199 mpfr_set_prec (y, 63);
200 mpfr_set_prec (z, 63);
201 mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
202 mpfr_sub_ui (z, y, 1541116494, MPFR_RNDN);
203 mpfr_set_str_binary (y, "-0.11111001001010010011010110101E-1");
204 if (mpfr_cmp (z, y))
205 {
206 printf ("Error in mpfr_sub (6)\n");
207 printf ("expected "); mpfr_dump (y);
208 printf ("got "); mpfr_dump (z);
209 exit (1);
210 }
211
212 mpfr_set_prec (x, 32);
213 mpfr_set_prec (y, 32);
214 mpfr_set_str_binary (x, "0.10110111101001110100100101111000E0");
215 mpfr_set_str_binary (y, "0.10001100100101000100110111000100E0");
216 if ((inexact = test_sub (x, x, y, MPFR_RNDN)))
217 {
218 printf ("Wrong inexact flag (2): got %d instead of 0\n", inexact);
219 exit (1);
220 }
221
222 mpfr_set_prec (x, 32);
223 mpfr_set_prec (y, 32);
224 mpfr_set_str_binary (x, "0.11111000110111011000100111011010E0");
225 mpfr_set_str_binary (y, "0.10011111101111000100001000000000E-3");
226 if ((inexact = test_sub (x, x, y, MPFR_RNDN)))
227 {
228 printf ("Wrong inexact flag (1): got %d instead of 0\n", inexact);
229 exit (1);
230 }
231
232 mpfr_set_prec (x, 33);
233 mpfr_set_prec (y, 33);
234 mpfr_set_ui (x, 1, MPFR_RNDN);
235 mpfr_set_str_binary (y, "-0.1E-32");
236 mpfr_add (x, x, y, MPFR_RNDN);
237 mpfr_set_str_binary (y, "0.111111111111111111111111111111111E0");
238 if (mpfr_cmp (x, y))
239 {
240 printf ("Error in mpfr_sub (1 - 1E-33) with prec=33\n");
241 printf ("Expected "); mpfr_dump (y);
242 printf ("got "); mpfr_dump (x);
243 exit (1);
244 }
245
246 mpfr_set_prec (x, 32);
247 mpfr_set_prec (y, 33);
248 mpfr_set_ui (x, 1, MPFR_RNDN);
249 mpfr_set_str_binary (y, "-0.1E-32");
250 mpfr_add (x, x, y, MPFR_RNDN);
251 if (mpfr_cmp_ui (x, 1))
252 {
253 printf ("Error in mpfr_sub (1 - 1E-33) with prec=32\n");
254 printf ("Expected 1.0, got "); mpfr_dump (x);
255 exit (1);
256 }
257
258 mpfr_set_prec (x, 65);
259 mpfr_set_prec (y, 65);
260 mpfr_set_prec (z, 64);
261 mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
262 mpfr_set_str_binary (y, "0.1110111011110001110111011111111111101000011001011100101100101100");
263 test_sub (z, x, y, MPFR_RNDZ);
264 if (mpfr_cmp_ui (z, 1))
265 {
266 printf ("Error in mpfr_sub (1)\n");
267 exit (1);
268 }
269 test_sub (z, x, y, MPFR_RNDU);
270 mpfr_set_str_binary (x, "1.000000000000000000000000000000000000000000000000000000000000001");
271 if (mpfr_cmp (z, x))
272 {
273 printf ("Error in mpfr_sub (2)\n");
274 printf ("Expected "); mpfr_dump (x);
275 printf ("Got "); mpfr_dump (z);
276 exit (1);
277 }
278 mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
279 test_sub (z, x, y, MPFR_RNDN);
280 if (mpfr_cmp_ui (z, 1))
281 {
282 printf ("Error in mpfr_sub (3)\n");
283 exit (1);
284 }
285 inexact = test_sub (z, x, y, MPFR_RNDA);
286 mpfr_set_str_binary (x, "1.000000000000000000000000000000000000000000000000000000000000001");
287 if (mpfr_cmp (z, x) || inexact <= 0)
288 {
289 printf ("Error in mpfr_sub (4)\n");
290 exit (1);
291 }
292 mpfr_set_prec (x, 66);
293 mpfr_set_str_binary (x, "1.11101110111100011101110111111111111010000110010111001011001010111");
294 test_sub (z, x, y, MPFR_RNDN);
295 if (mpfr_cmp_ui (z, 1))
296 {
297 printf ("Error in mpfr_sub (5)\n");
298 exit (1);
299 }
300
301 /* check in-place operations */
302 mpfr_set_ui (x, 1, MPFR_RNDN);
303 test_sub (x, x, x, MPFR_RNDN);
304 if (mpfr_cmp_ui(x, 0))
305 {
306 printf ("Error for mpfr_sub (x, x, x, MPFR_RNDN) with x=1.0\n");
307 exit (1);
308 }
309
310 mpfr_set_prec (x, 53);
311 mpfr_set_prec (y, 53);
312 mpfr_set_prec (z, 53);
313 mpfr_set_str1 (x, "1.229318102e+09");
314 mpfr_set_str1 (y, "2.32221184180698677665e+05");
315 test_sub (z, x, y, MPFR_RNDN);
316 if (mpfr_cmp_str1 (z, "1229085880.815819263458251953125"))
317 {
318 printf ("Error in mpfr_sub (1.22e9 - 2.32e5)\n");
319 printf ("expected 1229085880.815819263458251953125, got ");
320 mpfr_out_str(stdout, 10, 0, z, MPFR_RNDN);
321 putchar('\n');
322 exit (1);
323 }
324
325 mpfr_set_prec (x, 112);
326 mpfr_set_prec (y, 98);
327 mpfr_set_prec (z, 54);
328 mpfr_set_str_binary (x, "0.11111100100000000011000011100000101101010001000111E-401");
329 mpfr_set_str_binary (y, "0.10110000100100000101101100011111111011101000111000101E-464");
330 test_sub (z, x, y, MPFR_RNDN);
331 if (mpfr_cmp (z, x)) {
332 printf ("mpfr_sub(z, x, y) failed for prec(x)=112, prec(y)=98\n");
333 printf ("expected "); mpfr_dump (x);
334 printf ("got "); mpfr_dump (z);
335 exit (1);
336 }
337
338 mpfr_set_prec (x, 33);
339 mpfr_set_ui (x, 1, MPFR_RNDN);
340 mpfr_div_2ui (x, x, 32, MPFR_RNDN);
341 mpfr_sub_ui (x, x, 1, MPFR_RNDN);
342
343 mpfr_set_prec (x, 5);
344 mpfr_set_prec (y, 5);
345 mpfr_set_str_binary (x, "1e-12");
346 mpfr_set_ui (y, 1, MPFR_RNDN);
347 test_sub (x, y, x, MPFR_RNDD);
348 mpfr_set_str_binary (y, "0.11111");
349 if (mpfr_cmp (x, y))
350 {
351 printf ("Error in mpfr_sub (x, y, x, MPFR_RNDD) for x=2^(-12), y=1\n");
352 exit (1);
353 }
354
355 mpfr_set_prec (x, 24);
356 mpfr_set_prec (y, 24);
357 mpfr_set_str_binary (x, "-0.100010000000000000000000E19");
358 mpfr_set_str_binary (y, "0.100000000000000000000100E15");
359 mpfr_add (x, x, y, MPFR_RNDD);
360 mpfr_set_str_binary (y, "-0.1E19");
361 if (mpfr_cmp (x, y))
362 {
363 printf ("Error in mpfr_add (2)\n");
364 exit (1);
365 }
366
367 mpfr_set_prec (x, 2);
368 mpfr_set_prec (y, 10);
369 mpfr_set_prec (z, 10);
370 mpfr_set_ui (y, 0, MPFR_RNDN);
371 mpfr_set_str_binary (z, "0.10001");
372 if (test_sub (x, y, z, MPFR_RNDN) <= 0)
373 {
374 printf ("Wrong inexact flag in x=mpfr_sub(0,z) for prec(z)>prec(x)\n");
375 exit (1);
376 }
377 if (test_sub (x, z, y, MPFR_RNDN) >= 0)
378 {
379 printf ("Wrong inexact flag in x=mpfr_sub(z,0) for prec(z)>prec(x)\n");
380 exit (1);
381 }
382
383 mpfr_clear (x);
384 mpfr_clear (y);
385 mpfr_clear (z);
386 }
387
388 static void
bug_ddefour(void)389 bug_ddefour(void)
390 {
391 mpfr_t ex, ex1, ex2, ex3, tot, tot1;
392
393 mpfr_init2(ex, 53);
394 mpfr_init2(ex1, 53);
395 mpfr_init2(ex2, 53);
396 mpfr_init2(ex3, 53);
397 mpfr_init2(tot, 150);
398 mpfr_init2(tot1, 150);
399
400 mpfr_set_ui( ex, 1, MPFR_RNDN);
401 mpfr_mul_2ui( ex, ex, 906, MPFR_RNDN);
402 mpfr_log( tot, ex, MPFR_RNDN);
403 mpfr_set( ex1, tot, MPFR_RNDN); /* ex1 = high(tot) */
404 test_sub( ex2, tot, ex1, MPFR_RNDN); /* ex2 = high(tot - ex1) */
405 test_sub( tot1, tot, ex1, MPFR_RNDN); /* tot1 = tot - ex1 */
406 mpfr_set( ex3, tot1, MPFR_RNDN); /* ex3 = high(tot - ex1) */
407
408 if (mpfr_cmp(ex2, ex3))
409 {
410 printf ("Error in ddefour test.\n");
411 printf ("ex2="); mpfr_dump (ex2);
412 printf ("ex3="); mpfr_dump (ex3);
413 exit (1);
414 }
415
416 mpfr_clear (ex);
417 mpfr_clear (ex1);
418 mpfr_clear (ex2);
419 mpfr_clear (ex3);
420 mpfr_clear (tot);
421 mpfr_clear (tot1);
422 }
423
424 /* if u = o(x-y), v = o(u-x), w = o(v+y), then x-y = u-w */
425 static void
check_two_sum(mpfr_prec_t p)426 check_two_sum (mpfr_prec_t p)
427 {
428 mpfr_t x, y, u, v, w;
429 mpfr_rnd_t rnd;
430 int inexact;
431
432 mpfr_init2 (x, p);
433 mpfr_init2 (y, p);
434 mpfr_init2 (u, p);
435 mpfr_init2 (v, p);
436 mpfr_init2 (w, p);
437 mpfr_urandomb (x, RANDS);
438 mpfr_urandomb (y, RANDS);
439 if (mpfr_cmpabs (x, y) < 0)
440 mpfr_swap (x, y);
441 rnd = MPFR_RNDN;
442 inexact = test_sub (u, x, y, rnd);
443 test_sub (v, u, x, rnd);
444 mpfr_add (w, v, y, rnd);
445 /* as u = (x-y) - w, we should have inexact and w of opposite signs */
446 if (((inexact == 0) && mpfr_cmp_ui (w, 0)) ||
447 ((inexact > 0) && (mpfr_cmp_ui (w, 0) <= 0)) ||
448 ((inexact < 0) && (mpfr_cmp_ui (w, 0) >= 0)))
449 {
450 printf ("Wrong inexact flag for prec=%u, rnd=%s\n", (unsigned)p,
451 mpfr_print_rnd_mode (rnd));
452 printf ("x="); mpfr_dump (x);
453 printf ("y="); mpfr_dump (y);
454 printf ("u="); mpfr_dump (u);
455 printf ("v="); mpfr_dump (v);
456 printf ("w="); mpfr_dump (w);
457 printf ("inexact = %d\n", inexact);
458 exit (1);
459 }
460 mpfr_clear (x);
461 mpfr_clear (y);
462 mpfr_clear (u);
463 mpfr_clear (v);
464 mpfr_clear (w);
465 }
466
467 #define MAX_PREC 200
468
469 static void
check_inexact(void)470 check_inexact (void)
471 {
472 mpfr_t x, y, z, u;
473 mpfr_prec_t px, py, pu, pz;
474 int inexact, cmp;
475 mpfr_rnd_t rnd;
476
477 mpfr_init (x);
478 mpfr_init (y);
479 mpfr_init (z);
480 mpfr_init (u);
481
482 mpfr_set_prec (x, 2);
483 mpfr_set_ui (x, 6, MPFR_RNDN);
484 mpfr_div_2ui (x, x, 4, MPFR_RNDN); /* x = 6/16 */
485 mpfr_set_prec (y, 2);
486 mpfr_set_si (y, -1, MPFR_RNDN);
487 mpfr_div_2ui (y, y, 4, MPFR_RNDN); /* y = -1/16 */
488 inexact = test_sub (y, y, x, MPFR_RNDN); /* y = round(-7/16) = -1/2 */
489 if (inexact >= 0)
490 {
491 printf ("Error: wrong inexact flag for -1/16 - (6/16)\n");
492 exit (1);
493 }
494
495 for (px=2; px<MAX_PREC; px++)
496 {
497 mpfr_set_prec (x, px);
498 do
499 {
500 mpfr_urandomb (x, RANDS);
501 }
502 while (mpfr_cmp_ui (x, 0) == 0);
503 for (pu=2; pu<MAX_PREC; pu++)
504 {
505 mpfr_set_prec (u, pu);
506 do
507 {
508 mpfr_urandomb (u, RANDS);
509 }
510 while (mpfr_cmp_ui (u, 0) == 0);
511 {
512 py = 2 + (randlimb () % (MAX_PREC - 2));
513 mpfr_set_prec (y, py);
514 /* warning: MPFR_EXP is undefined for 0 */
515 pz = (mpfr_cmpabs (x, u) >= 0) ? MPFR_EXP(x) - MPFR_EXP(u)
516 : MPFR_EXP(u) - MPFR_EXP(x);
517 pz = pz + MAX(MPFR_PREC(x), MPFR_PREC(u));
518 mpfr_set_prec (z, pz);
519 rnd = RND_RAND_NO_RNDF ();
520 if (test_sub (z, x, u, rnd))
521 {
522 printf ("z <- x - u should be exact\n");
523 exit (1);
524 }
525 {
526 rnd = RND_RAND_NO_RNDF ();
527 inexact = test_sub (y, x, u, rnd);
528 cmp = mpfr_cmp (y, z);
529 if (((inexact == 0) && (cmp != 0)) ||
530 ((inexact > 0) && (cmp <= 0)) ||
531 ((inexact < 0) && (cmp >= 0)))
532 {
533 printf ("Wrong inexact flag for rnd=%s\n",
534 mpfr_print_rnd_mode(rnd));
535 printf ("expected %d, got %d\n", cmp, inexact);
536 printf ("x="); mpfr_dump (x);
537 printf ("u="); mpfr_dump (u);
538 printf ("y= "); mpfr_dump (y);
539 printf ("x-u="); mpfr_dump (z);
540 exit (1);
541 }
542 }
543 }
544 }
545 }
546
547 mpfr_clear (x);
548 mpfr_clear (y);
549 mpfr_clear (z);
550 mpfr_clear (u);
551 }
552
553 /* Bug found by Jakub Jelinek
554 * https://bugzilla.redhat.com/show_bug.cgi?id=643657
555 * https://gforge.inria.fr/tracker/index.php?func=detail&aid=11301
556 * The consequence can be either an assertion failure (i = 2 in the
557 * testcase below, in debug mode) or an incorrectly rounded value.
558 */
559 static void
bug20101017(void)560 bug20101017 (void)
561 {
562 mpfr_t a, b, c;
563 int inex;
564 int i;
565
566 mpfr_init2 (a, GMP_NUMB_BITS * 2);
567 mpfr_init2 (b, GMP_NUMB_BITS);
568 mpfr_init2 (c, GMP_NUMB_BITS);
569
570 /* a = 2^(2N) + k.2^(2N-1) + 2^N and b = 1
571 with N = GMP_NUMB_BITS and k = 0 or 1.
572 c = a - b should round to the same value as a. */
573
574 for (i = 2; i <= 3; i++)
575 {
576 mpfr_set_ui_2exp (a, i, GMP_NUMB_BITS - 1, MPFR_RNDN);
577 mpfr_add_ui (a, a, 1, MPFR_RNDN);
578 mpfr_mul_2ui (a, a, GMP_NUMB_BITS, MPFR_RNDN);
579 mpfr_set_ui (b, 1, MPFR_RNDN);
580 inex = mpfr_sub (c, a, b, MPFR_RNDN);
581 mpfr_set (b, a, MPFR_RNDN);
582 if (! mpfr_equal_p (c, b))
583 {
584 printf ("Error in bug20101017 for i = %d.\n", i);
585 printf ("Expected ");
586 mpfr_out_str (stdout, 16, 0, b, MPFR_RNDN);
587 putchar ('\n');
588 printf ("Got ");
589 mpfr_out_str (stdout, 16, 0, c, MPFR_RNDN);
590 putchar ('\n');
591 exit (1);
592 }
593 if (inex >= 0)
594 {
595 printf ("Error in bug20101017 for i = %d: bad inex value.\n", i);
596 printf ("Expected negative, got %d.\n", inex);
597 exit (1);
598 }
599 }
600
601 mpfr_set_prec (a, 64);
602 mpfr_set_prec (b, 129);
603 mpfr_set_prec (c, 2);
604 mpfr_set_str_binary (b, "0.100000000000000000000000000000000000000000000000000000000000000010000000000000000000000000000000000000000000000000000000000000001E65");
605 mpfr_set_str_binary (c, "0.10E1");
606 inex = mpfr_sub (a, b, c, MPFR_RNDN);
607 if (mpfr_cmp_ui_2exp (a, 1, 64) != 0 || inex >= 0)
608 {
609 printf ("Error in mpfr_sub for b-c for b=2^64+1+2^(-64), c=1\n");
610 printf ("Expected result 2^64 with inex < 0\n");
611 printf ("Got "); mpfr_dump (a);
612 printf ("with inex=%d\n", inex);
613 exit (1);
614 }
615
616 mpfr_clears (a, b, c, (mpfr_ptr) 0);
617 }
618
619 /* hard test of rounding */
620 static void
check_rounding(void)621 check_rounding (void)
622 {
623 mpfr_t a, b, c, res;
624 mpfr_prec_t p;
625 long k, l;
626 int i;
627
628 #define MAXKL (2 * GMP_NUMB_BITS)
629 for (p = MPFR_PREC_MIN; p <= GMP_NUMB_BITS; p++)
630 {
631 mpfr_init2 (a, p);
632 mpfr_init2 (res, p);
633 mpfr_init2 (b, p + 1 + MAXKL);
634 mpfr_init2 (c, MPFR_PREC_MIN);
635
636 /* b = 2^p + 1 + 2^(-k), c = 2^(-l) */
637 for (k = 0; k <= MAXKL; k++)
638 for (l = 0; l <= MAXKL; l++)
639 {
640 mpfr_set_ui_2exp (b, 1, p, MPFR_RNDN);
641 mpfr_add_ui (b, b, 1, MPFR_RNDN);
642 mpfr_mul_2ui (b, b, k, MPFR_RNDN);
643 mpfr_add_ui (b, b, 1, MPFR_RNDN);
644 mpfr_div_2ui (b, b, k, MPFR_RNDN);
645 mpfr_set_ui_2exp (c, 1, -l, MPFR_RNDN);
646 i = mpfr_sub (a, b, c, MPFR_RNDN);
647 /* b - c = 2^p + 1 + 2^(-k) - 2^(-l), should be rounded to
648 2^p for l <= k, and 2^p+2 for l < k, except when p=1 and
649 k=l, in which case b - c = 3, and the round-away rule implies
650 a = 4 = 2^p+2 = 2^(p+1) */
651 if (l < k || (l == k && p > 1))
652 {
653 if (mpfr_cmp_ui_2exp (a, 1, p) != 0)
654 {
655 printf ("Wrong result in check_rounding\n");
656 printf ("p=%lu k=%ld l=%ld\n", (unsigned long) p, k, l);
657 printf ("b="); mpfr_dump (b);
658 printf ("c="); mpfr_dump (c);
659 printf ("Expected 2^%lu\n", (unsigned long) p);
660 printf ("Got "); mpfr_dump (a);
661 exit (1);
662 }
663 if (i >= 0)
664 {
665 printf ("Wrong ternary value in check_rounding\n");
666 printf ("p=%lu k=%ld l=%ld\n", (unsigned long) p, k, l);
667 printf ("b="); mpfr_dump (b);
668 printf ("c="); mpfr_dump (c);
669 printf ("a="); mpfr_dump (a);
670 printf ("Expected < 0, got %d\n", i);
671 exit (1);
672 }
673 }
674 else /* l < k or (l = k and p = 1) */
675 {
676 mpfr_set_ui_2exp (res, 1, p, MPFR_RNDN);
677 mpfr_add_ui (res, res, 2, MPFR_RNDN);
678 if (mpfr_cmp (a, res) != 0)
679 {
680 printf ("Wrong result in check_rounding\n");
681 printf ("b="); mpfr_dump (b);
682 printf ("c="); mpfr_dump (c);
683 printf ("Expected "); mpfr_dump (res);
684 printf ("Got "); mpfr_dump (a);
685 exit (1);
686 }
687 if (i <= 0)
688 {
689 printf ("Wrong ternary value in check_rounding\n");
690 printf ("b="); mpfr_dump (b);
691 printf ("c="); mpfr_dump (c);
692 printf ("Expected > 0, got %d\n", i);
693 exit (1);
694 }
695 }
696 }
697
698 mpfr_clear (a);
699 mpfr_clear (res);
700 mpfr_clear (b);
701 mpfr_clear (c);
702 }
703 }
704
705 /* Check a = b - c, where the significand of b has all 1's, c is small
706 compared to b, and PREC(a) = PREC(b) - 1. Thus b is a midpoint for
707 the precision of the result a. The test is done with the extended
708 exponent range and with some reduced exponent range. Two choices
709 are made for the exponent of b: the maximum exponent - 1 (similar
710 to some normal case) and the maximum exponent (overflow case or
711 near overflow case, depending on the rounding mode).
712 This test is useful to trigger a bug in r10382: Since c is small,
713 the computation in sub1.c was done by first rounding b in the
714 precision of a, then correcting the result if b was a breakpoint
715 for this precision (exactly representable number for the directed
716 rounding modes, or midpoint for the round-to-nearest mode). The
717 problem was that for a midpoint in the round-to-nearest mode, the
718 rounding of b gave a spurious overflow; not only the overflow flag
719 was incorrect, but the result could not be corrected, since due to
720 this overflow, the "even rounding" information was lost.
721 In the case of reduced exponent range, an additional test is done
722 for consistency checks: the subtraction is done in the extended
723 exponent range (no overflow), then the result is converted to the
724 initial exponent range with mpfr_check_range. */
725 static void
check_max_almosteven(void)726 check_max_almosteven (void)
727 {
728 mpfr_exp_t old_emin, old_emax;
729 mpfr_exp_t emin[2] = { MPFR_EMIN_MIN, -1000 };
730 mpfr_exp_t emax[2] = { MPFR_EMAX_MAX, 1000 };
731 int i;
732
733 old_emin = mpfr_get_emin ();
734 old_emax = mpfr_get_emax ();
735
736 for (i = 0; i < 2; i++)
737 {
738 mpfr_t a1, a2, b, c;
739 mpfr_prec_t p;
740 int neg, j, rnd;
741
742 set_emin (emin[i]);
743 set_emax (emax[i]);
744
745 p = MPFR_PREC_MIN + randlimb () % 70;
746 mpfr_init2 (a1, p);
747 mpfr_init2 (a2, p);
748 mpfr_init2 (b, p+1);
749 mpfr_init2 (c, MPFR_PREC_MIN);
750
751 mpfr_setmax (b, 0);
752 mpfr_set_ui (c, 1, MPFR_RNDN);
753
754 for (neg = 0; neg < 2; neg++)
755 {
756 for (j = 1; j >= 0; j--)
757 {
758 mpfr_set_exp (b, __gmpfr_emax - j);
759 RND_LOOP_NO_RNDF (rnd)
760 {
761 mpfr_flags_t flags1, flags2;
762 int inex1, inex2;
763
764 /* Expected result. */
765 flags1 = MPFR_FLAGS_INEXACT;
766 if (rnd == MPFR_RNDN || MPFR_IS_LIKE_RNDZ (rnd, neg))
767 {
768 inex1 = neg ? 1 : -1;
769 mpfr_setmax (a1, __gmpfr_emax - j);
770 }
771 else
772 {
773 inex1 = neg ? -1 : 1;
774 if (j == 0)
775 {
776 flags1 |= MPFR_FLAGS_OVERFLOW;
777 mpfr_set_inf (a1, 1);
778 }
779 else
780 {
781 mpfr_setmin (a1, __gmpfr_emax);
782 }
783 }
784 MPFR_SET_SIGN (a1, neg ? -1 : 1);
785
786 /* Computed result. */
787 mpfr_clear_flags ();
788 inex2 = mpfr_sub (a2, b, c, (mpfr_rnd_t) rnd);
789 flags2 = __gmpfr_flags;
790
791 if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
792 mpfr_equal_p (a1, a2)))
793 {
794 printf ("Error 1 in check_max_almosteven for %s,"
795 " i = %d, j = %d, neg = %d\n",
796 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd),
797 i, j, neg);
798 printf (" b = ");
799 mpfr_dump (b);
800 printf ("Expected ");
801 mpfr_dump (a1);
802 printf (" with inex = %d, flags =", inex1);
803 flags_out (flags1);
804 printf ("Got ");
805 mpfr_dump (a2);
806 printf (" with inex = %d, flags =", inex2);
807 flags_out (flags2);
808 exit (1);
809 }
810
811 if (i == 0)
812 break;
813
814 /* Additional test for the reduced exponent range. */
815 mpfr_clear_flags ();
816 set_emin (MPFR_EMIN_MIN);
817 set_emax (MPFR_EMAX_MAX);
818 inex2 = mpfr_sub (a2, b, c, (mpfr_rnd_t) rnd);
819 set_emin (emin[i]);
820 set_emax (emax[i]);
821 inex2 = mpfr_check_range (a2, inex2, (mpfr_rnd_t) rnd);
822 flags2 = __gmpfr_flags;
823
824 if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
825 mpfr_equal_p (a1, a2)))
826 {
827 printf ("Error 2 in check_max_almosteven for %s,"
828 " i = %d, j = %d, neg = %d\n",
829 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd),
830 i, j, neg);
831 printf (" b = ");
832 mpfr_dump (b);
833 printf ("Expected ");
834 mpfr_dump (a1);
835 printf (" with inex = %d, flags =", inex1);
836 flags_out (flags1);
837 printf ("Got ");
838 mpfr_dump (a2);
839 printf (" with inex = %d, flags =", inex2);
840 flags_out (flags2);
841 exit (1);
842 }
843 }
844 } /* j */
845
846 mpfr_neg (b, b, MPFR_RNDN);
847 mpfr_neg (c, c, MPFR_RNDN);
848 } /* neg */
849
850 mpfr_clears (a1, a2, b, c, (mpfr_ptr) 0);
851 } /* i */
852
853 set_emin (old_emin);
854 set_emax (old_emax);
855 }
856
857 static void
test_rndf(void)858 test_rndf (void)
859 {
860 mpfr_t a, b, c, d;
861
862 mpfr_init2 (a, 7);
863 mpfr_init2 (b, 7);
864 mpfr_init2 (c, 7);
865 mpfr_init2 (d, 7);
866 mpfr_set_str_binary (b, "-1.000000e-7");
867 mpfr_set_str_binary (c, "-1.000000");
868 mpfr_sub (a, b, c, MPFR_RNDF);
869 MPFR_ASSERTN(MPFR_IS_NORMALIZED(a));
870 mpfr_sub (d, b, c, MPFR_RNDD);
871 if (!mpfr_equal_p (a, d))
872 {
873 mpfr_sub (d, b, c, MPFR_RNDU);
874 if (!mpfr_equal_p (a, d))
875 {
876 printf ("Error: mpfr_sub(a,b,c,RNDF) does not match RNDD/RNDU\n");
877 printf ("b="); mpfr_dump (b);
878 printf ("c="); mpfr_dump (c);
879 printf ("a="); mpfr_dump (a);
880 exit (1);
881 }
882 }
883 mpfr_clear (a);
884 mpfr_clear (b);
885 mpfr_clear (c);
886 mpfr_clear (d);
887 }
888
889 static void
testall_rndf(mpfr_prec_t pmax)890 testall_rndf (mpfr_prec_t pmax)
891 {
892 mpfr_t a, b, c, d;
893 mpfr_prec_t pa, pb, pc;
894 mpfr_exp_t eb;
895
896 for (pa = MPFR_PREC_MIN; pa <= pmax; pa++)
897 {
898 mpfr_init2 (a, pa);
899 mpfr_init2 (d, pa);
900 for (pb = MPFR_PREC_MIN; pb <= pmax; pb++)
901 {
902 mpfr_init2 (b, pb);
903 for (eb = 0; eb <= pmax + 3; eb ++)
904 {
905 mpfr_set_ui_2exp (b, 1, eb, MPFR_RNDN);
906 while (mpfr_cmp_ui_2exp (b, 1, eb + 1) < 0)
907 {
908 for (pc = MPFR_PREC_MIN; pc <= pmax; pc++)
909 {
910 mpfr_init2 (c, pc);
911 mpfr_set_ui (c, 1, MPFR_RNDN);
912 while (mpfr_cmp_ui (c, 2) < 0)
913 {
914 mpfr_sub (a, b, c, MPFR_RNDF);
915 mpfr_sub (d, b, c, MPFR_RNDD);
916 if (!mpfr_equal_p (a, d))
917 {
918 mpfr_sub (d, b, c, MPFR_RNDU);
919 if (!mpfr_equal_p (a, d))
920 {
921 printf ("Error: mpfr_sub(a,b,c,RNDF) does not "
922 "match RNDD/RNDU\n");
923 printf ("b="); mpfr_dump (b);
924 printf ("c="); mpfr_dump (c);
925 printf ("a="); mpfr_dump (a);
926 exit (1);
927 }
928 }
929 mpfr_nextabove (c);
930 }
931 mpfr_clear (c);
932 }
933 mpfr_nextabove (b);
934 }
935 }
936 mpfr_clear (b);
937 }
938 mpfr_clear (a);
939 mpfr_clear (d);
940 }
941 }
942
943 static void
test_rndf_exact(mp_size_t pmax)944 test_rndf_exact (mp_size_t pmax)
945 {
946 mpfr_t a, b, c, d;
947 mpfr_prec_t pa, pb, pc;
948 mpfr_exp_t eb;
949
950 for (pa = MPFR_PREC_MIN; pa <= pmax; pa++)
951 {
952 /* only check pa mod GMP_NUMB_BITS = -2, -1, 0, 1, 2 */
953 if ((pa + 2) % GMP_NUMB_BITS > 4)
954 continue;
955 mpfr_init2 (a, pa);
956 mpfr_init2 (d, pa);
957 for (pb = MPFR_PREC_MIN; pb <= pmax; pb++)
958 {
959 if ((pb + 2) % GMP_NUMB_BITS > 4)
960 continue;
961 mpfr_init2 (b, pb);
962 for (eb = 0; eb <= pmax + 3; eb ++)
963 {
964 mpfr_urandomb (b, RANDS);
965 mpfr_mul_2ui (b, b, eb, MPFR_RNDN);
966 for (pc = MPFR_PREC_MIN; pc <= pmax; pc++)
967 {
968 if ((pc + 2) % GMP_NUMB_BITS > 4)
969 continue;
970 mpfr_init2 (c, pc);
971 mpfr_urandomb (c, RANDS);
972 mpfr_sub (a, b, c, MPFR_RNDF);
973 mpfr_sub (d, b, c, MPFR_RNDD);
974 if (!mpfr_equal_p (a, d))
975 {
976 mpfr_sub (d, b, c, MPFR_RNDU);
977 if (!mpfr_equal_p (a, d))
978 {
979 printf ("Error: mpfr_sub(a,b,c,RNDF) does not "
980 "match RNDD/RNDU\n");
981 printf ("b="); mpfr_dump (b);
982 printf ("c="); mpfr_dump (c);
983 printf ("a="); mpfr_dump (a);
984 exit (1);
985 }
986 }
987
988 /* now make the low bits from c match those from b */
989 mpfr_add (c, b, d, MPFR_RNDN);
990 mpfr_sub (a, b, c, MPFR_RNDF);
991 mpfr_sub (d, b, c, MPFR_RNDD);
992 if (!mpfr_equal_p (a, d))
993 {
994 mpfr_sub (d, b, c, MPFR_RNDU);
995 if (!mpfr_equal_p (a, d))
996 {
997 printf ("Error: mpfr_sub(a,b,c,RNDF) does not "
998 "match RNDD/RNDU\n");
999 printf ("b="); mpfr_dump (b);
1000 printf ("c="); mpfr_dump (c);
1001 printf ("a="); mpfr_dump (a);
1002 exit (1);
1003 }
1004 }
1005
1006 mpfr_clear (c);
1007 }
1008 }
1009 mpfr_clear (b);
1010 }
1011 mpfr_clear (a);
1012 mpfr_clear (d);
1013 }
1014 }
1015
1016 /* Bug in the case 2 <= d < p in generic code mpfr_sub1sp() introduced
1017 * in r12242. Before this change, the special case that is failing was
1018 * handled by the MPFR_UNLIKELY(ap[n-1] == MPFR_LIMB_HIGHBIT) in the
1019 * "truncate:" code.
1020 */
1021 static void
bug20180215(void)1022 bug20180215 (void)
1023 {
1024 mpfr_t x, y, z1, z2;
1025 mpfr_rnd_t r[] = { MPFR_RNDN, MPFR_RNDU, MPFR_RNDA };
1026 int i, p;
1027
1028 for (p = 3; p <= 3 + 4 * GMP_NUMB_BITS; p++)
1029 {
1030 mpfr_inits2 (p, x, y, z1, z2, (mpfr_ptr) 0);
1031 mpfr_set_ui_2exp (x, 1, p - 1, MPFR_RNDN);
1032 mpfr_nextabove (x);
1033 mpfr_set_ui_2exp (y, 3, -1, MPFR_RNDN);
1034 mpfr_set (z1, x, MPFR_RNDN);
1035 mpfr_nextbelow (z1);
1036 mpfr_nextbelow (z1);
1037 for (i = 0; i < numberof (r); i++)
1038 {
1039 test_sub (z2, x, y, r[i]);
1040 if (! mpfr_equal_p (z1, z2))
1041 {
1042 printf ("Error in bug20180215 in precision %d, %s\n",
1043 p, mpfr_print_rnd_mode (r[i]));
1044 printf ("expected "); mpfr_dump (z1);
1045 printf ("got "); mpfr_dump (z2);
1046 exit (1);
1047 }
1048 }
1049 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
1050 }
1051 }
1052
1053 static void
bug20180216(void)1054 bug20180216 (void)
1055 {
1056 mpfr_t x, y, z1, z2;
1057 int r, p, d, inex;
1058
1059 for (p = 3; p <= 3 + 4 * GMP_NUMB_BITS; p++)
1060 {
1061 mpfr_inits2 (p, x, y, z1, z2, (mpfr_ptr) 0);
1062 for (d = 1; d <= p-2; d++)
1063 {
1064 inex = mpfr_set_ui_2exp (z1, 1, d, MPFR_RNDN); /* z1 = 2^d */
1065 MPFR_ASSERTN (inex == 0);
1066 inex = mpfr_add_ui (x, z1, 1, MPFR_RNDN);
1067 MPFR_ASSERTN (inex == 0);
1068 mpfr_nextabove (x); /* x = 2^d + 1 + epsilon */
1069 inex = mpfr_sub (y, x, z1, MPFR_RNDN); /* y = 1 + epsilon */
1070 MPFR_ASSERTN (inex == 0);
1071 inex = mpfr_add (z2, y, z1, MPFR_RNDN);
1072 MPFR_ASSERTN (inex == 0);
1073 MPFR_ASSERTN (mpfr_equal_p (z2, x)); /* consistency check */
1074 RND_LOOP (r)
1075 {
1076 inex = test_sub (z2, x, y, (mpfr_rnd_t) r);
1077 if (! mpfr_equal_p (z1, z2) || inex != 0)
1078 {
1079 printf ("Error in bug20180216 with p=%d, d=%d, %s\n",
1080 p, d, mpfr_print_rnd_mode ((mpfr_rnd_t) r));
1081 printf ("expected "); mpfr_dump (z1);
1082 printf (" with inex = 0\n");
1083 printf ("got "); mpfr_dump (z2);
1084 printf (" with inex = %d\n", inex);
1085 exit (1);
1086 }
1087 }
1088 }
1089 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
1090 }
1091 }
1092
1093 /* Fails with r12281: "reuse of c error in b - c in MPFR_RNDN".
1094 *
1095 * If the fix in r10697 (2016-07-29) is reverted, this test also fails
1096 * (there were no non-regression tests for this bug until this one);
1097 * note that if --enable-assert=full is used, the error message is:
1098 * "sub1 & sub1sp return different values for MPFR_RNDN".
1099 */
1100 static void
bug20180217(void)1101 bug20180217 (void)
1102 {
1103 mpfr_t x, y, z1, z2;
1104 int r, p, d, i, inex1, inex2;
1105
1106 for (p = 3; p <= 3 + 4 * GMP_NUMB_BITS; p++)
1107 {
1108 mpfr_inits2 (p, x, y, z1, z2, (mpfr_ptr) 0);
1109 for (d = p; d <= p+4; d++)
1110 {
1111 mpfr_set_ui (x, 1, MPFR_RNDN);
1112 mpfr_set_ui_2exp (y, 1, -d, MPFR_RNDN);
1113 for (i = 0; i < 3; i++)
1114 {
1115 RND_LOOP_NO_RNDF (r)
1116 {
1117 mpfr_set (z1, x, MPFR_RNDN);
1118 if (d == p)
1119 {
1120 mpfr_nextbelow (z1);
1121 if (i == 0)
1122 inex1 = 0;
1123 else if (r == MPFR_RNDD || r == MPFR_RNDZ ||
1124 (r == MPFR_RNDN && i > 1))
1125 {
1126 mpfr_nextbelow (z1);
1127 inex1 = -1;
1128 }
1129 else
1130 inex1 = 1;
1131 }
1132 else if (r == MPFR_RNDD || r == MPFR_RNDZ ||
1133 (r == MPFR_RNDN && d == p+1 && i > 0))
1134 {
1135 mpfr_nextbelow (z1);
1136 inex1 = -1;
1137 }
1138 else
1139 inex1 = 1;
1140 inex2 = test_sub (z2, x, y, (mpfr_rnd_t) r);
1141 if (!(mpfr_equal_p (z1, z2) && SAME_SIGN (inex1, inex2)))
1142 {
1143 printf ("Error in bug20180217 with "
1144 "p=%d, d=%d, i=%d, %s\n", p, d, i,
1145 mpfr_print_rnd_mode ((mpfr_rnd_t) r));
1146 printf ("x = ");
1147 mpfr_dump (x);
1148 printf ("y = ");
1149 mpfr_dump (y);
1150 printf ("Expected "); mpfr_dump (z1);
1151 printf (" with inex = %d\n", inex1);
1152 printf ("Got "); mpfr_dump (z2);
1153 printf (" with inex = %d\n", inex2);
1154 exit (1);
1155 }
1156 }
1157 if (i == 0)
1158 mpfr_nextabove (y);
1159 else
1160 {
1161 if (p < 6)
1162 break;
1163 mpfr_nextbelow (y);
1164 mpfr_mul_ui (y, y, 25, MPFR_RNDD);
1165 mpfr_div_2ui (y, y, 4, MPFR_RNDN);
1166 }
1167 }
1168 }
1169 mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
1170 }
1171 }
1172
1173 /* Tests on UBF.
1174 *
1175 * Note: mpfr_sub1sp will never be used as it does not support UBF.
1176 * Thus there is no need to generate tests for both mpfr_sub1 and
1177 * mpfr_sub1sp.
1178 *
1179 * Note that mpfr_sub1 has a special branch "c small", where the second
1180 * argument c is sufficiently smaller than the ulp of the first argument
1181 * and the ulp of the result: MAX (aq, bq) + 2 <= diff_exp.
1182 * Tests should be done for both the main branch and this special branch
1183 * when this makes sense.
1184 */
1185 #define REXP 1024
1186
test_ubf_aux(void)1187 static void test_ubf_aux (void)
1188 {
1189 mpfr_ubf_t x[11];
1190 mpfr_ptr p[11];
1191 int ex[11];
1192 mpfr_t ee, y, z, w;
1193 int i, j, k, neg, inexact, rnd;
1194 const int kn = 2;
1195 mpfr_exp_t e[] =
1196 { MPFR_EXP_MIN, MPFR_EMIN_MIN, -REXP, 0,
1197 REXP, MPFR_EMAX_MAX, MPFR_EXP_MAX };
1198
1199 mpfr_init2 (ee, sizeof (mpfr_exp_t) * CHAR_BIT);
1200 mpfr_inits2 (64, y, z, (mpfr_ptr) 0);
1201 mpfr_init2 (w, 2);
1202
1203 for (i = 0; i < 11; i++)
1204 p[i] = (mpfr_ptr) x[i];
1205
1206 /* exact zero result, with small and large exponents */
1207 for (i = 0; i < 2; i++)
1208 {
1209 mpfr_init2 (p[i], 5 + (randlimb () % 128));
1210 mpfr_set_ui (p[i], 17, MPFR_RNDN);
1211 mpz_init (MPFR_ZEXP (p[i]));
1212 MPFR_SET_UBF (p[i]);
1213 }
1214 for (j = 0; j < numberof (e); j++)
1215 {
1216 inexact = mpfr_set_exp_t (ee, e[j], MPFR_RNDN);
1217 MPFR_ASSERTN (inexact == 0);
1218 inexact = mpfr_get_z (MPFR_ZEXP (p[0]), ee, MPFR_RNDN);
1219 MPFR_ASSERTN (inexact == 0);
1220 mpz_sub_ui (MPFR_ZEXP (p[0]), MPFR_ZEXP (p[0]), kn);
1221
1222 for (k = -kn; k <= kn; k++)
1223 {
1224 /* exponent: e[j] + k, with |k| <= kn */
1225 mpz_set (MPFR_ZEXP (p[1]), MPFR_ZEXP (p[0]));
1226
1227 for (neg = 0; neg <= 1; neg++)
1228 {
1229 RND_LOOP (rnd)
1230 {
1231 /* Note: x[0] and x[1] are equal MPFR numbers, but do not
1232 test mpfr_sub with arg2 == arg3 as pointers in order to
1233 skip potentially optimized mpfr_sub code. */
1234 inexact = mpfr_sub (z, p[0], p[1], (mpfr_rnd_t) rnd);
1235 if (inexact != 0 || MPFR_NOTZERO (z) ||
1236 (rnd != MPFR_RNDD ? MPFR_IS_NEG (z) : MPFR_IS_POS (z)))
1237 {
1238 printf ("Error 1 in test_ubf for exact zero result: "
1239 "j=%d k=%d neg=%d, rnd=%s\nGot ", j, k, neg,
1240 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
1241 mpfr_dump (z);
1242 printf ("inexact = %d\n", inexact);
1243 exit (1);
1244 }
1245 }
1246
1247 for (i = 0; i < 2; i++)
1248 MPFR_CHANGE_SIGN (p[i]);
1249 }
1250
1251 mpz_add_ui (MPFR_ZEXP (p[0]), MPFR_ZEXP (p[0]), 1);
1252 }
1253 }
1254 for (i = 0; i < 2; i++)
1255 {
1256 MPFR_UBF_CLEAR_EXP (p[i]);
1257 mpfr_clear (p[i]);
1258 }
1259
1260 /* Up to a given exponent (for the result) and sign, test:
1261 * (t + .11010) - (t + .00001) = .11001
1262 * (t + 8) - (t + 111.00111) = .11001
1263 * where t = 0 or a power of 2, e.g. 2^200. Test various exponents
1264 * (including those near the underflow/overflow boundaries) so that
1265 * the subtraction yields a normal number, an overflow or an underflow.
1266 * In MPFR_RNDA, also test with a 2-bit precision target, as this
1267 * yields an exponent change.
1268 *
1269 * Also test the "MAX (aq, bq) + 2 <= diff_exp" branch of sub1.c with
1270 * .1 - epsilon (possible decrease of the exponent) and .111 - epsilon
1271 * in precision 2 (possible increase of the exponent). The first test
1272 * triggers a possible decrease of the exponent (see bug fixed in r13806).
1273 * The second test triggers a possible increase of the exponent (see the
1274 * "exp_a != MPFR_EXP_MAX" test to avoid an integer overflow).
1275 */
1276 for (i = 0; i < 8; i++)
1277 {
1278 static int v[4] = { 26, 1, 256, 231 };
1279
1280 mpfr_init2 (p[i], i < 4 ? 8 + (randlimb () % 128) : 256);
1281 if (i < 4)
1282 {
1283 inexact = mpfr_set_si_2exp (p[i], v[i], -5, MPFR_RNDN);
1284 MPFR_ASSERTN (inexact == 0);
1285 }
1286 else
1287 {
1288 inexact = mpfr_set_si_2exp (p[i], 1, 200, MPFR_RNDN);
1289 MPFR_ASSERTN (inexact == 0);
1290 inexact = mpfr_add (p[i], p[i], p[i-4], MPFR_RNDN);
1291 MPFR_ASSERTN (inexact == 0);
1292 }
1293 ex[i] = mpfr_get_exp (p[i]) + 5;
1294 MPFR_ASSERTN (ex[i] >= 0);
1295 }
1296 mpfr_inits2 (3, p[8], p[9], p[10], (mpfr_ptr) 0);
1297 inexact = mpfr_set_si_2exp (p[8], 1, 0, MPFR_RNDN);
1298 MPFR_ASSERTN (inexact == 0);
1299 ex[8] = 5;
1300 inexact = mpfr_set_si_2exp (p[9], 1, 0, MPFR_RNDN); /* will be epsilon */
1301 MPFR_ASSERTN (inexact == 0);
1302 ex[9] = 0;
1303 inexact = mpfr_set_si_2exp (p[10], 7, 0, MPFR_RNDN);
1304 MPFR_ASSERTN (inexact == 0);
1305 ex[10] = 5;
1306
1307 for (i = 0; i < 11; i++)
1308 {
1309 mpz_init (MPFR_ZEXP (p[i]));
1310 MPFR_SET_UBF (p[i]);
1311 }
1312
1313 for (j = 0; j < numberof (e); j++)
1314 {
1315 inexact = mpfr_set_exp_t (ee, e[j], MPFR_RNDN);
1316 MPFR_ASSERTN (inexact == 0);
1317 inexact = mpfr_get_z (MPFR_ZEXP (p[0]), ee, MPFR_RNDN);
1318 MPFR_ASSERTN (inexact == 0);
1319 for (i = 1; i < 11; i++)
1320 mpz_set (MPFR_ZEXP (p[i]), MPFR_ZEXP (p[0]));
1321 for (i = 0; i < 11; i++)
1322 {
1323 mpz_add_ui (MPFR_ZEXP (p[i]), MPFR_ZEXP (p[i]), ex[i]);
1324 mpz_sub_ui (MPFR_ZEXP (p[i]), MPFR_ZEXP (p[i]), 5 + kn);
1325 }
1326 mpz_sub_ui (MPFR_ZEXP (p[9]), MPFR_ZEXP (p[9]), 256);
1327 for (k = -kn; k <= kn; k++)
1328 {
1329 for (neg = 0; neg <= 1; neg++)
1330 {
1331 int sign = neg ? -1 : 1;
1332
1333 RND_LOOP (rnd)
1334 for (i = 0; i <= 10; i += 2)
1335 {
1336 mpfr_exp_t e0;
1337 mpfr_flags_t flags, flags_y;
1338 int inex_y;
1339
1340 if (i >= 8)
1341 {
1342 int d;
1343
1344 e0 = MPFR_UBF_GET_EXP (p[i]);
1345 if (e0 < MPFR_EXP_MIN + 3)
1346 e0 += 3;
1347
1348 if (rnd == MPFR_RNDN)
1349 d = i == 8 ? (e0 == __gmpfr_emin - 1 ? 3 : 4) : 6;
1350 else if (MPFR_IS_LIKE_RNDZ (rnd, neg))
1351 d = i == 8 ? 3 : 6;
1352 else
1353 d = i == 8 ? 4 : 8;
1354
1355 mpfr_clear_flags ();
1356 inex_y = mpfr_set_si_2exp (w, sign * d, e0 - 3,
1357 (mpfr_rnd_t) rnd);
1358 flags_y = __gmpfr_flags | MPFR_FLAGS_INEXACT;
1359 if (inex_y == 0)
1360 inex_y = rnd == MPFR_RNDN ?
1361 sign * (i == 8 ? 1 : -1) :
1362 MPFR_IS_LIKE_RNDD ((mpfr_rnd_t) rnd, sign) ?
1363 -1 : 1;
1364 mpfr_set (y, w, MPFR_RNDN);
1365
1366 mpfr_clear_flags ();
1367 inexact = mpfr_sub (w, p[i], p[9], (mpfr_rnd_t) rnd);
1368 flags = __gmpfr_flags;
1369
1370 /* For MPFR_RNDF, only do a basic test. */
1371 MPFR_ASSERTN (mpfr_check (w));
1372 if (rnd == MPFR_RNDF)
1373 continue;
1374
1375 goto testw;
1376 }
1377
1378 mpfr_clear_flags ();
1379 inexact = mpfr_sub (z, p[i], p[i+1], (mpfr_rnd_t) rnd);
1380 flags = __gmpfr_flags;
1381
1382 /* For MPFR_RNDF, only do a basic test. */
1383 MPFR_ASSERTN (mpfr_check (z));
1384 if (rnd == MPFR_RNDF)
1385 continue;
1386
1387 e0 = MPFR_UBF_GET_EXP (p[0]);
1388
1389 if (e0 < __gmpfr_emin)
1390 {
1391 mpfr_rnd_t r =
1392 rnd == MPFR_RNDN && e0 < __gmpfr_emin - 1 ?
1393 MPFR_RNDZ : (mpfr_rnd_t) rnd;
1394 flags_y = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
1395 inex_y = mpfr_underflow (y, r, sign);
1396 }
1397 else if (e0 > __gmpfr_emax)
1398 {
1399 flags_y = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT;
1400 inex_y = mpfr_overflow (y, (mpfr_rnd_t) rnd, sign);
1401 }
1402 else
1403 {
1404 mpfr_set_si_2exp (y, sign * 25, e0 - 5, MPFR_RNDN);
1405 flags_y = 0;
1406 inex_y = 0;
1407 }
1408
1409 if (flags != flags_y ||
1410 ! SAME_SIGN (inexact, inex_y) ||
1411 ! mpfr_equal_p (y, z))
1412 {
1413 printf ("Error 2 in test_ubf with "
1414 "j=%d k=%d neg=%d i=%d rnd=%s\n",
1415 j, k, neg, i,
1416 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
1417 printf ("emin=%" MPFR_EXP_FSPEC "d "
1418 "emax=%" MPFR_EXP_FSPEC "d\n",
1419 (mpfr_eexp_t) __gmpfr_emin,
1420 (mpfr_eexp_t) __gmpfr_emax);
1421 printf ("b = ");
1422 mpfr_dump (p[i]);
1423 printf ("c = ");
1424 mpfr_dump (p[i+1]);
1425 printf ("Expected ");
1426 mpfr_dump (y);
1427 printf ("with inex = %d and flags =", inex_y);
1428 flags_out (flags_y);
1429 printf ("Got ");
1430 mpfr_dump (z);
1431 printf ("with inex = %d and flags =", inexact);
1432 flags_out (flags);
1433 exit (1);
1434 }
1435
1436 /* Do the following 2-bit precision test only in RNDA. */
1437 if (rnd != MPFR_RNDA)
1438 continue;
1439
1440 mpfr_clear_flags ();
1441 inexact = mpfr_sub (w, p[i], p[i+1], MPFR_RNDA);
1442 flags = __gmpfr_flags;
1443 if (e0 < MPFR_EXP_MAX)
1444 e0++;
1445
1446 if (e0 < __gmpfr_emin)
1447 {
1448 flags_y = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
1449 inex_y = mpfr_underflow (y, MPFR_RNDA, sign);
1450 }
1451 else if (e0 > __gmpfr_emax)
1452 {
1453 flags_y = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT;
1454 inex_y = mpfr_overflow (y, MPFR_RNDA, sign);
1455 }
1456 else
1457 {
1458 mpfr_set_si_2exp (y, sign, e0 - 1, MPFR_RNDN);
1459 flags_y = MPFR_FLAGS_INEXACT;
1460 inex_y = sign;
1461 }
1462
1463 testw:
1464 if (flags != flags_y ||
1465 ! SAME_SIGN (inexact, inex_y) ||
1466 ! mpfr_equal_p (y, w))
1467 {
1468 printf ("Error 3 in test_ubf with "
1469 "j=%d k=%d neg=%d i=%d rnd=%s\n",
1470 j, k, neg, i,
1471 mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
1472 printf ("emin=%" MPFR_EXP_FSPEC "d "
1473 "emax=%" MPFR_EXP_FSPEC "d\n",
1474 (mpfr_eexp_t) __gmpfr_emin,
1475 (mpfr_eexp_t) __gmpfr_emax);
1476 printf ("b = ");
1477 mpfr_dump (p[i]);
1478 printf ("c = ");
1479 mpfr_dump (p[i <= 8 ? i+1 : 9]);
1480 printf ("Expected ");
1481 /* Set y to a 2-bit precision just for the output.
1482 Since we exit, this will have no other effect. */
1483 mpfr_prec_round (y, 2, MPFR_RNDA);
1484 mpfr_dump (y);
1485 printf ("with inex = %d and flags =", inex_y);
1486 flags_out (flags_y);
1487 printf ("Got ");
1488 mpfr_dump (w);
1489 printf ("with inex = %d and flags =", inexact);
1490 flags_out (flags);
1491 exit (1);
1492 }
1493 }
1494
1495 for (i = 0; i < 11; i++)
1496 MPFR_CHANGE_SIGN (p[i]);
1497 }
1498
1499 for (i = 0; i < 11; i++)
1500 mpz_add_ui (MPFR_ZEXP (p[i]), MPFR_ZEXP (p[i]), 1);
1501 }
1502 }
1503 for (i = 0; i < 11; i++)
1504 {
1505 MPFR_UBF_CLEAR_EXP (p[i]);
1506 mpfr_clear (p[i]);
1507 }
1508
1509 mpfr_clears (ee, y, z, w, (mpfr_ptr) 0);
1510 }
1511
1512 /* Run the tests on UBF with the maximum exponent range and with a
1513 reduced exponent range. */
test_ubf(void)1514 static void test_ubf (void)
1515 {
1516 mpfr_exp_t emin, emax;
1517
1518 emin = mpfr_get_emin ();
1519 emax = mpfr_get_emax ();
1520
1521 set_emin (MPFR_EMIN_MIN);
1522 set_emax (MPFR_EMAX_MAX);
1523 test_ubf_aux ();
1524
1525 set_emin (-REXP);
1526 set_emax (REXP);
1527 test_ubf_aux ();
1528
1529 set_emin (emin);
1530 set_emax (emax);
1531 }
1532
1533 #define TEST_FUNCTION test_sub
1534 #define TWO_ARGS
1535 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
1536 #include "tgeneric.c"
1537
1538 int
main(void)1539 main (void)
1540 {
1541 mpfr_prec_t p;
1542 unsigned int i;
1543
1544 tests_start_mpfr ();
1545
1546 test_rndf ();
1547 testall_rndf (7);
1548 test_rndf_exact (200);
1549 bug20101017 ();
1550 check_rounding ();
1551 check_diverse ();
1552 check_inexact ();
1553 check_max_almosteven ();
1554 bug_ddefour ();
1555 bug20180215 ();
1556 bug20180216 ();
1557 bug20180217 ();
1558 for (p=2; p<200; p++)
1559 for (i=0; i<50; i++)
1560 check_two_sum (p);
1561 test_generic (MPFR_PREC_MIN, 800, 100);
1562 test_ubf ();
1563
1564 tests_end_mpfr ();
1565 return 0;
1566 }
1567