/* Test file for mpfr_sub.

Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
Contributed by the AriC and Caramel projects, INRIA.

This file is part of the GNU MPFR Library.

The GNU MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

The GNU MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */

#include <stdio.h>
#include <stdlib.h>

#include "mpfr-test.h"

#ifdef CHECK_EXTERNAL
static int
test_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
  int res;
  int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c);
  if (ok)
    {
      mpfr_print_raw (b);
      printf (" ");
      mpfr_print_raw (c);
    }
  res = mpfr_sub (a, b, c, rnd_mode);
  if (ok)
    {
      printf (" ");
      mpfr_print_raw (a);
      printf ("\n");
    }
  return res;
}
#else
#define test_sub mpfr_sub
#endif

static void
check_diverse (void)
{
  mpfr_t x, y, z;
  int inexact;

  mpfr_init (x);
  mpfr_init (y);
  mpfr_init (z);

  /* check corner case cancel=0, but add_exp=1 */
  mpfr_set_prec (x, 2);
  mpfr_set_prec (y, 4);
  mpfr_set_prec (z, 2);
  mpfr_setmax (y, __gmpfr_emax);
  mpfr_set_str_binary (z, "0.1E-10"); /* tiny */
  test_sub (x, y, z, MPFR_RNDN); /* should round to 2^emax, i.e. overflow */
  if (!mpfr_inf_p (x) || mpfr_sgn (x) < 0)
    {
      printf ("Error in mpfr_sub(a,b,c,RNDN) for b=maxfloat, prec(a)<prec(b), c tiny\n");
      exit (1);
    }

  /* other coverage test */
  mpfr_set_prec (x, 2);
  mpfr_set_prec (y, 2);
  mpfr_set_prec (z, 2);
  mpfr_set_ui (y, 1, MPFR_RNDN);
  mpfr_set_si (z, -2, MPFR_RNDN);
  test_sub (x, y, z, MPFR_RNDD);
  if (mpfr_cmp_ui (x, 3))
    {
      printf ("Error in mpfr_sub(1,-2,RNDD)\n");
      exit (1);
    }

  mpfr_set_prec (x, 288);
  mpfr_set_prec (y, 288);
  mpfr_set_prec (z, 288);
  mpfr_set_str_binary (y, "0.111000110011000001000111101010111011110011101001101111111110000011100101000001001010110010101010011001010100000001110011110001010101101010001011101110100100001011110100110000101101100011010001001011011010101010000010001101001000110010010111111011110001111101001000101101001100101100101000E80");
  mpfr_set_str_binary (z, "0.100001111111101001011010001100110010100111001110000110011101001011010100001000000100111011010110110010000000000010101101011000010000110001110010100001100101011100100100001011000100011110000001010101000100011101001000010111100000111000111011001000100100011000100000010010111000000100100111E-258");
  inexact = test_sub (x, y, z, MPFR_RNDN);
  if (inexact <= 0)
    {
      printf ("Wrong inexact flag for prec=288\n");
      exit (1);
    }

  mpfr_set_prec (x, 32);
  mpfr_set_prec (y, 63);
  mpfr_set_prec (z, 63);
  mpfr_set_str_binary (x, "0.101101111011011100100100100111E31");
  mpfr_set_str_binary (y, "0.111110010010100100110101101010001001100101110001000101110111111E-1");
  test_sub (z, x, y, MPFR_RNDN);
  mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
  if (mpfr_cmp (z, y))
    {
      printf ("Error in mpfr_sub (5)\n");
      printf ("expected "); mpfr_print_binary (y); puts ("");
      printf ("got      "); mpfr_print_binary (z); puts ("");
      exit (1);
    }

  mpfr_set_prec (y, 63);
  mpfr_set_prec (z, 63);
  mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
  mpfr_sub_ui (z, y, 1541116494, MPFR_RNDN);
  mpfr_set_str_binary (y, "-0.11111001001010010011010110101E-1");
  if (mpfr_cmp (z, y))
    {
      printf ("Error in mpfr_sub (7)\n");
      printf ("expected "); mpfr_print_binary (y); puts ("");
      printf ("got      "); mpfr_print_binary (z); puts ("");
      exit (1);
    }

  mpfr_set_prec (y, 63);
  mpfr_set_prec (z, 63);
  mpfr_set_str_binary (y, "0.1011011110110111001001001001101100000110110101101100101001011E31");
  mpfr_sub_ui (z, y, 1541116494, MPFR_RNDN);
  mpfr_set_str_binary (y, "-0.11111001001010010011010110101E-1");
  if (mpfr_cmp (z, y))
    {
      printf ("Error in mpfr_sub (6)\n");
      printf ("expected "); mpfr_print_binary (y); puts ("");
      printf ("got      "); mpfr_print_binary (z); puts ("");
      exit (1);
    }

  mpfr_set_prec (x, 32);
  mpfr_set_prec (y, 32);
  mpfr_set_str_binary (x, "0.10110111101001110100100101111000E0");
  mpfr_set_str_binary (y, "0.10001100100101000100110111000100E0");
  if ((inexact = test_sub (x, x, y, MPFR_RNDN)))
    {
      printf ("Wrong inexact flag (2): got %d instead of 0\n", inexact);
      exit (1);
    }

  mpfr_set_prec (x, 32);
  mpfr_set_prec (y, 32);
  mpfr_set_str_binary (x, "0.11111000110111011000100111011010E0");
  mpfr_set_str_binary (y, "0.10011111101111000100001000000000E-3");
  if ((inexact = test_sub (x, x, y, MPFR_RNDN)))
    {
      printf ("Wrong inexact flag (1): got %d instead of 0\n", inexact);
      exit (1);
    }

  mpfr_set_prec (x, 33);
  mpfr_set_prec (y, 33);
  mpfr_set_ui (x, 1, MPFR_RNDN);
  mpfr_set_str_binary (y, "-0.1E-32");
  mpfr_add (x, x, y, MPFR_RNDN);
  mpfr_set_str_binary (y, "0.111111111111111111111111111111111E0");
  if (mpfr_cmp (x, y))
    {
      printf ("Error in mpfr_sub (1 - 1E-33) with prec=33\n");
      printf ("Expected "); mpfr_print_binary (y); puts ("");
      printf ("got      "); mpfr_print_binary (x); puts ("");
      exit (1);
    }

  mpfr_set_prec (x, 32);
  mpfr_set_prec (y, 33);
  mpfr_set_ui (x, 1, MPFR_RNDN);
  mpfr_set_str_binary (y, "-0.1E-32");
  mpfr_add (x, x, y, MPFR_RNDN);
  if (mpfr_cmp_ui (x, 1))
    {
      printf ("Error in mpfr_sub (1 - 1E-33) with prec=32\n");
      printf ("Expected 1.0, got "); mpfr_print_binary (x); puts ("");
      exit (1);
    }

  mpfr_set_prec (x, 65);
  mpfr_set_prec (y, 65);
  mpfr_set_prec (z, 64);
  mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
  mpfr_set_str_binary (y, "0.1110111011110001110111011111111111101000011001011100101100101100");
  test_sub (z, x, y, MPFR_RNDZ);
  if (mpfr_cmp_ui (z, 1))
    {
      printf ("Error in mpfr_sub (1)\n");
      exit (1);
    }
  test_sub (z, x, y, MPFR_RNDU);
  mpfr_set_str_binary (x, "1.000000000000000000000000000000000000000000000000000000000000001");
  if (mpfr_cmp (z, x))
    {
      printf ("Error in mpfr_sub (2)\n");
      printf ("Expected "); mpfr_print_binary (x); puts ("");
      printf ("Got      "); mpfr_print_binary (z); puts ("");
      exit (1);
    }
  mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
  test_sub (z, x, y, MPFR_RNDN);
  if (mpfr_cmp_ui (z, 1))
    {
      printf ("Error in mpfr_sub (3)\n");
      exit (1);
    }
  inexact = test_sub (z, x, y, MPFR_RNDA);
  mpfr_set_str_binary (x, "1.000000000000000000000000000000000000000000000000000000000000001");
  if (mpfr_cmp (z, x) || inexact <= 0)
    {
      printf ("Error in mpfr_sub (4)\n");
      exit (1);
    }
  mpfr_set_prec (x, 66);
  mpfr_set_str_binary (x, "1.11101110111100011101110111111111111010000110010111001011001010111");
  test_sub (z, x, y, MPFR_RNDN);
  if (mpfr_cmp_ui (z, 1))
    {
      printf ("Error in mpfr_sub (5)\n");
      exit (1);
    }

  /* check in-place operations */
  mpfr_set_ui (x, 1, MPFR_RNDN);
  test_sub (x, x, x, MPFR_RNDN);
  if (mpfr_cmp_ui(x, 0))
    {
      printf ("Error for mpfr_sub (x, x, x, MPFR_RNDN) with x=1.0\n");
      exit (1);
    }

  mpfr_set_prec (x, 53);
  mpfr_set_prec (y, 53);
  mpfr_set_prec (z, 53);
  mpfr_set_str1 (x, "1.229318102e+09");
  mpfr_set_str1 (y, "2.32221184180698677665e+05");
  test_sub (z, x, y, MPFR_RNDN);
  if (mpfr_cmp_str1 (z, "1229085880.815819263458251953125"))
    {
      printf ("Error in mpfr_sub (1.22e9 - 2.32e5)\n");
      printf ("expected 1229085880.815819263458251953125, got ");
      mpfr_out_str(stdout, 10, 0, z, MPFR_RNDN);
      putchar('\n');
      exit (1);
    }

  mpfr_set_prec (x, 112);
  mpfr_set_prec (y, 98);
  mpfr_set_prec (z, 54);
  mpfr_set_str_binary (x, "0.11111100100000000011000011100000101101010001000111E-401");
  mpfr_set_str_binary (y, "0.10110000100100000101101100011111111011101000111000101E-464");
  test_sub (z, x, y, MPFR_RNDN);
  if (mpfr_cmp (z, x)) {
    printf ("mpfr_sub(z, x, y) failed for prec(x)=112, prec(y)=98\n");
    printf ("expected "); mpfr_print_binary (x); puts ("");
    printf ("got      "); mpfr_print_binary (z); puts ("");
    exit (1);
  }

  mpfr_set_prec (x, 33);
  mpfr_set_ui (x, 1, MPFR_RNDN);
  mpfr_div_2exp (x, x, 32, MPFR_RNDN);
  mpfr_sub_ui (x, x, 1, MPFR_RNDN);

  mpfr_set_prec (x, 5);
  mpfr_set_prec (y, 5);
  mpfr_set_str_binary (x, "1e-12");
  mpfr_set_ui (y, 1, MPFR_RNDN);
  test_sub (x, y, x, MPFR_RNDD);
  mpfr_set_str_binary (y, "0.11111");
  if (mpfr_cmp (x, y))
    {
      printf ("Error in mpfr_sub (x, y, x, MPFR_RNDD) for x=2^(-12), y=1\n");
      exit (1);
    }

  mpfr_set_prec (x, 24);
  mpfr_set_prec (y, 24);
  mpfr_set_str_binary (x, "-0.100010000000000000000000E19");
  mpfr_set_str_binary (y, "0.100000000000000000000100E15");
  mpfr_add (x, x, y, MPFR_RNDD);
  mpfr_set_str_binary (y, "-0.1E19");
  if (mpfr_cmp (x, y))
    {
      printf ("Error in mpfr_add (2)\n");
      exit (1);
    }

  mpfr_set_prec (x, 2);
  mpfr_set_prec (y, 10);
  mpfr_set_prec (z, 10);
  mpfr_set_ui (y, 0, MPFR_RNDN);
  mpfr_set_str_binary (z, "0.10001");
  if (test_sub (x, y, z, MPFR_RNDN) <= 0)
    {
      printf ("Wrong inexact flag in x=mpfr_sub(0,z) for prec(z)>prec(x)\n");
      exit (1);
    }
  if (test_sub (x, z, y, MPFR_RNDN) >= 0)
    {
      printf ("Wrong inexact flag in x=mpfr_sub(z,0) for prec(z)>prec(x)\n");
      exit (1);
    }

  mpfr_clear (x);
  mpfr_clear (y);
  mpfr_clear (z);
}

static void
bug_ddefour(void)
{
    mpfr_t ex, ex1, ex2, ex3, tot, tot1;

    mpfr_init2(ex, 53);
    mpfr_init2(ex1, 53);
    mpfr_init2(ex2, 53);
    mpfr_init2(ex3, 53);
    mpfr_init2(tot, 150);
    mpfr_init2(tot1, 150);

    mpfr_set_ui( ex, 1, MPFR_RNDN);
    mpfr_mul_2exp( ex, ex, 906, MPFR_RNDN);
    mpfr_log( tot, ex, MPFR_RNDN);
    mpfr_set( ex1, tot, MPFR_RNDN); /* ex1 = high(tot) */
    test_sub( ex2, tot, ex1, MPFR_RNDN); /* ex2 = high(tot - ex1) */
    test_sub( tot1, tot, ex1, MPFR_RNDN); /* tot1 = tot - ex1 */
    mpfr_set( ex3, tot1, MPFR_RNDN); /* ex3 = high(tot - ex1) */

    if (mpfr_cmp(ex2, ex3))
      {
        printf ("Error in ddefour test.\n");
        printf ("ex2="); mpfr_print_binary (ex2); puts ("");
        printf ("ex3="); mpfr_print_binary (ex3); puts ("");
        exit (1);
      }

    mpfr_clear (ex);
    mpfr_clear (ex1);
    mpfr_clear (ex2);
    mpfr_clear (ex3);
    mpfr_clear (tot);
    mpfr_clear (tot1);
}

/* if u = o(x-y), v = o(u-x), w = o(v+y), then x-y = u-w */
static void
check_two_sum (mpfr_prec_t p)
{
  mpfr_t x, y, u, v, w;
  mpfr_rnd_t rnd;
  int inexact;

  mpfr_init2 (x, p);
  mpfr_init2 (y, p);
  mpfr_init2 (u, p);
  mpfr_init2 (v, p);
  mpfr_init2 (w, p);
  mpfr_urandomb (x, RANDS);
  mpfr_urandomb (y, RANDS);
  if (mpfr_cmpabs (x, y) < 0)
    mpfr_swap (x, y);
  rnd = MPFR_RNDN;
  inexact = test_sub (u, x, y, rnd);
  test_sub (v, u, x, rnd);
  mpfr_add (w, v, y, rnd);
  /* as u = (x-y) - w, we should have inexact and w of opposite signs */
  if (((inexact == 0) && mpfr_cmp_ui (w, 0)) ||
      ((inexact > 0) && (mpfr_cmp_ui (w, 0) <= 0)) ||
      ((inexact < 0) && (mpfr_cmp_ui (w, 0) >= 0)))
    {
      printf ("Wrong inexact flag for prec=%u, rnd=%s\n", (unsigned)p,
               mpfr_print_rnd_mode (rnd));
      printf ("x="); mpfr_print_binary(x); puts ("");
      printf ("y="); mpfr_print_binary(y); puts ("");
      printf ("u="); mpfr_print_binary(u); puts ("");
      printf ("v="); mpfr_print_binary(v); puts ("");
      printf ("w="); mpfr_print_binary(w); puts ("");
      printf ("inexact = %d\n", inexact);
      exit (1);
    }
  mpfr_clear (x);
  mpfr_clear (y);
  mpfr_clear (u);
  mpfr_clear (v);
  mpfr_clear (w);
}

#define MAX_PREC 200

static void
check_inexact (void)
{
  mpfr_t x, y, z, u;
  mpfr_prec_t px, py, pu, pz;
  int inexact, cmp;
  mpfr_rnd_t rnd;

  mpfr_init (x);
  mpfr_init (y);
  mpfr_init (z);
  mpfr_init (u);

  mpfr_set_prec (x, 2);
  mpfr_set_ui (x, 6, MPFR_RNDN);
  mpfr_div_2exp (x, x, 4, MPFR_RNDN); /* x = 6/16 */
  mpfr_set_prec (y, 2);
  mpfr_set_si (y, -1, MPFR_RNDN);
  mpfr_div_2exp (y, y, 4, MPFR_RNDN); /* y = -1/16 */
  inexact = test_sub (y, y, x, MPFR_RNDN); /* y = round(-7/16) = -1/2 */
  if (inexact >= 0)
    {
      printf ("Error: wrong inexact flag for -1/16 - (6/16)\n");
      exit (1);
    }

  for (px=2; px<MAX_PREC; px++)
    {
      mpfr_set_prec (x, px);
      do
        {
          mpfr_urandomb (x, RANDS);
        }
      while (mpfr_cmp_ui (x, 0) == 0);
      for (pu=2; pu<MAX_PREC; pu++)
        {
          mpfr_set_prec (u, pu);
          do
            {
              mpfr_urandomb (u, RANDS);
            }
          while (mpfr_cmp_ui (u, 0) == 0);
          {
              py = 2 + (randlimb () % (MAX_PREC - 2));
              mpfr_set_prec (y, py);
              /* warning: MPFR_EXP is undefined for 0 */
              pz =  (mpfr_cmpabs (x, u) >= 0) ? MPFR_EXP(x) - MPFR_EXP(u)
                : MPFR_EXP(u) - MPFR_EXP(x);
              pz = pz + MAX(MPFR_PREC(x), MPFR_PREC(u));
              mpfr_set_prec (z, pz);
              rnd = RND_RAND ();
              if (test_sub (z, x, u, rnd))
                {
                  printf ("z <- x - u should be exact\n");
                  exit (1);
                }
                {
                  rnd = RND_RAND ();
                  inexact = test_sub (y, x, u, rnd);
                  cmp = mpfr_cmp (y, z);
                  if (((inexact == 0) && (cmp != 0)) ||
                      ((inexact > 0) && (cmp <= 0)) ||
                      ((inexact < 0) && (cmp >= 0)))
                    {
                      printf ("Wrong inexact flag for rnd=%s\n",
                              mpfr_print_rnd_mode(rnd));
                      printf ("expected %d, got %d\n", cmp, inexact);
                      printf ("x="); mpfr_print_binary (x); puts ("");
                      printf ("u="); mpfr_print_binary (u); puts ("");
                      printf ("y=  "); mpfr_print_binary (y); puts ("");
                      printf ("x-u="); mpfr_print_binary (z); puts ("");
                      exit (1);
                    }
                }
            }
        }
    }

  mpfr_clear (x);
  mpfr_clear (y);
  mpfr_clear (z);
  mpfr_clear (u);
}

/* Bug found by Jakub Jelinek
 * http://bugzilla.redhat.com/643657
 * https://gforge.inria.fr/tracker/index.php?func=detail&aid=11301
 * The consequence can be either an assertion failure (i = 2 in the
 * testcase below, in debug mode) or an incorrectly rounded value.
 */
static void
bug20101017 (void)
{
  mpfr_t a, b, c;
  int inex;
  int i;

  mpfr_init2 (a, GMP_NUMB_BITS * 2);
  mpfr_init2 (b, GMP_NUMB_BITS);
  mpfr_init2 (c, GMP_NUMB_BITS);

  /* a = 2^(2N) + k.2^(2N-1) + 2^N and b = 1
     with N = GMP_NUMB_BITS and k = 0 or 1.
     c = a - b should round to the same value as a. */

  for (i = 2; i <= 3; i++)
    {
      mpfr_set_ui_2exp (a, i, GMP_NUMB_BITS - 1, MPFR_RNDN);
      mpfr_add_ui (a, a, 1, MPFR_RNDN);
      mpfr_mul_2ui (a, a, GMP_NUMB_BITS, MPFR_RNDN);
      mpfr_set_ui (b, 1, MPFR_RNDN);
      inex = mpfr_sub (c, a, b, MPFR_RNDN);
      mpfr_set (b, a, MPFR_RNDN);
      if (! mpfr_equal_p (c, b))
        {
          printf ("Error in bug20101017 for i = %d.\n", i);
          printf ("Expected ");
          mpfr_out_str (stdout, 16, 0, b, MPFR_RNDN);
          putchar ('\n');
          printf ("Got      ");
          mpfr_out_str (stdout, 16, 0, c, MPFR_RNDN);
          putchar ('\n');
          exit (1);
        }
      if (inex >= 0)
        {
          printf ("Error in bug20101017 for i = %d: bad inex value.\n", i);
          printf ("Expected negative, got %d.\n", inex);
          exit (1);
        }
    }

  mpfr_set_prec (a, 64);
  mpfr_set_prec (b, 129);
  mpfr_set_prec (c, 2);
  mpfr_set_str_binary (b, "0.100000000000000000000000000000000000000000000000000000000000000010000000000000000000000000000000000000000000000000000000000000001E65");
  mpfr_set_str_binary (c, "0.10E1");
  inex = mpfr_sub (a, b, c, MPFR_RNDN);
  if (mpfr_cmp_ui_2exp (a, 1, 64) != 0 || inex >= 0)
    {
      printf ("Error in mpfr_sub for b-c for b=2^64+1+2^(-64), c=1\n");
      printf ("Expected result 2^64 with inex < 0\n");
      printf ("Got "); mpfr_print_binary (a);
      printf (" with inex=%d\n", inex);
      exit (1);
    }

  mpfr_clears (a, b, c, (mpfr_ptr) 0);
}

/* hard test of rounding */
static void
check_rounding (void)
{
  mpfr_t a, b, c, res;
  mpfr_prec_t p;
  long k, l;
  int i;

#define MAXKL (2 * GMP_NUMB_BITS)
  for (p = MPFR_PREC_MIN; p <= GMP_NUMB_BITS; p++)
    {
      mpfr_init2 (a, p);
      mpfr_init2 (res, p);
      mpfr_init2 (b, p + 1 + MAXKL);
      mpfr_init2 (c, MPFR_PREC_MIN);

      /* b = 2^p + 1 + 2^(-k), c = 2^(-l) */
      for (k = 0; k <= MAXKL; k++)
        for (l = 0; l <= MAXKL; l++)
          {
            mpfr_set_ui_2exp (b, 1, p, MPFR_RNDN);
            mpfr_add_ui (b, b, 1, MPFR_RNDN);
            mpfr_mul_2ui (b, b, k, MPFR_RNDN);
            mpfr_add_ui (b, b, 1, MPFR_RNDN);
            mpfr_div_2ui (b, b, k, MPFR_RNDN);
            mpfr_set_ui_2exp (c, 1, -l, MPFR_RNDN);
            i = mpfr_sub (a, b, c, MPFR_RNDN);
            /* b - c = 2^p + 1 + 2^(-k) - 2^(-l), should be rounded to
               2^p for l <= k, and 2^p+2 for l < k */
            if (l <= k)
              {
                if (mpfr_cmp_ui_2exp (a, 1, p) != 0)
                  {
                    printf ("Wrong result in check_rounding\n");
                    printf ("p=%lu k=%ld l=%ld\n", (unsigned long) p, k, l);
                    printf ("b="); mpfr_print_binary (b); puts ("");
                    printf ("c="); mpfr_print_binary (c); puts ("");
                    printf ("Expected 2^%lu\n", (unsigned long) p);
                    printf ("Got      "); mpfr_print_binary (a); puts ("");
                    exit (1);
                  }
                if (i >= 0)
                  {
                    printf ("Wrong ternary value in check_rounding\n");
                    printf ("p=%lu k=%ld l=%ld\n", (unsigned long) p, k, l);
                    printf ("b="); mpfr_print_binary (b); puts ("");
                    printf ("c="); mpfr_print_binary (c); puts ("");
                    printf ("a="); mpfr_print_binary (a); puts ("");
                    printf ("Expected < 0, got %d\n", i);
                    exit (1);
                  }
              }
            else /* l < k */
              {
                mpfr_set_ui_2exp (res, 1, p, MPFR_RNDN);
                mpfr_add_ui (res, res, 2, MPFR_RNDN);
                if (mpfr_cmp (a, res) != 0)
                  {
                    printf ("Wrong result in check_rounding\n");
                    printf ("b="); mpfr_print_binary (b); puts ("");
                    printf ("c="); mpfr_print_binary (c); puts ("");
                    printf ("Expected "); mpfr_print_binary (res); puts ("");
                    printf ("Got      "); mpfr_print_binary (a); puts ("");
                    exit (1);
                  }
                if (i <= 0)
                  {
                    printf ("Wrong ternary value in check_rounding\n");
                    printf ("b="); mpfr_print_binary (b); puts ("");
                    printf ("c="); mpfr_print_binary (c); puts ("");
                    printf ("Expected > 0, got %d\n", i);
                    exit (1);
                  }
              }
          }

      mpfr_clear (a);
      mpfr_clear (res);
      mpfr_clear (b);
      mpfr_clear (c);
    }
}

#define TEST_FUNCTION test_sub
#define TWO_ARGS
#define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
#include "tgeneric.c"

int
main (void)
{
  mpfr_prec_t p;
  unsigned int i;

  tests_start_mpfr ();

  bug20101017 ();
  check_rounding ();
  check_diverse ();
  check_inexact ();
  bug_ddefour ();
  for (p=2; p<200; p++)
    for (i=0; i<50; i++)
      check_two_sum (p);
  test_generic (2, 800, 100);

  tests_end_mpfr ();
  return 0;
}