static void
check_overflow (void)
{
  mpfr_t x, y, z1, z2;
  mpfr_exp_t emin, emax;

  emin = mpfr_get_emin ();
  emax = mpfr_get_emax ();

  set_emin (-1021);
  set_emax (1024);

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

  mpfr_set_str1 (x, "8.00468257869324898448e+307");
  mpfr_set_str1 (y, "7.44784712422708645156e+307");
  mpfr_add1sp (z1, x, y, MPFR_RNDN);
  mpfr_add1   (z2, x, y, MPFR_RNDN);
  if (mpfr_cmp (z1, z2))
    {
      printf ("Overflow bug in add1sp.\n");
      exit (1);
    }
  mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);

  set_emin (emin);
  set_emax (emax);
}

static void
special_overflow (void)
{
  mpfr_t x, y;
  mpfr_exp_t emin, emax;

  emin = mpfr_get_emin ();
  emax = mpfr_get_emax ();

  set_emin (-125);
  set_emax (128);
  mpfr_init2 (x, 24);
  mpfr_init2 (y, 48);
  mpfr_set_str_binary (x, "0.101100100000000000110100E0");
  mpfr_asin (y, x, MPFR_RNDN);
  if (mpfr_cmp_str (y, "0.110001001101001111110000010110001000111011001000E0",
                    2, MPFR_RNDN))
    {
      printf("Special Overflow error.\n");
      mpfr_dump (y);
      exit (1);
    }
  mpfr_clear (y);
  mpfr_clear (x);
  set_emin (emin);
  set_emax (emax);
}

/* This bug occurs in mpfr_exp_2 on a Linux-64 machine, r5475. */
static void
bug20080731 (void)
{
    mpfr_exp_t emin;
    mpfr_t x, y1, y2;
    mpfr_prec_t prec = 64;

    emin = mpfr_get_emin ();
    set_emin (MPFR_EMIN_MIN);

    mpfr_init2 (x, 200);
    mpfr_set_str (x, "[email protected]",
                  16, MPFR_RNDN);

    mpfr_init2 (y1, prec);
    mpfr_exp (y1, x, MPFR_RNDU);

    /* Compute the result with a higher internal precision. */
    mpfr_init2 (y2, 300);
    mpfr_exp (y2, x, MPFR_RNDU);
    mpfr_prec_round (y2, prec, MPFR_RNDU);

    if (mpfr_cmp0 (y1, y2) != 0)
    {
        printf ("Error in bug20080731\nExpected ");
        mpfr_out_str (stdout, 16, 0, y2, MPFR_RNDN);
        printf ("\nGot      ");
        mpfr_out_str (stdout, 16, 0, y1, MPFR_RNDN);
        printf ("\n");
        exit (1);
    }

    mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
    set_emin (emin);
}

static void
special_overflow (void)
{
  mpfr_t x, y;
  int inex;
  mpfr_exp_t emin, emax;

  emin = mpfr_get_emin ();
  emax = mpfr_get_emax ();

  set_emin (-125);
  set_emax (128);

  mpfr_init2 (x, 24);
  mpfr_init2 (y, 24);

  mpfr_set_str_binary (x, "0.101100100000000000110100E15");
  inex = mpfr_exp2 (y, x, MPFR_RNDN);
  if (!mpfr_inf_p (y) || inex <= 0)
    {
      printf ("Overflow error.\n");
      mpfr_dump (y);
      printf ("inex = %d\n", inex);
      exit (1);
    }

  mpfr_clear (y);
  mpfr_clear (x);
  set_emin (emin);
  set_emax (emax);
}

static void
check_emin_aux (mpfr_exp_t e)
{
  mpfr_t x;
  char *s1, s2[256];
  int i;
  mpfr_exp_t emin;
  mpz_t ee;

  MPFR_ASSERTN (e >= LONG_MIN);
  emin = mpfr_get_emin ();
  set_emin (e);

  mpfr_init2 (x, 16);
  mpz_init (ee);

  mpfr_setmin (x, e);
  mpz_set_si (ee, e);
  mpz_sub_ui (ee, ee, 1);

  i = mpfr_asprintf (&s1, "%Ra", x);
  MPFR_ASSERTN (i > 0);

  gmp_snprintf (s2, 256, "0x1p%Zd", ee);

  if (strcmp (s1, s2) != 0)
    {
      printf ("Error in check_emin_aux for emin = %ld\n", (long) e);
      printf ("Expected %s\n", s2);
      printf ("Got      %s\n", s1);
      exit (1);
    }

  mpfr_free_str (s1);

  i = mpfr_asprintf (&s1, "%Rb", x);
  MPFR_ASSERTN (i > 0);

  gmp_snprintf (s2, 256, "1p%Zd", ee);

  if (strcmp (s1, s2) != 0)
    {
      printf ("Error in check_emin_aux for emin = %ld\n", (long) e);
      printf ("Expected %s\n", s2);
      printf ("Got      %s\n", s1);
      exit (1);
    }

  mpfr_free_str (s1);

  mpfr_clear (x);
  mpz_clear (ee);
  set_emin (emin);
}

static void
special_overflow (void)
{
  /* Check for overflow in 3 cases:
     1. cosh(x) is representable, but not exp(x)
     2. cosh(x) is not representable in the selected range of exp.
     3. cosh(x) exp overflow even with the largest range of exp */
  mpfr_t x, y;
  mp_exp_t emin, emax;

  emin = mpfr_get_emin ();
  emax = mpfr_get_emax ();

  set_emin (-125);
  set_emax (128);

  mpfr_init2 (x, 24);
  mpfr_init2 (y, 24);

  mpfr_set_str_binary (x, "0.101100100000000000110100E7");
  mpfr_cosh (y, x, GMP_RNDN);
  if (mpfr_cmp_str (y, "0.101010001111001010001110E128", 2, GMP_RNDN))
    {
      printf("Special overflow error 1.\n");
      mpfr_dump (y);
      exit (1);
    }

  mpfr_set_str_binary (x, "0.101100100000000000110100E8");
  mpfr_cosh (y, x, GMP_RNDN);
  if (!mpfr_inf_p(y))
    {
      printf("Special overflow error 2.\n");
      mpfr_dump (y);
      exit (1);
    }

  set_emin (emin);
  set_emax (emax);

  mpfr_set_str_binary (x, "0.101100100000000000110100E1000000");
  mpfr_cosh (y, x, GMP_RNDN);
  if (!mpfr_inf_p(y))
    {
      printf("Special overflow error 3.\n");
      mpfr_dump (y);
      exit (1);
    }

  mpfr_clear (y);
  mpfr_clear (x);
}

static void
check_overflow (void)
{
  mpfr_t sum1, sum2, x, y;
  mpfr_ptr t[2 * NOVFL];
  mpfr_exp_t emin, emax;
  int i, r;

  emin = mpfr_get_emin ();
  emax = mpfr_get_emax ();
  set_emin (MPFR_EMIN_MIN);
  set_emax (MPFR_EMAX_MAX);

  mpfr_inits2 (32, sum1, sum2, x, y, (mpfr_ptr) 0);
  mpfr_setmax (x, mpfr_get_emax ());
  mpfr_neg (y, x, MPFR_RNDN);

  for (i = 0; i < 2 * NOVFL; i++)
    t[i] = i < NOVFL ? x : y;

  /* Two kinds of test:
   *   i = 1: overflow.
   *   i = 2: intermediate overflow (exact sum is 0).
   */
  for (i = 1; i <= 2; i++)
    RND_LOOP(r)
      {
        int inex1, inex2;

        inex1 = mpfr_add (sum1, x, i == 1 ? x : y, (mpfr_rnd_t) r);
        inex2 = mpfr_sum (sum2, t, i * NOVFL, (mpfr_rnd_t) r);
        MPFR_ASSERTN (mpfr_check (sum1));
        MPFR_ASSERTN (mpfr_check (sum2));
        if (!(mpfr_equal_p (sum1, sum2) && SAME_SIGN (inex1, inex2)))
          {
            printf ("Error in check_overflow on %s, i = %d\n",
                    mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
            printf ("Expected ");
            mpfr_dump (sum1);
            printf ("with inex = %d\n", inex1);
            printf ("Got      ");
            mpfr_dump (sum2);
            printf ("with inex = %d\n", inex2);
            exit (1);
          }
      }

  mpfr_clears (sum1, sum2, x, y, (mpfr_ptr) 0);

  set_emin (emin);
  set_emax (emax);
}

static void
check_large (void)
{
  mpz_t z;
  mpfr_t x, y;
  mpfr_exp_t emax, emin;

  mpz_init (z);
  mpfr_init2 (x, 160);
  mpfr_init2 (y, 160);

  mpz_set_str (z, "77031627725494291259359895954016675357279104942148788042", 10);
  mpfr_set_z (x, z, MPFR_RNDN);
  mpfr_set_str_binary (y, "0.1100100100001111110110101010001000100001011010001100001000110100110001001100011001100010100010111000000011011100000111001101000100101001000000100100111000001001E186");
  if (mpfr_cmp (x, y))
    {
      printf ("Error in mpfr_set_z on large input\n");
      exit (1);
    }

  /* check overflow */
  emax = mpfr_get_emax ();
  set_emax (2);
  mpz_set_str (z, "7", 10);
  mpfr_set_z (x, z, MPFR_RNDU);
  MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);
  set_emax (3);
  mpfr_set_prec (x, 2);
  mpz_set_str (z, "7", 10);
  mpfr_set_z (x, z, MPFR_RNDU);
  MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);
  set_emax (emax);

  /* check underflow */
  emin = mpfr_get_emin ();
  set_emin (3);
  mpz_set_str (z, "1", 10);
  mpfr_set_z (x, z, MPFR_RNDZ);
  MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS(x));
  set_emin (2);
  mpfr_set_z (x, z, MPFR_RNDN);
  MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS(x));
  set_emin (emin);

  mpz_clear (z);

  mpfr_clear (x);
  mpfr_clear (y);
}

static void
underflow (void)
{
    mpfr_exp_t emin;

    underflow_up (0);

    emin = mpfr_get_emin ();
    set_emin (MPFR_EMIN_MIN);
    if (mpfr_get_emin () != emin)
    {
        underflow_up (1);
        set_emin (emin);
    }
}

static void
check_emin_emax (void)
{
    mpfr_exp_t old_emin, old_emax;

    old_emin = mpfr_get_emin ();
    old_emax = mpfr_get_emax ();

    /* Check the functions not the macros ! */
    if ((mpfr_set_emin)(MPFR_EMIN_MIN) != 0)
        ERROR("set_emin failed!");
    if ((mpfr_get_emin)() != MPFR_EMIN_MIN)
        ERROR("get_emin FAILED!");
    if ((mpfr_set_emin)(MPFR_EMIN_MIN-1) == 0)
        ERROR("set_emin failed! (2)");

    if ((mpfr_set_emax)(MPFR_EMAX_MAX) != 0)
        ERROR("set_emax failed!");
    if ((mpfr_get_emax)() != MPFR_EMAX_MAX)
        ERROR("get_emax FAILED!");
    if ((mpfr_set_emax)(MPFR_EMAX_MAX+1) == 0)
        ERROR("set_emax failed! (2)");

    if ((mpfr_get_emin_min) () != MPFR_EMIN_MIN)
        ERROR ("get_emin_min");
    if ((mpfr_get_emin_max) () != MPFR_EMIN_MAX)
        ERROR ("get_emin_max");
    if ((mpfr_get_emax_min) () != MPFR_EMAX_MIN)
        ERROR ("get_emax_min");
    if ((mpfr_get_emax_max) () != MPFR_EMAX_MAX)
        ERROR ("get_emax_max");

    set_emin (old_emin);
    set_emax (old_emax);
}

static void
mpfr_set_double_range (void)
{
    mpfr_set_default_prec (54);
    if (mpfr_get_default_prec () != 54)
        ERROR ("get_default_prec failed (1)");
    mpfr_set_default_prec (53);
    if ((mpfr_get_default_prec) () != 53)
        ERROR ("get_default_prec failed (2)");

    /* in double precision format, the unbiased exponent is between 0 and
       2047, where 0 is used for subnormal numbers, and 2047 for special
       numbers (infinities, NaN), and the bias is 1023, thus "normal" numbers
       have an exponent between -1022 and 1023, corresponding to numbers
       between 2^(-1022) and previous(2^(1024)).
       (The smallest subnormal number is 0.(0^51)1*2^(-1022)= 2^(-1074).)

       The smallest normal power of two is 1.0*2^(-1022).
       The largest normal power of two is 2^1023.
       (We have to add one for mpfr since mantissa are between 1/2 and 1.)
    */

    set_emin (-1021);
    set_emax (1024);
}

/* bug found by Kevin P. Rauch on 22 Oct 2007 */
static void
check2 (void)
{
  mpfr_t x, y, z;
  int tern;
  mpfr_exp_t emin;

  emin = mpfr_get_emin ();

  mpfr_init2 (x, 32);
  mpfr_init2 (y, 32);
  mpfr_init2 (z, 32);

  mpfr_set_ui (x, 0xC0000000U, MPFR_RNDN);
  mpfr_neg (x, x, MPFR_RNDN);
  mpfr_set_ui (y, 0xFFFFFFFEU, MPFR_RNDN);
  mpfr_set_exp (x, 0);
  mpfr_set_exp (y, 0);
  mpfr_set_emin (-29);

  tern = mpfr_mul (z, x, y, MPFR_RNDN);
  /* z = -0.BFFFFFFE, tern > 0 */

  tern = mpfr_subnormalize (z, tern, MPFR_RNDN);
  /* z should be -0.75 */
  MPFR_ASSERTN (tern < 0 && mpfr_cmp_si_2exp (z, -3, -2) == 0);

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

  MPFR_ASSERTN (mpfr_get_emin () == -29);

  set_emin (emin);
}

/* https://sympa.inria.fr/sympa/arc/mpfr/2011-05/msg00008.html
 * Incorrect flags (in debug mode on a 32-bit machine, assertion failure).
 */
static void
reduced_expo_range (void)
{
    mpfr_exp_t emin, emax;
    mpfr_t x, y, ex_y;
    int inex, ex_inex;
    unsigned int flags, ex_flags;

    emin = mpfr_get_emin ();
    emax = mpfr_get_emax ();

    mpfr_inits2 (12, x, y, ex_y, (mpfr_ptr) 0);
    mpfr_set_str (x, "0.1e-5", 2, MPFR_RNDN);

    set_emin (-5);
    set_emax (-5);
    mpfr_clear_flags ();
    inex = mpfr_atan (y, x, MPFR_RNDN);
    flags = __gmpfr_flags;
    set_emin (emin);
    set_emax (emax);

    mpfr_set_str (ex_y, "0.1e-5", 2, MPFR_RNDN);
    ex_inex = 1;
    ex_flags = MPFR_FLAGS_INEXACT;

    if (SIGN (inex) != ex_inex || flags != ex_flags ||
            ! mpfr_equal_p (y, ex_y))
    {
        printf ("Error in reduced_expo_range\non x = ");
        mpfr_dump (x);
        printf ("Expected y = ");
        mpfr_out_str (stdout, 2, 0, ex_y, MPFR_RNDN);
        printf ("\n         inex = %d, flags = %u\n", ex_inex, ex_flags);
        printf ("Got      y = ");
        mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
        printf ("\n         inex = %d, flags = %u\n", SIGN (inex), flags);
        exit (1);
    }

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

/* bug found by Kevin P. Rauch on 22 Oct 2007 */
static void
check3 (void)
{
  mpfr_t x, y, z;
  int tern;
  mpfr_exp_t emin;

  emin = mpfr_get_emin ();

  mpfr_init2 (x, 32);
  mpfr_init2 (y, 32);
  mpfr_init2 (z, 32);

  mpfr_set_ui (x, 0xBFFFFFFFU, MPFR_RNDN); /* 3221225471/2^32 */
  mpfr_set_ui (y, 0x80000001U, MPFR_RNDN); /* 2147483649/2^32 */
  mpfr_set_exp (x, 0);
  mpfr_set_exp (y, 0);
  mpfr_set_emin (-1);

  /* the exact product is 6917529028714823679/2^64, which is rounded to
     3/8 = 0.375, which is smaller, thus tern < 0 */
  tern = mpfr_mul (z, x, y, MPFR_RNDN);
  MPFR_ASSERTN (tern < 0 && mpfr_cmp_ui_2exp (z, 3, -3) == 0);

  tern = mpfr_subnormalize (z, tern, MPFR_RNDN);
  /* since emin = -1, and EXP(z)=-1, z should be rounded to precision
     EXP(z)-emin+1 = 1, i.e., z should be a multiple of the smallest possible
     positive representable value with emin=-1, which is 1/4. The two
     possible values are 1/4 and 2/4, which are at equal distance of z.
     But since tern < 0, we should choose the largest value, i.e., 2/4. */
  MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 1, -1) == 0);

  /* here is another test for the alternate case, where z was rounded up
     first, thus we have to round down */
  mpfr_set_str_binary (x, "0.11111111111010110101011011011011");
  mpfr_set_str_binary (y, "0.01100000000001111100000000001110");
  tern = mpfr_mul (z, x, y, MPFR_RNDN);
  MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 3, -3) == 0);
  tern = mpfr_subnormalize (z, tern, MPFR_RNDN);
  MPFR_ASSERTN (tern < 0 && mpfr_cmp_ui_2exp (z, 1, -2) == 0);

  /* finally the case where z was exact, which we simulate here */
  mpfr_set_ui_2exp (z, 3, -3, MPFR_RNDN);
  tern = mpfr_subnormalize (z, 0, MPFR_RNDN);
  MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 1, -1) == 0);

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

  MPFR_ASSERTN (mpfr_get_emin () == -1);

  set_emin (emin);
}

int
main (int argc, char *argv[])
{
  mpfr_t x;
  int ret;
  mpfr_exp_t emin, emax;

  tests_start_mpfr ();

  emin = mpfr_get_emin ();
  emax = mpfr_get_emax ();

  mpfr_init (x);

  mpfr_set_ui (x, 1, MPFR_RNDN);
  ret = mpfr_set_exp (x, 2);
  MPFR_ASSERTN(ret == 0 && mpfr_cmp_ui (x, 2) == 0);

  set_emin (-1);
  ret = mpfr_set_exp (x, -1);
  MPFR_ASSERTN(ret == 0 && mpfr_cmp_ui_2exp (x, 1, -2) == 0);

  set_emax (1);
  ret = mpfr_set_exp (x, 1);
  MPFR_ASSERTN(ret == 0 && mpfr_cmp_ui (x, 1) == 0);

  ret = mpfr_set_exp (x, -2);
  MPFR_ASSERTN(ret != 0 && mpfr_cmp_ui (x, 1) == 0);

  ret = mpfr_set_exp (x, 2);
  MPFR_ASSERTN(ret != 0 && mpfr_cmp_ui (x, 1) == 0);

  mpfr_clear (x);

  set_emin (emin);
  set_emax (emax);

  tests_end_mpfr ();
  return 0;
}

static void
exp_range (void)
{
  mpfr_t x;
  mpfr_exp_t emin;

  emin = mpfr_get_emin ();
  set_emin (3);
  mpfr_init2 (x, 8);
  mpfr_set_ui (x, 5, MPFR_RNDN);
  mpfr_exp2 (x, x, MPFR_RNDN);
  set_emin (emin);
  if (mpfr_nan_p (x) || mpfr_cmp_ui (x, 32) != 0)
    {
      printf ("Error in mpfr_exp2 for x = 5, with emin = 3\n");
      printf ("Expected 32, got ");
      mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
      printf ("\n");
      exit (1);
    }
  mpfr_clear (x);
}

static void
special_overflow (void)
{
  mpfr_t x, y;
  mpfr_exp_t emin, emax;

  emin = mpfr_get_emin ();
  emax = mpfr_get_emax ();

  mpfr_init2 (x, 24);
  mpfr_init2 (y, 73);

  /* Check special case: An overflow in const_pi could occurs! */
  set_emin (-125);
  set_emax (128);
  mpfr_set_str_binary (x, "0.111101010110110011101101E6");
  test_cos (y, x, MPFR_RNDZ);
  set_emin (emin);
  set_emax (emax);

  mpfr_clear (x);
  mpfr_clear (y);
}

static void
check_flags (void)
{
    mpfr_t x;
    mpfr_exp_t old_emin, old_emax;

    old_emin = mpfr_get_emin ();
    old_emax = mpfr_get_emax ();
    mpfr_init (x);

    /* Check the functions not the macros ! */
    (mpfr_clear_flags)();
    mpfr_set_double_range ();

    mpfr_set_ui (x, 1, MPFR_RNDN);
    (mpfr_clear_overflow)();
    mpfr_mul_2exp (x, x, 1024, MPFR_RNDN);
    if (!(mpfr_overflow_p)())
        ERROR("ERROR: No overflow detected!\n");

    (mpfr_clear_underflow)();
    mpfr_set_ui (x, 1, MPFR_RNDN);
    mpfr_div_2exp (x, x, 1025, MPFR_RNDN);
    if (!(mpfr_underflow_p)())
        ERROR("ERROR: No underflow detected!\n");

    (mpfr_clear_nanflag)();
    MPFR_SET_NAN(x);
    mpfr_add (x, x, x, MPFR_RNDN);
    if (!(mpfr_nanflag_p)())
        ERROR("ERROR: No NaN flag!\n");

    (mpfr_clear_inexflag)();
    mpfr_set_ui(x, 2, MPFR_RNDN);
    mpfr_cos(x, x, MPFR_RNDN);
    if (!(mpfr_inexflag_p)())
        ERROR("ERROR: No inexact flag!\n");

    (mpfr_clear_erangeflag) ();
    mpfr_set_ui (x, 1, MPFR_RNDN);
    mpfr_mul_2exp (x, x, 1024, MPFR_RNDN);
    mpfr_get_ui (x, MPFR_RNDN);
    if (!(mpfr_erangeflag_p)())
        ERROR ("ERROR: No erange flag!\n");

    mpfr_clear (x);
    set_emin (old_emin);
    set_emax (old_emax);
}

static void
alltst (void)
{
  mpfr_exp_t emin, emax;

  ext = 0;
  test_small ();
  test_large_small ();
  check_overflow ();

  emin = mpfr_get_emin ();
  emax = mpfr_get_emax ();
  set_emin (MPFR_EMIN_MIN);
  set_emax (MPFR_EMAX_MAX);
  if (mpfr_get_emin () != emin || mpfr_get_emax () != emax)
    {
      ext = 1;
      test_small ();
      test_large_small ();
      check_overflow ();
      set_emin (emin);
      set_emax (emax);
    }
}

static void
alltst (void)
{
  mpfr_exp_t emin, emax;

  ext = 0;
  tst ();
  underflow_up ();
  overflow_inv ();

  emin = mpfr_get_emin ();
  emax = mpfr_get_emax ();
  set_emin (MPFR_EMIN_MIN);
  set_emax (MPFR_EMAX_MAX);
  if (mpfr_get_emin () != emin || mpfr_get_emax () != emax)
    {
      ext = 1;
      tst ();
      underflow_up ();
      overflow_inv ();
      set_emin (emin);
      set_emax (emax);
    }
}