tests/tsub.c

   1 /* Test file for mpfr_sub.
   2
   3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramel projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
  20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  22
  23 #include <stdio.h>
  24 #include <stdlib.h>
  25
  26 #include "mpfr-test.h"
  27
  28 #ifdef CHECK_EXTERNAL
  29 static int
  30 test_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
  31 {
  32   int res;
  33   int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c);
  34   if (ok)
  35     {
  36       mpfr_print_raw (b);
  37       printf (" ");
  38       mpfr_print_raw (c);
  39     }
  40   res = mpfr_sub (a, b, c, rnd_mode);
  41   if (ok)
  42     {
  43       printf (" ");
  44       mpfr_print_raw (a);
  45       printf ("\n");
  46     }
  47   return res;
  48 }
  49 #else
  50 #define test_sub mpfr_sub
  51 #endif
  52
  53 static void
  54 check_diverse (void)
  55 {
  56   mpfr_t x, y, z;
  57   int inexact;
  58
  59   mpfr_init (x);
  60   mpfr_init (y);
  61   mpfr_init (z);
  62
  63   /* check corner case cancel=0, but add_exp=1 */
  64   mpfr_set_prec (x, 2);
  65   mpfr_set_prec (y, 4);
  66   mpfr_set_prec (z, 2);
  67   mpfr_setmax (y, __gmpfr_emax);
  68   mpfr_set_str_binary (z, "0.1E-10"); /* tiny */
  69   test_sub (x, y, z, MPFR_RNDN); /* should round to 2^emax, i.e. overflow */
  70   if (!mpfr_inf_p (x) || mpfr_sgn (x) < 0)
  71     {
  72       printf ("Error in mpfr_sub(a,b,c,RNDN) for b=maxfloat, prec(a)<prec(b), c tiny\n");
  73       exit (1);
  74     }
  75
  76   /* other coverage test */
  77   mpfr_set_prec (x, 2);
  78   mpfr_set_prec (y, 2);
  79   mpfr_set_prec (z, 2);
  80   mpfr_set_ui (y, 1, MPFR_RNDN);
  81   mpfr_set_si (z, -2, MPFR_RNDN);
  82   test_sub (x, y, z, MPFR_RNDD);
  83   if (mpfr_cmp_ui (x, 3))
  84     {
  85       printf ("Error in mpfr_sub(1,-2,RNDD)\n");
  86       exit (1);
  87     }
  88
  89   mpfr_set_prec (x, 288);
  90   mpfr_set_prec (y, 288);
  91   mpfr_set_prec (z, 288);
  92   mpfr_set_str_binary (y, "0.111000110011000001000111101010111011110011101001101111111110000011100101000001001010110010101010011001010100000001110011110001010101101010001011101110100100001011110100110000101101100011010001001011011010101010000010001101001000110010010111111011110001111101001000101101001100101100101000E80");
  93   mpfr_set_str_binary (z, "0.100001111111101001011010001100110010100111001110000110011101001011010100001000000100111011010110110010000000000010101101011000010000110001110010100001100101011100100100001011000100011110000001010101000100011101001000010111100000111000111011001000100100011000100000010010111000000100100111E-258");
  94   inexact = test_sub (x, y, z, MPFR_RNDN);
  95   if (inexact <= 0)
  96     {
  97       printf ("Wrong inexact flag for prec=288\n");
  98       exit (1);
  99     }
 100
 101   mpfr_set_prec (x, 32);
 102   mpfr_set_prec (y, 63);
 103   mpfr_set_prec (z, 63);
 104   mpfr_set_str_binary (x, "0.101101111011011100100100100111E31");
 105   mpfr_set_str_binary (y, "0.111110010010100100110101101010001001100101110001000101110111111E-1");
 106   test_sub (z, x, y, MPFR_RNDN);
 107   mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
 108   if (mpfr_cmp (z, y))
 109     {
 110       printf ("Error in mpfr_sub (5)\n");
 111       printf ("expected "); mpfr_print_binary (y); puts ("");
 112       printf ("got      "); mpfr_print_binary (z); puts ("");
 113       exit (1);
 114     }
 115
 116   mpfr_set_prec (y, 63);
 117   mpfr_set_prec (z, 63);
 118   mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
 119   mpfr_sub_ui (z, y, 1541116494, MPFR_RNDN);
 120   mpfr_set_str_binary (y, "-0.11111001001010010011010110101E-1");
 121   if (mpfr_cmp (z, y))
 122     {
 123       printf ("Error in mpfr_sub (7)\n");
 124       printf ("expected "); mpfr_print_binary (y); puts ("");
 125       printf ("got      "); mpfr_print_binary (z); puts ("");
 126       exit (1);
 127     }
 128
 129   mpfr_set_prec (y, 63);
 130   mpfr_set_prec (z, 63);
 131   mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
 132   mpfr_sub_ui (z, y, 1541116494, MPFR_RNDN);
 133   mpfr_set_str_binary (y, "-0.11111001001010010011010110101E-1");
 134   if (mpfr_cmp (z, y))
 135     {
 136       printf ("Error in mpfr_sub (6)\n");
 137       printf ("expected "); mpfr_print_binary (y); puts ("");
 138       printf ("got      "); mpfr_print_binary (z); puts ("");
 139       exit (1);
 140     }
 141
 142   mpfr_set_prec (x, 32);
 143   mpfr_set_prec (y, 32);
 144   mpfr_set_str_binary (x, "0.10110111101001110100100101111000E0");
 145   mpfr_set_str_binary (y, "0.10001100100101000100110111000100E0");
 146   if ((inexact = test_sub (x, x, y, MPFR_RNDN)))
 147     {
 148       printf ("Wrong inexact flag (2): got %d instead of 0\n", inexact);
 149       exit (1);
 150     }
 151
 152   mpfr_set_prec (x, 32);
 153   mpfr_set_prec (y, 32);
 154   mpfr_set_str_binary (x, "0.11111000110111011000100111011010E0");
 155   mpfr_set_str_binary (y, "0.10011111101111000100001000000000E-3");
 156   if ((inexact = test_sub (x, x, y, MPFR_RNDN)))
 157     {
 158       printf ("Wrong inexact flag (1): got %d instead of 0\n", inexact);
 159       exit (1);
 160     }
 161
 162   mpfr_set_prec (x, 33);
 163   mpfr_set_prec (y, 33);
 164   mpfr_set_ui (x, 1, MPFR_RNDN);
 165   mpfr_set_str_binary (y, "-0.1E-32");
 166   mpfr_add (x, x, y, MPFR_RNDN);
 167   mpfr_set_str_binary (y, "0.111111111111111111111111111111111E0");
 168   if (mpfr_cmp (x, y))
 169     {
 170       printf ("Error in mpfr_sub (1 - 1E-33) with prec=33\n");
 171       printf ("Expected "); mpfr_print_binary (y); puts ("");
 172       printf ("got      "); mpfr_print_binary (x); puts ("");
 173       exit (1);
 174     }
 175
 176   mpfr_set_prec (x, 32);
 177   mpfr_set_prec (y, 33);
 178   mpfr_set_ui (x, 1, MPFR_RNDN);
 179   mpfr_set_str_binary (y, "-0.1E-32");
 180   mpfr_add (x, x, y, MPFR_RNDN);
 181   if (mpfr_cmp_ui (x, 1))
 182     {
 183       printf ("Error in mpfr_sub (1 - 1E-33) with prec=32\n");
 184       printf ("Expected 1.0, got "); mpfr_print_binary (x); puts ("");
 185       exit (1);
 186     }
 187
 188   mpfr_set_prec (x, 65);
 189   mpfr_set_prec (y, 65);
 190   mpfr_set_prec (z, 64);
 191   mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
 192   mpfr_set_str_binary (y, "0.1110111011110001110111011111111111101000011001011100101100101100");
 193   test_sub (z, x, y, MPFR_RNDZ);
 194   if (mpfr_cmp_ui (z, 1))
 195     {
 196       printf ("Error in mpfr_sub (1)\n");
 197       exit (1);
 198     }
 199   test_sub (z, x, y, MPFR_RNDU);
 200   mpfr_set_str_binary (x, "1.000000000000000000000000000000000000000000000000000000000000001");
 201   if (mpfr_cmp (z, x))
 202     {
 203       printf ("Error in mpfr_sub (2)\n");
 204       printf ("Expected "); mpfr_print_binary (x); puts ("");
 205       printf ("Got      "); mpfr_print_binary (z); puts ("");
 206       exit (1);
 207     }
 208   mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
 209   test_sub (z, x, y, MPFR_RNDN);
 210   if (mpfr_cmp_ui (z, 1))
 211     {
 212       printf ("Error in mpfr_sub (3)\n");
 213       exit (1);
 214     }
 215   inexact = test_sub (z, x, y, MPFR_RNDA);
 216   mpfr_set_str_binary (x, "1.000000000000000000000000000000000000000000000000000000000000001");
 217   if (mpfr_cmp (z, x) || inexact <= 0)
 218     {
 219       printf ("Error in mpfr_sub (4)\n");
 220       exit (1);
 221     }
 222   mpfr_set_prec (x, 66);
 223   mpfr_set_str_binary (x, "1.11101110111100011101110111111111111010000110010111001011001010111");
 224   test_sub (z, x, y, MPFR_RNDN);
 225   if (mpfr_cmp_ui (z, 1))
 226     {
 227       printf ("Error in mpfr_sub (5)\n");
 228       exit (1);
 229     }
 230
 231   /* check in-place operations */
 232   mpfr_set_ui (x, 1, MPFR_RNDN);
 233   test_sub (x, x, x, MPFR_RNDN);
 234   if (mpfr_cmp_ui(x, 0))
 235     {
 236       printf ("Error for mpfr_sub (x, x, x, MPFR_RNDN) with x=1.0\n");
 237       exit (1);
 238     }
 239
 240   mpfr_set_prec (x, 53);
 241   mpfr_set_prec (y, 53);
 242   mpfr_set_prec (z, 53);
 243   mpfr_set_str1 (x, "1.229318102e+09");
 244   mpfr_set_str1 (y, "2.32221184180698677665e+05");
 245   test_sub (z, x, y, MPFR_RNDN);
 246   if (mpfr_cmp_str1 (z, "1229085880.815819263458251953125"))
 247     {
 248       printf ("Error in mpfr_sub (1.22e9 - 2.32e5)\n");
 249       printf ("expected 1229085880.815819263458251953125, got ");
 250       mpfr_out_str(stdout, 10, 0, z, MPFR_RNDN);
 251       putchar('\n');
 252       exit (1);
 253     }
 254
 255   mpfr_set_prec (x, 112);
 256   mpfr_set_prec (y, 98);
 257   mpfr_set_prec (z, 54);
 258   mpfr_set_str_binary (x, "0.11111100100000000011000011100000101101010001000111E-401");
 259   mpfr_set_str_binary (y, "0.10110000100100000101101100011111111011101000111000101E-464");
 260   test_sub (z, x, y, MPFR_RNDN);
 261   if (mpfr_cmp (z, x)) {
 262     printf ("mpfr_sub(z, x, y) failed for prec(x)=112, prec(y)=98\n");
 263     printf ("expected "); mpfr_print_binary (x); puts ("");
 264     printf ("got      "); mpfr_print_binary (z); puts ("");
 265     exit (1);
 266   }
 267
 268   mpfr_set_prec (x, 33);
 269   mpfr_set_ui (x, 1, MPFR_RNDN);
 270   mpfr_div_2exp (x, x, 32, MPFR_RNDN);
 271   mpfr_sub_ui (x, x, 1, MPFR_RNDN);
 272
 273   mpfr_set_prec (x, 5);
 274   mpfr_set_prec (y, 5);
 275   mpfr_set_str_binary (x, "1e-12");
 276   mpfr_set_ui (y, 1, MPFR_RNDN);
 277   test_sub (x, y, x, MPFR_RNDD);
 278   mpfr_set_str_binary (y, "0.11111");
 279   if (mpfr_cmp (x, y))
 280     {
 281       printf ("Error in mpfr_sub (x, y, x, MPFR_RNDD) for x=2^(-12), y=1\n");
 282       exit (1);
 283     }
 284
 285   mpfr_set_prec (x, 24);
 286   mpfr_set_prec (y, 24);
 287   mpfr_set_str_binary (x, "-0.100010000000000000000000E19");
 288   mpfr_set_str_binary (y, "0.100000000000000000000100E15");
 289   mpfr_add (x, x, y, MPFR_RNDD);
 290   mpfr_set_str_binary (y, "-0.1E19");
 291   if (mpfr_cmp (x, y))
 292     {
 293       printf ("Error in mpfr_add (2)\n");
 294       exit (1);
 295     }
 296
 297   mpfr_set_prec (x, 2);
 298   mpfr_set_prec (y, 10);
 299   mpfr_set_prec (z, 10);
 300   mpfr_set_ui (y, 0, MPFR_RNDN);
 301   mpfr_set_str_binary (z, "0.10001");
 302   if (test_sub (x, y, z, MPFR_RNDN) <= 0)
 303     {
 304       printf ("Wrong inexact flag in x=mpfr_sub(0,z) for prec(z)>prec(x)\n");
 305       exit (1);
 306     }
 307   if (test_sub (x, z, y, MPFR_RNDN) >= 0)
 308     {
 309       printf ("Wrong inexact flag in x=mpfr_sub(z,0) for prec(z)>prec(x)\n");
 310       exit (1);
 311     }
 312
 313   mpfr_clear (x);
 314   mpfr_clear (y);
 315   mpfr_clear (z);
 316 }
 317
 318 static void
 319 bug_ddefour(void)
 320 {
 321     mpfr_t ex, ex1, ex2, ex3, tot, tot1;
 322
 323     mpfr_init2(ex, 53);
 324     mpfr_init2(ex1, 53);
 325     mpfr_init2(ex2, 53);
 326     mpfr_init2(ex3, 53);
 327     mpfr_init2(tot, 150);
 328     mpfr_init2(tot1, 150);
 329
 330     mpfr_set_ui( ex, 1, MPFR_RNDN);
 331     mpfr_mul_2exp( ex, ex, 906, MPFR_RNDN);
 332     mpfr_log( tot, ex, MPFR_RNDN);
 333     mpfr_set( ex1, tot, MPFR_RNDN); /* ex1 = high(tot) */
 334     test_sub( ex2, tot, ex1, MPFR_RNDN); /* ex2 = high(tot - ex1) */
 335     test_sub( tot1, tot, ex1, MPFR_RNDN); /* tot1 = tot - ex1 */
 336     mpfr_set( ex3, tot1, MPFR_RNDN); /* ex3 = high(tot - ex1) */
 337
 338     if (mpfr_cmp(ex2, ex3))
 339       {
 340         printf ("Error in ddefour test.\n");
 341         printf ("ex2="); mpfr_print_binary (ex2); puts ("");
 342         printf ("ex3="); mpfr_print_binary (ex3); puts ("");
 343         exit (1);
 344       }
 345
 346     mpfr_clear (ex);
 347     mpfr_clear (ex1);
 348     mpfr_clear (ex2);
 349     mpfr_clear (ex3);
 350     mpfr_clear (tot);
 351     mpfr_clear (tot1);
 352 }
 353
 354 /* if u = o(x-y), v = o(u-x), w = o(v+y), then x-y = u-w */
 355 static void
 356 check_two_sum (mpfr_prec_t p)
 357 {
 358   mpfr_t x, y, u, v, w;
 359   mpfr_rnd_t rnd;
 360   int inexact;
 361
 362   mpfr_init2 (x, p);
 363   mpfr_init2 (y, p);
 364   mpfr_init2 (u, p);
 365   mpfr_init2 (v, p);
 366   mpfr_init2 (w, p);
 367   mpfr_urandomb (x, RANDS);
 368   mpfr_urandomb (y, RANDS);
 369   if (mpfr_cmpabs (x, y) < 0)
 370     mpfr_swap (x, y);
 371   rnd = MPFR_RNDN;
 372   inexact = test_sub (u, x, y, rnd);
 373   test_sub (v, u, x, rnd);
 374   mpfr_add (w, v, y, rnd);
 375   /* as u = (x-y) - w, we should have inexact and w of opposite signs */
 376   if (((inexact == 0) && mpfr_cmp_ui (w, 0)) ||
 377       ((inexact > 0) && (mpfr_cmp_ui (w, 0) <= 0)) ||
 378       ((inexact < 0) && (mpfr_cmp_ui (w, 0) >= 0)))
 379     {
 380       printf ("Wrong inexact flag for prec=%u, rnd=%s\n", (unsigned)p,
 381                mpfr_print_rnd_mode (rnd));
 382       printf ("x="); mpfr_print_binary(x); puts ("");
 383       printf ("y="); mpfr_print_binary(y); puts ("");
 384       printf ("u="); mpfr_print_binary(u); puts ("");
 385       printf ("v="); mpfr_print_binary(v); puts ("");
 386       printf ("w="); mpfr_print_binary(w); puts ("");
 387       printf ("inexact = %d\n", inexact);
 388       exit (1);
 389     }
 390   mpfr_clear (x);
 391   mpfr_clear (y);
 392   mpfr_clear (u);
 393   mpfr_clear (v);
 394   mpfr_clear (w);
 395 }
 396
 397 #define MAX_PREC 200
 398
 399 static void
 400 check_inexact (void)
 401 {
 402   mpfr_t x, y, z, u;
 403   mpfr_prec_t px, py, pu, pz;
 404   int inexact, cmp;
 405   mpfr_rnd_t rnd;
 406
 407   mpfr_init (x);
 408   mpfr_init (y);
 409   mpfr_init (z);
 410   mpfr_init (u);
 411
 412   mpfr_set_prec (x, 2);
 413   mpfr_set_ui (x, 6, MPFR_RNDN);
 414   mpfr_div_2exp (x, x, 4, MPFR_RNDN); /* x = 6/16 */
 415   mpfr_set_prec (y, 2);
 416   mpfr_set_si (y, -1, MPFR_RNDN);
 417   mpfr_div_2exp (y, y, 4, MPFR_RNDN); /* y = -1/16 */
 418   inexact = test_sub (y, y, x, MPFR_RNDN); /* y = round(-7/16) = -1/2 */
 419   if (inexact >= 0)
 420     {
 421       printf ("Error: wrong inexact flag for -1/16 - (6/16)\n");
 422       exit (1);
 423     }
 424
 425   for (px=2; px<MAX_PREC; px++)
 426     {
 427       mpfr_set_prec (x, px);
 428       do
 429         {
 430           mpfr_urandomb (x, RANDS);
 431         }
 432       while (mpfr_cmp_ui (x, 0) == 0);
 433       for (pu=2; pu<MAX_PREC; pu++)
 434         {
 435           mpfr_set_prec (u, pu);
 436           do
 437             {
 438               mpfr_urandomb (u, RANDS);
 439             }
 440           while (mpfr_cmp_ui (u, 0) == 0);
 441           {
 442               py = 2 + (randlimb () % (MAX_PREC - 2));
 443               mpfr_set_prec (y, py);
 444               /* warning: MPFR_EXP is undefined for 0 */
 445               pz =  (mpfr_cmpabs (x, u) >= 0) ? MPFR_EXP(x) - MPFR_EXP(u)
 446                 : MPFR_EXP(u) - MPFR_EXP(x);
 447               pz = pz + MAX(MPFR_PREC(x), MPFR_PREC(u));
 448               mpfr_set_prec (z, pz);
 449               rnd = RND_RAND ();
 450               if (test_sub (z, x, u, rnd))
 451                 {
 452                   printf ("z <- x - u should be exact\n");
 453                   exit (1);
 454                 }
 455                 {
 456                   rnd = RND_RAND ();
 457                   inexact = test_sub (y, x, u, rnd);
 458                   cmp = mpfr_cmp (y, z);
 459                   if (((inexact == 0) && (cmp != 0)) ||
 460                       ((inexact > 0) && (cmp <= 0)) ||
 461                       ((inexact < 0) && (cmp >= 0)))
 462                     {
 463                       printf ("Wrong inexact flag for rnd=%s\n",
 464                               mpfr_print_rnd_mode(rnd));
 465                       printf ("expected %d, got %d\n", cmp, inexact);
 466                       printf ("x="); mpfr_print_binary (x); puts ("");
 467                       printf ("u="); mpfr_print_binary (u); puts ("");
 468                       printf ("y=  "); mpfr_print_binary (y); puts ("");
 469                       printf ("x-u="); mpfr_print_binary (z); puts ("");
 470                       exit (1);
 471                     }
 472                 }
 473             }
 474         }
 475     }
 476
 477   mpfr_clear (x);
 478   mpfr_clear (y);
 479   mpfr_clear (z);
 480   mpfr_clear (u);
 481 }
 482
 483 /* Bug found by Jakub Jelinek
 484  * http://bugzilla.redhat.com/643657
 485  * https://gforge.inria.fr/tracker/index.php?func=detail&aid=11301
 486  * The consequence can be either an assertion failure (i = 2 in the
 487  * testcase below, in debug mode) or an incorrectly rounded value.
 488  */
 489 static void
 490 bug20101017 (void)
 491 {
 492   mpfr_t a, b, c;
 493   int inex;
 494   int i;
 495
 496   mpfr_init2 (a, GMP_NUMB_BITS * 2);
 497   mpfr_init2 (b, GMP_NUMB_BITS);
 498   mpfr_init2 (c, GMP_NUMB_BITS);
 499
 500   /* a = 2^(2N) + k.2^(2N-1) + 2^N and b = 1
 501      with N = GMP_NUMB_BITS and k = 0 or 1.
 502      c = a - b should round to the same value as a. */
 503
 504   for (i = 2; i <= 3; i++)
 505     {
 506       mpfr_set_ui_2exp (a, i, GMP_NUMB_BITS - 1, MPFR_RNDN);
 507       mpfr_add_ui (a, a, 1, MPFR_RNDN);
 508       mpfr_mul_2ui (a, a, GMP_NUMB_BITS, MPFR_RNDN);
 509       mpfr_set_ui (b, 1, MPFR_RNDN);
 510       inex = mpfr_sub (c, a, b, MPFR_RNDN);
 511       mpfr_set (b, a, MPFR_RNDN);
 512       if (! mpfr_equal_p (c, b))
 513         {
 514           printf ("Error in bug20101017 for i = %d.\n", i);
 515           printf ("Expected ");
 516           mpfr_out_str (stdout, 16, 0, b, MPFR_RNDN);
 517           putchar ('\n');
 518           printf ("Got      ");
 519           mpfr_out_str (stdout, 16, 0, c, MPFR_RNDN);
 520           putchar ('\n');
 521           exit (1);
 522         }
 523       if (inex >= 0)
 524         {
 525           printf ("Error in bug20101017 for i = %d: bad inex value.\n", i);
 526           printf ("Expected negative, got %d.\n", inex);
 527           exit (1);
 528         }
 529     }
 530
 531   mpfr_set_prec (a, 64);
 532   mpfr_set_prec (b, 129);
 533   mpfr_set_prec (c, 2);
 534   mpfr_set_str_binary (b, "0.100000000000000000000000000000000000000000000000000000000000000010000000000000000000000000000000000000000000000000000000000000001E65");
 535   mpfr_set_str_binary (c, "0.10E1");
 536   inex = mpfr_sub (a, b, c, MPFR_RNDN);
 537   if (mpfr_cmp_ui_2exp (a, 1, 64) != 0 || inex >= 0)
 538     {
 539       printf ("Error in mpfr_sub for b-c for b=2^64+1+2^(-64), c=1\n");
 540       printf ("Expected result 2^64 with inex < 0\n");
 541       printf ("Got "); mpfr_print_binary (a);
 542       printf (" with inex=%d\n", inex);
 543       exit (1);
 544     }
 545
 546   mpfr_clears (a, b, c, (mpfr_ptr) 0);
 547 }
 548
 549 /* hard test of rounding */
 550 static void
 551 check_rounding (void)
 552 {
 553   mpfr_t a, b, c, res;
 554   mpfr_prec_t p;
 555   long k, l;
 556   int i;
 557
 558 #define MAXKL (2 * GMP_NUMB_BITS)
 559   for (p = MPFR_PREC_MIN; p <= GMP_NUMB_BITS; p++)
 560     {
 561       mpfr_init2 (a, p);
 562       mpfr_init2 (res, p);
 563       mpfr_init2 (b, p + 1 + MAXKL);
 564       mpfr_init2 (c, MPFR_PREC_MIN);
 565
 566       /* b = 2^p + 1 + 2^(-k), c = 2^(-l) */
 567       for (k = 0; k <= MAXKL; k++)
 568         for (l = 0; l <= MAXKL; l++)
 569           {
 570             mpfr_set_ui_2exp (b, 1, p, MPFR_RNDN);
 571             mpfr_add_ui (b, b, 1, MPFR_RNDN);
 572             mpfr_mul_2ui (b, b, k, MPFR_RNDN);
 573             mpfr_add_ui (b, b, 1, MPFR_RNDN);
 574             mpfr_div_2ui (b, b, k, MPFR_RNDN);
 575             mpfr_set_ui_2exp (c, 1, -l, MPFR_RNDN);
 576             i = mpfr_sub (a, b, c, MPFR_RNDN);
 577             /* b - c = 2^p + 1 + 2^(-k) - 2^(-l), should be rounded to
 578                2^p for l <= k, and 2^p+2 for l < k */
 579             if (l <= k)
 580               {
 581                 if (mpfr_cmp_ui_2exp (a, 1, p) != 0)
 582                   {
 583                     printf ("Wrong result in check_rounding\n");
 584                     printf ("p=%lu k=%ld l=%ld\n", (unsigned long) p, k, l);
 585                     printf ("b="); mpfr_print_binary (b); puts ("");
 586                     printf ("c="); mpfr_print_binary (c); puts ("");
 587                     printf ("Expected 2^%lu\n", (unsigned long) p);
 588                     printf ("Got      "); mpfr_print_binary (a); puts ("");
 589                     exit (1);
 590                   }
 591                 if (i >= 0)
 592                   {
 593                     printf ("Wrong ternary value in check_rounding\n");
 594                     printf ("p=%lu k=%ld l=%ld\n", (unsigned long) p, k, l);
 595                     printf ("b="); mpfr_print_binary (b); puts ("");
 596                     printf ("c="); mpfr_print_binary (c); puts ("");
 597                     printf ("a="); mpfr_print_binary (a); puts ("");
 598                     printf ("Expected < 0, got %d\n", i);
 599                     exit (1);
 600                   }
 601               }
 602             else /* l < k */
 603               {
 604                 mpfr_set_ui_2exp (res, 1, p, MPFR_RNDN);
 605                 mpfr_add_ui (res, res, 2, MPFR_RNDN);
 606                 if (mpfr_cmp (a, res) != 0)
 607                   {
 608                     printf ("Wrong result in check_rounding\n");
 609                     printf ("b="); mpfr_print_binary (b); puts ("");
 610                     printf ("c="); mpfr_print_binary (c); puts ("");
 611                     printf ("Expected "); mpfr_print_binary (res); puts ("");
 612                     printf ("Got      "); mpfr_print_binary (a); puts ("");
 613                     exit (1);
 614                   }
 615                 if (i <= 0)
 616                   {
 617                     printf ("Wrong ternary value in check_rounding\n");
 618                     printf ("b="); mpfr_print_binary (b); puts ("");
 619                     printf ("c="); mpfr_print_binary (c); puts ("");
 620                     printf ("Expected > 0, got %d\n", i);
 621                     exit (1);
 622                   }
 623               }
 624           }
 625
 626       mpfr_clear (a);
 627       mpfr_clear (res);
 628       mpfr_clear (b);
 629       mpfr_clear (c);
 630     }
 631 }
 632
 633 #define TEST_FUNCTION test_sub
 634 #define TWO_ARGS
 635 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
 636 #include "tgeneric.c"
 637
 638 int
 639 main (void)
 640 {
 641   mpfr_prec_t p;
 642   unsigned int i;
 643
 644   tests_start_mpfr ();
 645
 646   bug20101017 ();
 647   check_rounding ();
 648   check_diverse ();
 649   check_inexact ();
 650   bug_ddefour ();
 651   for (p=2; p<200; p++)
 652     for (i=0; i<50; i++)
 653       check_two_sum (p);
 654   test_generic (2, 800, 100);
 655
 656   tests_end_mpfr ();
 657   return 0;
 658 }