static void
check_round_p (void)
{
  mp_limb_t buf[MAX_LIMB_SIZE];
  mp_size_t n, i;
  mp_prec_t p;
  mp_exp_t err;
  int r1, r2;

  for (n = 2 ; n <= MAX_LIMB_SIZE ; n++)
    {
      /* avoid mpn_random which leaks memory */
      for (i = 0; i < n; i++)
        buf[i] = randlimb ();
      p = (mp_prec_t) randlimb() % ((n-1) * BITS_PER_MP_LIMB) + MPFR_PREC_MIN;
      err = p + randlimb () % BITS_PER_MP_LIMB;
      r1 = mpfr_round_p (buf, n, err, p);
      r2 = mpfr_can_round_raw (buf, n, MPFR_SIGN_POS, err,
                               GMP_RNDN, GMP_RNDZ, p);
      if (r1 != r2)
        {
          printf ("mpfr_round_p(%d) != mpfr_can_round(%d)!\n"
                  "bn = %ld, err0 = %ld, prec = %lu\nbp = ",
                  r1, r2, n, err, (unsigned long) p);
          gmp_printf ("%NX\n", buf, n);
          exit (1);
        }
    }
}

/* pos is 512 times the proportion of negative numbers.
   If pos=256, half of the numbers are negative.
   If pos=0, all generated numbers are positive.
*/
void
tests_default_random (mpfr_ptr x, int pos, mpfr_exp_t emin, mpfr_exp_t emax)
{
  MPFR_ASSERTN (emin <= emax);
  MPFR_ASSERTN (emin >= MPFR_EMIN_MIN);
  MPFR_ASSERTN (emax <= MPFR_EMAX_MAX);
  /* but it isn't required that emin and emax are in the current
     exponent range (see below), so that underflow/overflow checks
     can be done on 64-bit machines. */

  mpfr_urandomb (x, RANDS);
  if (MPFR_IS_PURE_FP (x) && (emin >= 1 || (randlimb () & 1)))
    {
      mpfr_exp_t e;
      e = MPFR_GET_EXP (x) +
        (emin + (long) (randlimb () % (emax - emin + 1)));
      /* Note: There should be no overflow here because both terms are
         between MPFR_EMIN_MIN and MPFR_EMAX_MAX, but the sum e isn't
         necessarily between MPFR_EMIN_MIN and MPFR_EMAX_MAX. */
      if (mpfr_set_exp (x, e))
        {
          /* The random number doesn't fit in the current exponent range.
             In this case, test the function in the extended exponent range,
             which should be restored by the caller. */
          mpfr_set_emin (MPFR_EMIN_MIN);
          mpfr_set_emax (MPFR_EMAX_MAX);
          mpfr_set_exp (x, e);
        }
    }
  if (randlimb () % 512 < pos)
    mpfr_neg (x, x, MPFR_RNDN);
}

static void
cmp_tests (void)
{
  mpfr_t x, y;
  long i;

  mpfr_inits (x, y, (mpfr_ptr) 0);
  for (i = 0; i < 80000; i++)
    {
      mpfr_prec_t precx, precy;
      int signx, signy, cmp;
      unsigned int cmpbool = 0;

      precx = (randlimb () % 17) * 11 + MPFR_PREC_MIN;
      precy = (randlimb () % 17) * 11 + MPFR_PREC_MIN;
      mpfr_set_prec (x, precx);
      mpfr_set_prec (y, precy);
      mpfr_urandomb (x, RANDS);
      mpfr_urandomb (y, RANDS);
      signx = randlimb () & 1;
      signy = randlimb () % 256 ? signx : 1 - signx;
      /* signy = signx most of the time (most interesting case) */
      if (signx)
        mpfr_neg (x, x, MPFR_RNDN);
      if (signy)
        mpfr_neg (y, y, MPFR_RNDN);
      if (i <= 1)
        mpfr_set_nan (x);
      if (i == 0 || i == 2)
        mpfr_set_nan (y);
      if (mpfr_greater_p (x, y))
        cmpbool |= 0x01;
      if (mpfr_greaterequal_p (x, y))
        cmpbool |= 0x02;
      if (mpfr_less_p (x, y))
        cmpbool |= 0x04;
      if (mpfr_lessequal_p (x, y))
        cmpbool |= 0x08;
      if (mpfr_lessgreater_p (x, y))
        cmpbool |= 0x10;
      if (mpfr_equal_p (x, y))
        cmpbool |= 0x20;
      if (mpfr_unordered_p (x, y))
        cmpbool |= 0x40;
      if ((i <= 2 && cmpbool != 0x40) ||
          (i > 2 && (cmp = mpfr_cmp (x, y),
                     (cmp == 0 && cmpbool != 0x2a) ||
                     (cmp < 0 && cmpbool != 0x1c) ||
                     (cmp > 0 && cmpbool != 0x13))))
        {
          printf ("Error in cmp_tests for\nx = ");
          mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
          printf (" and\ny = ");
          mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
          printf ("\n");
          exit (1);
        }
    }
  mpfr_clears (x, y, (mpfr_ptr) 0);
}

/* pos is 512 times the proportion of negative numbers.
   If pos=256, half of the numbers are negative.
   If pos=0, all generated numbers are positive.
*/
void
tests_default_random (mpfr_ptr x, int pos, mpfr_exp_t emin, mpfr_exp_t emax,
                      int always_scale)
{
  MPFR_ASSERTN (emin <= emax);
  MPFR_ASSERTN (emin >= MPFR_EMIN_MIN);
  MPFR_ASSERTN (emax <= MPFR_EMAX_MAX);
  /* but it isn't required that emin and emax are in the current
     exponent range (see below), so that underflow/overflow checks
     can be done on 64-bit machines without a manual change of the
     exponent range (well, this is a bit ugly...). */

  mpfr_urandomb (x, RANDS);
  if (MPFR_IS_PURE_FP (x) && (emin >= 1 || always_scale || (randlimb () & 1)))
    {
      mpfr_exp_t e;
      e = emin + (mpfr_exp_t) (randlimb () % (emax - emin + 1));
      /* Note: There should be no overflow here because both terms are
         between MPFR_EMIN_MIN and MPFR_EMAX_MAX. */
      MPFR_ASSERTD (e >= emin && e <= emax);
      if (mpfr_set_exp (x, e))
        {
          /* The random number doesn't fit in the current exponent range.
             In this case, test the function in the extended exponent range,
             which should be restored by the caller. */
          mpfr_set_emin (MPFR_EMIN_MIN);
          mpfr_set_emax (MPFR_EMAX_MAX);
          mpfr_set_exp (x, e);
        }
    }
  if (randlimb () % 512 < pos)
    mpfr_neg (x, x, MPFR_RNDN);
}

static void
eq_tests (void)
{
  mpfr_t x, y;
  long i;

  mpfr_inits (x, y, (mpfr_ptr) 0);
  for (i = 0; i < 20000; i++)
    {
      mpfr_prec_t precx;

      precx = (randlimb () % 17) * 11 + MPFR_PREC_MIN;
      mpfr_set_prec (x, precx);
      mpfr_set_prec (y, precx + (randlimb () % 64));
      mpfr_urandomb (x, RANDS);
      if (randlimb () & 1)
        mpfr_neg (x, x, MPFR_RNDN);
      mpfr_set (y, x, MPFR_RNDN);  /* exact -> x = y */
      if (mpfr_greater_p (x, y) || !mpfr_greaterequal_p (x, y) ||
          mpfr_less_p (x, y) || !mpfr_lessequal_p (x, y) ||
          mpfr_lessgreater_p (x, y) || !mpfr_equal_p (x, y) ||
          mpfr_unordered_p (x, y))
        {
          printf ("Error in eq_tests for x = ");
          mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
          printf ("\n");
          exit (1);
        }
    }
  mpfr_clears (x, y, (mpfr_ptr) 0);
}

static void
test_cmp_q (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
{
  mpfr_t x, z;
  mpq_t  y;
  mpfr_prec_t p;
  int res1, res2;
  int n;

  mpfr_init (x);
  mpfr_init2 (z, MPFR_PREC_MIN);
  mpq_init (y);

  /* check the erange flag when x is NaN */
  mpfr_set_nan (x);
  mpq_set_ui (y, 17, 1);
  mpfr_clear_erangeflag ();
  res1 = mpfr_cmp_q (x, y);
  if (res1 != 0 || mpfr_erangeflag_p () == 0)
    {
      printf ("Error for mpfr_cmp_q (NaN, 17)\n");
      printf ("Return value: expected 0, got %d\n", res1);
      printf ("Erange flag: expected set, got %d\n", mpfr_erangeflag_p ());
      exit (1);
    }

  for(p=pmin ; p < pmax ; p++)
    {
      mpfr_set_prec (x, p);
      for (n = 0 ; n < nmax ; n++)
        {
          mpfr_urandomb (x, RANDS);
          mpq_set_ui (y, randlimb (), randlimb() );
          if (!MPFR_IS_SINGULAR (x))
            {
              mpfr_sub_q (z, x, y, MPFR_RNDN);
              res1 = mpfr_sgn (z);
              res2 = mpfr_cmp_q (x, y);
              if (res1 != res2)
                {
                  printf("Error for mpfr_cmp_q: res=%d sub_z gives %d\n",
                         res2, res1);
                  exit (1);
                }
            }
        }
    }
  mpq_clear (y);
  mpfr_clear (x);
  mpfr_clear (z);
}

/* 1. compare expected string with the string BUFFER returned by
   mpfr_sprintf(buffer, fmt, x)
   2. then test mpfr_snprintf (buffer, p, fmt, x) with a random p. */
static int
check_sprintf (const char *expected, const char *fmt, mpfr_srcptr x)
{
    int n0, n1, p;
    char buffer[BUF_SIZE];

    /* test mpfr_sprintf */
    n0 = mpfr_sprintf (buffer, fmt, x);
    if (strcmp (buffer, expected) != 0)
    {
        printf ("Error in mpfr_sprintf (s, \"%s\", x);\n", fmt);
        printf ("expected: \"%s\"\ngot:      \"%s\"\n", expected, buffer);

        exit (1);
    }

    /* test mpfr_snprintf */
    p = (int) (randlimb () % n0);
    if (p == 0 && (randlimb () & 1) == 0)
    {
        n1 = mpfr_snprintf (NULL, 0, fmt, x);
    }
    else
    {
        buffer[p] = 17;
        n1 = mpfr_snprintf (buffer, p, fmt, x);
        if (buffer[p] != 17)
        {
            printf ("Buffer overflow in mpfr_snprintf for p = %d!\n", p);
            exit (1);
        }
    }
    if (n0 != n1)
    {
        printf ("Error in mpfr_snprintf (s, %d, \"%s\", x) return value\n",
                p, fmt);
        printf ("expected: %d\ngot:      %d\n", n0, n1);
        exit (1);
    }
    if ((p > 1 && strncmp (expected, buffer, p-1) != 0)
            || (p == 1 && buffer[0] != '\0'))
    {
        char part_expected[BUF_SIZE];
        strncpy (part_expected, expected, p);
        part_expected[p-1] = '\0';
        printf ("Error in mpfr_vsnprintf (s, %d, \"%s\", ...);\n", p, fmt);
        printf ("expected: \"%s\"\ngot:      \"%s\"\n", part_expected, buffer);
        exit (1);
    }
    return n0;
}

static void
test_specialq (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N,
               int (*mpfr_func)(mpfr_ptr, mpfr_srcptr, mpq_srcptr, mpfr_rnd_t),
               void (*mpq_func)(mpq_ptr, mpq_srcptr, mpq_srcptr),
               const char *op)
{
  mpfr_t fra, frb, frq;
  mpq_t  q1, q2, qr;
  unsigned int n;
  mpfr_prec_t prec;

  for (prec = p0 ; prec < p1 ; prec++)
    {
      mpfr_inits2 (prec, fra, frb, frq, (mpfr_ptr) 0);
      mpq_init (q1); mpq_init(q2); mpq_init (qr);

      for( n = 0 ; n < N ; n++)
        {
          mpq_set_ui(q1, randlimb(), randlimb() );
          mpq_set_ui(q2, randlimb(), randlimb() );
          mpq_canonicalize (q1);
          mpq_canonicalize (q2);
          mpq_func (qr, q1, q2);
          mpfr_set_q (fra, q1, MPFR_RNDD);
          mpfr_func (fra, fra, q2, MPFR_RNDD);
          mpfr_set_q (frb, q1, MPFR_RNDU);
          mpfr_func (frb, frb, q2, MPFR_RNDU);
          mpfr_set_q (frq, qr, MPFR_RNDN);
          /* We should have fra <= qr <= frb */
          if ( (mpfr_cmp(fra, frq) > 0) || (mpfr_cmp (frq, frb) > 0))
            {
              printf("Range error for prec=%lu and %s",
                     (unsigned long) prec, op);
              printf ("\nq1="); mpq_out_str(stdout, 2, q1);
              printf ("\nq2="); mpq_out_str(stdout, 2, q2);
              printf ("\nfr_dn="); mpfr_print_binary (fra);
              printf ("\nfr_q ="); mpfr_print_binary (frq);
              printf ("\nfr_up="); mpfr_print_binary (frb);
              putchar('\n');
              exit (1);
            }
        }

      mpq_clear (q1); mpq_clear (q2); mpq_clear (qr);
      mpfr_clears (fra, frb, frq, (mpfr_ptr) 0);
    }
}

static void
consistency (void)
{
  mpfr_t x, y, z1, z2;
  int i;

  mpfr_inits (x, y, z1, z2, (mpfr_ptr) 0);

  for (i = 0; i < 10000; i++)
    {
      mpfr_rnd_t rnd;
      mpfr_prec_t px, py, pz, p;
      int inex1, inex2;

      rnd = RND_RAND ();
      px = (randlimb () % 256) + 2;
      py = (randlimb () % 128) + 2;
      pz = (randlimb () % 256) + 2;
      mpfr_set_prec (x, px);
      mpfr_set_prec (y, py);
      mpfr_set_prec (z1, pz);
      mpfr_set_prec (z2, pz);
      mpfr_urandomb (x, RANDS);
      do
        mpfr_urandomb (y, RANDS);
      while (mpfr_zero_p (y));
      inex1 = mpfr_div (z1, x, y, rnd);
      MPFR_ASSERTN (!MPFR_IS_NAN (z1));
      p = MAX (MAX (px, py), pz);
      if (mpfr_prec_round (x, p, MPFR_RNDN) != 0 ||
          mpfr_prec_round (y, p, MPFR_RNDN) != 0)
        {
          printf ("mpfr_prec_round error for i = %d\n", i);
          exit (1);
        }
      inex2 = mpfr_div (z2, x, y, rnd);
      MPFR_ASSERTN (!MPFR_IS_NAN (z2));
      if (inex1 != inex2 || mpfr_cmp (z1, z2) != 0)
        {
          printf ("Consistency error for i = %d\n", i);
          exit (1);
        }
    }

  mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
}

static void
test_sum (mpfr_prec_t f, unsigned long n)
{
  mpfr_t sum, real_sum, real_non_rounded;
  mpfr_t *tab;
  unsigned long i;
  int rnd_mode;

  /* Init */
  tab = (mpfr_t *) (*__gmp_allocate_func) (n * sizeof(mpfr_t));
  for (i = 0; i < n; i++)
    mpfr_init2 (tab[i], f);
  mpfr_inits2 (f, sum, real_sum, real_non_rounded, (mpfr_ptr) 0);

  /* First Uniform */
  for (i = 0; i < n; i++)
    mpfr_urandomb (tab[i], RANDS);
  algo_exact (real_non_rounded, tab, n, f);
  for (rnd_mode = 0; rnd_mode < MPFR_RND_MAX; rnd_mode++)
    {
      sum_tab (sum, tab, n, (mpfr_rnd_t) rnd_mode);
      mpfr_set (real_sum, real_non_rounded, (mpfr_rnd_t) rnd_mode);
      if (mpfr_cmp (real_sum, sum) != 0)
        {
          printf ("mpfr_sum incorrect.\n");
          mpfr_dump (real_sum);
          mpfr_dump (sum);
          exit (1);
        }
    }

  /* Then non uniform */
  for (i = 0; i < n; i++)
    {
      mpfr_urandomb (tab[i], RANDS);
      if (! mpfr_zero_p (tab[i]))
        mpfr_set_exp (tab[i], randlimb () % 1000);
    }
  algo_exact (real_non_rounded, tab, n, f);
  for (rnd_mode = 0; rnd_mode < MPFR_RND_MAX; rnd_mode++)
    {
      sum_tab (sum, tab, n, (mpfr_rnd_t) rnd_mode);
      mpfr_set (real_sum, real_non_rounded, (mpfr_rnd_t) rnd_mode);
      if (mpfr_cmp (real_sum, sum) != 0)
        {
          printf ("mpfr_sum incorrect.\n");
          mpfr_dump (real_sum);
          mpfr_dump (sum);
          exit (1);
        }
    }

  /* Clear stuff */
  for (i = 0; i < n; i++)
    mpfr_clear (tab[i]);
  mpfr_clears (sum, real_sum, real_non_rounded, (mpfr_ptr) 0);
  (*__gmp_free_func) (tab, n * sizeof(mpfr_t));
}

/* check case when c does not overlap with a, but both b and c count
   for rounding */
static void
check_case_1b (void)
{
    mpfr_t a, b, c;
    unsigned int prec_a, prec_b, prec_c, dif;

    mpfr_init (a);
    mpfr_init (b);
    mpfr_init (c);

    {
        prec_a = MPFR_PREC_MIN + (randlimb () % 63);
        mpfr_set_prec (a, prec_a);
        for (prec_b = prec_a + 2; prec_b <= 64; prec_b++)
        {
            dif = prec_b - prec_a;
            mpfr_set_prec (b, prec_b);
            /* b = 1 - 2^(-prec_a) + 2^(-prec_b) */
            mpfr_set_ui (b, 1, MPFR_RNDN);
            mpfr_div_2exp (b, b, dif, MPFR_RNDN);
            mpfr_sub_ui (b, b, 1, MPFR_RNDN);
            mpfr_div_2exp (b, b, prec_a, MPFR_RNDN);
            mpfr_add_ui (b, b, 1, MPFR_RNDN);
            for (prec_c = dif; prec_c <= 64; prec_c++)
            {
                /* c = 2^(-prec_a) - 2^(-prec_b) */
                mpfr_set_prec (c, prec_c);
                mpfr_set_si (c, -1, MPFR_RNDN);
                mpfr_div_2exp (c, c, dif, MPFR_RNDN);
                mpfr_add_ui (c, c, 1, MPFR_RNDN);
                mpfr_div_2exp (c, c, prec_a, MPFR_RNDN);
                test_add (a, b, c, MPFR_RNDN);
                if (mpfr_cmp_ui (a, 1) != 0)
                {
                    printf ("case (1b) failed for prec_a=%u, prec_b=%u,"
                            " prec_c=%u\n", prec_a, prec_b, prec_c);
                    printf ("b=");
                    mpfr_print_binary(b);
                    puts ("");
                    printf ("c=");
                    mpfr_print_binary(c);
                    puts ("");
                    printf ("a=");
                    mpfr_print_binary(a);
                    puts ("");
                    exit (1);
                }
            }
        }
    }

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

static void
check_inexact (void)
{
  mpfr_t x, y, z;
  mp_prec_t px, py;
  int inexact, cmp;
  unsigned long int u;
  int rnd;

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

  for (px=2; px<300; px++)
    {
      mpfr_set_prec (x, px);
      do
        {
          mpfr_urandomb (x, RANDS);
        }
      while (mpfr_cmp_ui (x, 0) == 0);
      u = randlimb ();
      for (py=2; py<300; py++)
        {
          mpfr_set_prec (y, py);
          mpfr_set_prec (z, py + px);
          for (rnd = 0; rnd < GMP_RND_MAX; rnd++)
            {
              inexact = mpfr_ui_div (y, u, x, (mp_rnd_t) rnd);
              if (mpfr_mul (z, y, x, (mp_rnd_t) rnd))
                {
                  printf ("z <- y * x should be exact\n");
                  exit (1);
                }
              cmp = mpfr_cmp_ui (z, u);
              if (((inexact == 0) && (cmp != 0)) ||
                  ((inexact > 0) && (cmp <= 0)) ||
                  ((inexact < 0) && (cmp >= 0)))
                {
                  printf ("Wrong inexact flag for u=%lu, rnd=%s\n",
                          u, mpfr_print_rnd_mode ((mp_rnd_t) rnd));
                  printf ("expected %d, got %d\n", cmp, inexact);
                  printf ("x="); mpfr_print_binary (x); puts ("");
                  printf ("y="); mpfr_print_binary (y); puts ("");
                  printf ("y*x="); mpfr_print_binary (z); puts ("");
                  exit (1);
                }
            }
        }
    }

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

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

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

  for (px=2; px<300; px++)
    {
      mpfr_set_prec (x, px);
      mpfr_urandomb (x, RANDS);
      do
        {
          u = randlimb ();
        }
      while (u == 0);
      for (py=2; py<300; py++)
        {
          mpfr_set_prec (y, py);
          mpfr_set_prec (z, py + mp_bits_per_limb);
          for (rnd = 0; rnd < MPFR_RND_MAX; rnd++)
            {
              inexact = mpfr_div_ui (y, x, u, (mpfr_rnd_t) rnd);
              if (mpfr_mul_ui (z, y, u, (mpfr_rnd_t) rnd))
                {
                  printf ("z <- y * u should be exact for u=%lu\n", u);
                  printf ("y="); mpfr_print_binary (y); puts ("");
                  printf ("z="); mpfr_print_binary (z); puts ("");
                  exit (1);
                }
              cmp = mpfr_cmp (z, x);
              if (((inexact == 0) && (cmp != 0)) ||
                  ((inexact > 0) && (cmp <= 0)) ||
                  ((inexact < 0) && (cmp >= 0)))
                {
                  printf ("Wrong inexact flag for u=%lu, rnd=%s\n", u,
                          mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
                  printf ("x="); mpfr_print_binary (x); puts ("");
                  printf ("y="); mpfr_print_binary (y); puts ("");
                  exit (1);
                }
            }
        }
    }

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

static void
check_set_uj (mpfr_prec_t pmin, mpfr_prec_t pmax, int N)
{
    mpfr_t x, y;
    mpfr_prec_t p;
    int inex1, inex2, n;
    mp_limb_t limb;

    mpfr_inits2 (pmax, x, y, (mpfr_ptr) 0);

    for ( p = pmin ; p < pmax ; p++)
    {
        mpfr_set_prec (x, p);
        mpfr_set_prec (y, p);
        for (n = 0 ; n < N ; n++)
        {
            /* mp_limb_t may be unsigned long long */
            limb = (unsigned long) randlimb ();
            inex1 = mpfr_set_uj (x, limb, MPFR_RNDN);
            inex2 = mpfr_set_ui (y, limb, MPFR_RNDN);
            if (mpfr_cmp (x, y))
            {
                printf ("ERROR for mpfr_set_uj and j=%lu and p=%lu\n",
                        (unsigned long) limb, (unsigned long) p);
                printf ("X=");
                mpfr_dump (x);
                printf ("Y=");
                mpfr_dump (y);
                exit (1);
            }
            if (inexact_sign (inex1) != inexact_sign (inex2))
            {
                printf ("ERROR for inexact(set_uj): j=%lu p=%lu\n"
                        "Inexact1= %d Inexact2= %d\n",
                        (unsigned long) limb, (unsigned long) p, inex1, inex2);
                exit (1);
            }
        }
    }
    /* Special case */
    mpfr_set_prec (x, sizeof(uintmax_t)*CHAR_BIT);
    inex1 = mpfr_set_uj (x, MPFR_UINTMAX_MAX, MPFR_RNDN);
    if (inex1 != 0 || mpfr_sgn(x) <= 0)
        ERROR ("inexact / UINTMAX_MAX");
    inex1 = mpfr_add_ui (x, x, 1, MPFR_RNDN);
    if (inex1 != 0 || !mpfr_powerof2_raw (x)
            || MPFR_EXP (x) != (sizeof(uintmax_t)*CHAR_BIT+1) )
        ERROR ("power of 2");
    mpfr_set_uj (x, 0, MPFR_RNDN);
    if (!MPFR_IS_ZERO (x))
        ERROR ("Setting 0");

    mpfr_clears (x, y, (mpfr_ptr) 0);
}

static void
set_special (mpfr_ptr x, unsigned int select)
{
  MPFR_ASSERTN (select < SPECIAL_MAX);
  switch (select)
    {
    case 0:
      MPFR_SET_NAN (x);
      break;
    case 1:
      MPFR_SET_INF (x);
      MPFR_SET_POS (x);
      break;
    case 2:
      MPFR_SET_INF (x);
      MPFR_SET_NEG (x);
      break;
    case 3:
      MPFR_SET_ZERO (x);
      MPFR_SET_POS  (x);
      break;
    case 4:
      MPFR_SET_ZERO (x);
      MPFR_SET_NEG  (x);
      break;
    case 5:
      mpfr_set_str_binary (x, "1");
      break;
    case 6:
      mpfr_set_str_binary (x, "-1");
      break;
    case 7:
      mpfr_set_str_binary (x, "1e-1");
      break;
    case 8:
      mpfr_set_str_binary (x, "1e+1");
      break;
    case 9:
      mpfr_const_pi (x, MPFR_RNDN);
      break;
    case 10:
      mpfr_const_pi (x, MPFR_RNDN);
      MPFR_SET_EXP (x, MPFR_GET_EXP (x)-1);
      break;
    default:
      mpfr_urandomb (x, RANDS);
      if (randlimb () & 1)
        mpfr_neg (x, x, MPFR_RNDN);
      break;
    }
}

static void
check_inexact (void)
{
  mpfr_prec_t p, q;
  mpfr_t x, y, absx;
  int rnd;
  int inexact, cmp;

  mpfr_init (x);
  mpfr_init (y);
  mpfr_init (absx);

  for (p=2; p<500; p++)
    {
      mpfr_set_prec (x, p);
      mpfr_set_prec (absx, p);
      mpfr_urandomb (x, RANDS);
      if (randlimb () % 2)
        {
          mpfr_set (absx, x, MPFR_RNDN);
          mpfr_neg (x, x, MPFR_RNDN);
        }
      else
        mpfr_set (absx, x, MPFR_RNDN);
      for (q=2; q<2*p; q++)
        {
          mpfr_set_prec (y, q);
          RND_LOOP (rnd)
            {
              inexact = mpfr_abs (y, x, (mpfr_rnd_t) rnd);
              cmp = mpfr_cmp (y, absx);
              if (((inexact == 0) && (cmp != 0)) ||
                  ((inexact > 0) && (cmp <= 0)) ||
                  ((inexact < 0) && (cmp >= 0)))
                {
                  printf ("Wrong inexact flag: expected %d, got %d\n",
                          cmp, inexact);
                  printf ("x="); mpfr_dump (x);
                  printf ("absx="); mpfr_dump (absx);
                  printf ("y="); mpfr_dump (y);
                  exit (1);
                }
            }
        }
    }

  mpfr_clear (x);
  mpfr_clear (y);
  mpfr_clear (absx);
}

static void
test_cmp_q (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
{
  mpfr_t x, z;
  mpq_t  y;
  mpfr_prec_t p;
  int res1, res2;
  int n;

  mpfr_init (x);
  mpfr_init2 (z, MPFR_PREC_MIN);
  mpq_init (y);
  for(p=pmin ; p < pmax ; p++)
    {
      mpfr_set_prec (x, p);
      for (n = 0 ; n < nmax ; n++)
        {
          mpfr_urandomb (x, RANDS);
          mpq_set_ui (y, randlimb (), randlimb() );
          if (!MPFR_IS_SINGULAR (x))
            {
              mpfr_sub_q (z, x, y, MPFR_RNDN);
              res1 = mpfr_sgn (z);
              res2 = mpfr_cmp_q (x, y);
              if (res1 != res2)
                {
                  printf("Error for mpfr_cmp_q: res=%d sub_z gives %d\n",
                         res2, res1);
                  exit (1);
                }
            }
        }
    }
  mpq_clear (y);
  mpfr_clear (x);
  mpfr_clear (z);
}

int
main (int argc, char *argv[])
{
  long j;

  tests_start_mpfr ();

  check (0, MPFR_RNDN);
  for (j = 0; j < 200000; j++)
    check (randlimb () & LONG_MAX, RND_RAND ());
  check0 ();

  tests_end_mpfr ();

  return 0;
}

static void
pow_int (mpfr_rnd_t rnd)
{
  mpfr_t ref1, ref2, ref3;
  mpfr_t res1;
  int i;

#ifdef DEBUG
  printf("pow_int\n");
#endif
  mpfr_inits2 ((randlimb () % 200) + MPFR_PREC_MIN,
               ref1, ref2, res1, (mpfr_ptr) 0);
  mpfr_init2 (ref3, 1005);

  for (i = 0; i <= 15; i++)
    {
      mpfr_urandomb (ref2, RANDS);
      if (i & 1)
        mpfr_neg (ref2, ref2, MPFR_RNDN);
      mpfr_set_ui (ref3, 20, MPFR_RNDN);
      /* We need to test huge integers because different algorithms/codes
         are used for not-too-large integers (mpfr_pow_z) and for general
         cases, in particular huge integers (mpfr_pow_general). [r7606] */
      if (i & 2)
        mpfr_mul_2ui (ref3, ref3, 1000, MPFR_RNDN);
      if (i & 4)
        mpfr_add_ui (ref3, ref3, 1, MPFR_RNDN);  /* odd integer */

      /* reference call: pow(a, b, c) */
      mpfr_pow (ref1, ref2, ref3, rnd);

      /* pow(a, a, c) */
      mpfr_set (res1, ref2, rnd); /* exact operation */
      mpfr_pow (res1, res1, ref3, rnd);

      if (mpfr_compare (res1, ref1))
        {
          printf ("Error for pow_int(a, a, c) for ");
          DISP("a=",ref2); DISP2(", c=",ref3);
          printf ("expected "); mpfr_print_binary (ref1); puts ("");
          printf ("got      "); mpfr_print_binary (res1); puts ("");
          exit (1);
        }
    }

  mpfr_clears (ref1, ref2, ref3, res1, (mpfr_ptr) 0);
}

static void
test2ui (int (*testfunc)(mpfr_ptr, mpfr_srcptr, unsigned long int, mp_rnd_t),
         char *foo, mp_prec_t prec, mp_rnd_t rnd)
{
  mpfr_t ref1, ref2;
  unsigned int ref3;
  mpfr_t res1;
  int i;

#ifdef DEBUG
  printf("checking %s\n", foo);
#endif
  mpfr_init2 (ref1, prec);
  mpfr_init2 (ref2, prec);
  mpfr_init2 (res1, prec);

  /* ref2 can be NaN, +Inf, -Inf, +0, -0 or any number
     ref3 can be 0 or any number */
  for (i=0; i<SPECIAL_MAX*2; i++)
    {
      set_special (ref2, i%SPECIAL_MAX);
      ref3 = i/SPECIAL_MAX == 0 ? 0 : randlimb ();

      /* reference call: foo(a, b, c) */
      testfunc (ref1, ref2, ref3, rnd);

      /* foo(a, a, c) */
      mpfr_set (res1, ref2, rnd); /* exact operation */
      testfunc (res1, res1, ref3, rnd);

      if (mpfr_compare (res1, ref1))
        {
          printf ("Error for %s(a, a, c) for c=%u\n", foo, ref3);
          DISP2("a=",ref2);
          printf ("expected "); mpfr_print_binary (ref1); puts ("");
          printf ("got      "); mpfr_print_binary (res1); puts ("");
          exit (1);
        }
    }

  mpfr_clear (ref1);
  mpfr_clear (ref2);
  mpfr_clear (res1);
}