/* Test file for mpfr_sqrt.

Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007 Free Software Foundation, Inc.
Contributed by the Arenaire and Cacao projects, INRIA.

This file is part of the MPFR Library.

The 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 2.1 of the License, or (at your
option) any later version.

The 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 MPFR Library; see the file COPYING.LIB.  If not, 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_sqrt (mpfr_ptr a, mpfr_srcptr b, mp_rnd_t rnd_mode)
{
  int res;
  int ok = rnd_mode == GMP_RNDN && mpfr_number_p (b);
  if (ok)
    {
      mpfr_print_raw (b);
    }
  res = mpfr_sqrt (a, b, rnd_mode);
  if (ok)
    {
      printf (" ");
      mpfr_print_raw (a);
      printf ("\n");
    }
  return res;
}
#else
#define test_sqrt mpfr_sqrt
#endif

static void
check3 (const char *as, mp_rnd_t rnd_mode, const char *qs)
{
  mpfr_t q;

  mpfr_init2 (q, 53);
  mpfr_set_str1 (q, as);
  test_sqrt (q, q, rnd_mode);
  if (mpfr_cmp_str1 (q, qs) )
    {
      printf ("mpfr_sqrt failed for a=%s, rnd_mode=%s\n",
              as, mpfr_print_rnd_mode (rnd_mode));
      printf ("expected sqrt is %s, got ",qs);
      mpfr_out_str (stdout, 10, 0, q, GMP_RNDN);
      putchar ('\n');
      exit (1);
    }
  mpfr_clear (q);
}

static void
check4 (const char *as, mp_rnd_t rnd_mode, const char *Qs)
{
  mpfr_t q;

  mpfr_init2 (q, 53);
  mpfr_set_str1 (q, as);
  test_sqrt (q, q, rnd_mode);
  if (mpfr_cmp_str (q, Qs, 16, GMP_RNDN))
    {
      printf ("mpfr_sqrt failed for a=%s, rnd_mode=%s\n",
              as, mpfr_print_rnd_mode(rnd_mode));
      printf ("expected ");
      mpfr_out_str (stdout, 16, 0, q, GMP_RNDN);
      printf ("\ngot      %s\n", Qs);
      mpfr_clear (q);
      exit (1);
    }
  mpfr_clear (q);
}

static void
check24 (const char *as, mp_rnd_t rnd_mode, const char *qs)
{
  mpfr_t q;

  mpfr_init2 (q, 24);
  mpfr_set_str1 (q, as);
  test_sqrt (q, q, rnd_mode);
  if (mpfr_cmp_str1 (q, qs))
    {
      printf ("mpfr_sqrt failed for a=%s, prec=24, rnd_mode=%s\n",
              as, mpfr_print_rnd_mode(rnd_mode));
      printf ("expected sqrt is %s, got ",qs);
      mpfr_out_str (stdout, 10, 0, q, GMP_RNDN);
      printf ("\n");
      exit (1);
    }
  mpfr_clear (q);
}

static void
check_diverse (const char *as, mp_prec_t p, const char *qs)
{
  mpfr_t q;

  mpfr_init2 (q, p);
  mpfr_set_str1 (q, as);
  test_sqrt (q, q, GMP_RNDN);
  if (mpfr_cmp_str1 (q, qs))
    {
      printf ("mpfr_sqrt failed for a=%s, prec=%lu, rnd_mode=%s\n",
              as, p, mpfr_print_rnd_mode (GMP_RNDN));
      printf ("expected sqrt is %s, got ", qs);
      mpfr_out_str (stdout, 10, 0, q, GMP_RNDN);
      printf ("\n");
      exit (1);
    }
  mpfr_clear (q);
}

/* the following examples come from the paper "Number-theoretic Test
   Generation for Directed Rounding" from Michael Parks, Table 3 */
static void
check_float (void)
{
  check24("70368760954880.0", GMP_RNDN, "8.388609e6");
  check24("281474943156224.0", GMP_RNDN, "1.6777215e7");
  check24("70368777732096.0", GMP_RNDN, "8.388610e6");
  check24("281474909601792.0", GMP_RNDN, "1.6777214e7");
  check24("100216216748032.0", GMP_RNDN, "1.0010805e7");
  check24("120137273311232.0", GMP_RNDN, "1.0960715e7");
  check24("229674600890368.0", GMP_RNDN, "1.5155019e7");
  check24("70368794509312.0", GMP_RNDN, "8.388611e6");
  check24("281474876047360.0", GMP_RNDN, "1.6777213e7");
  check24("91214552498176.0", GMP_RNDN, "9.550631e6");

  check24("70368760954880.0", GMP_RNDZ, "8.388608e6");
  check24("281474943156224.0", GMP_RNDZ, "1.6777214e7");
  check24("70368777732096.0", GMP_RNDZ, "8.388609e6");
  check24("281474909601792.0", GMP_RNDZ, "1.6777213e7");
  check24("100216216748032.0", GMP_RNDZ, "1.0010805e7");
  check24("120137273311232.0", GMP_RNDZ, "1.0960715e7");
  check24("229674600890368.0", GMP_RNDZ, "1.5155019e7");
  check24("70368794509312.0", GMP_RNDZ, "8.38861e6");
  check24("281474876047360.0", GMP_RNDZ, "1.6777212e7");
  check24("91214552498176.0", GMP_RNDZ, "9.550631e6");

  check24("70368760954880.0", GMP_RNDU, "8.388609e6");
  check24("281474943156224.0",GMP_RNDU, "1.6777215e7");
  check24("70368777732096.0", GMP_RNDU, "8.388610e6");
  check24("281474909601792.0", GMP_RNDU, "1.6777214e7");
  check24("100216216748032.0", GMP_RNDU, "1.0010806e7");
  check24("120137273311232.0", GMP_RNDU, "1.0960716e7");
  check24("229674600890368.0", GMP_RNDU, "1.515502e7");
  check24("70368794509312.0", GMP_RNDU, "8.388611e6");
  check24("281474876047360.0", GMP_RNDU, "1.6777213e7");
  check24("91214552498176.0", GMP_RNDU, "9.550632e6");

  check24("70368760954880.0", GMP_RNDD, "8.388608e6");
  check24("281474943156224.0", GMP_RNDD, "1.6777214e7");
  check24("70368777732096.0", GMP_RNDD, "8.388609e6");
  check24("281474909601792.0", GMP_RNDD, "1.6777213e7");
  check24("100216216748032.0", GMP_RNDD, "1.0010805e7");
  check24("120137273311232.0", GMP_RNDD, "1.0960715e7");
  check24("229674600890368.0", GMP_RNDD, "1.5155019e7");
  check24("70368794509312.0", GMP_RNDD, "8.38861e6");
  check24("281474876047360.0", GMP_RNDD, "1.6777212e7");
  check24("91214552498176.0", GMP_RNDD, "9.550631e6");
}

static void
special (void)
{
  mpfr_t x, y, z;
  int inexact;
  mp_prec_t p;

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

  mpfr_set_prec (x, 64);
  mpfr_set_str_binary (x, "1010000010100011011001010101010010001100001101011101110001011001E-1");
  mpfr_set_prec (y, 32);
  test_sqrt (y, x, GMP_RNDN);
  if (mpfr_cmp_ui (y, 2405743844UL))
    {
      printf ("Error for n^2+n+1/2 with n=2405743843\n");
      exit (1);
    }

  mpfr_set_prec (x, 65);
  mpfr_set_str_binary (x, "10100000101000110110010101010100100011000011010111011100010110001E-2");
  mpfr_set_prec (y, 32);
  test_sqrt (y, x, GMP_RNDN);
  if (mpfr_cmp_ui (y, 2405743844UL))
    {
      printf ("Error for n^2+n+1/4 with n=2405743843\n");
      mpfr_dump (y);
      exit (1);
    }

  mpfr_set_prec (x, 66);
  mpfr_set_str_binary (x, "101000001010001101100101010101001000110000110101110111000101100011E-3");
  mpfr_set_prec (y, 32);
  test_sqrt (y, x, GMP_RNDN);
  if (mpfr_cmp_ui (y, 2405743844UL))
    {
      printf ("Error for n^2+n+1/4+1/8 with n=2405743843\n");
      mpfr_dump (y);
      exit (1);
    }

  mpfr_set_prec (x, 66);
  mpfr_set_str_binary (x, "101000001010001101100101010101001000110000110101110111000101100001E-3");
  mpfr_set_prec (y, 32);
  test_sqrt (y, x, GMP_RNDN);
  if (mpfr_cmp_ui (y, 2405743843UL))
    {
      printf ("Error for n^2+n+1/8 with n=2405743843\n");
      mpfr_dump (y);
      exit (1);
    }

  mpfr_set_prec (x, 27);
  mpfr_set_str_binary (x, "0.110100111010101000010001011");
  if ((inexact = test_sqrt (x, x, GMP_RNDZ)) >= 0)
    {
      printf ("Wrong inexact flag: expected -1, got %d\n", inexact);
      exit (1);
    }

  mpfr_set_prec (x, 2);
  for (p=2; p<1000; p++)
    {
      mpfr_set_prec (z, p);
      mpfr_set_ui (z, 1, GMP_RNDN);
      mpfr_nexttoinf (z);
      test_sqrt (x, z, GMP_RNDU);
      if (mpfr_cmp_ui_2exp(x, 3, -1))
        {
          printf ("Error: sqrt(1+ulp(1), up) should give 1.5 (prec=%u)\n",
                  (unsigned int) p);
          printf ("got "); mpfr_print_binary (x); puts ("");
          exit (1);
        }
    }

  /* check inexact flag */
  mpfr_set_prec (x, 5);
  mpfr_set_str_binary (x, "1.1001E-2");
  if ((inexact = test_sqrt (x, x, GMP_RNDN)))
    {
      printf ("Wrong inexact flag: expected 0, got %d\n", inexact);
      exit (1);
    }

  mpfr_set_prec (x, 2);
  mpfr_set_prec (z, 2);

  /* checks the sign is correctly set */
  mpfr_set_si (x, 1,  GMP_RNDN);
  mpfr_set_si (z, -1, GMP_RNDN);
  test_sqrt (z, x, GMP_RNDN);
  if (mpfr_cmp_ui (z, 0) < 0)
    {
      printf ("Error: square root of 1 gives ");
      mpfr_print_binary(z);
      putchar('\n');
      exit (1);
    }

  mpfr_set_prec (x, 192);
  mpfr_set_prec (z, 160);
  mpfr_set_str_binary (z, "0.1011010100000100100100100110011001011100100100000011000111011001011101101101110000110100001000100001100001011000E1");
  mpfr_set_prec (x, 160);
  test_sqrt(x, z, GMP_RNDN);
  test_sqrt(z, x, GMP_RNDN);

  mpfr_set_prec (x, 53);
  mpfr_set_str (x, "8093416094703476.0", 10, GMP_RNDN);
  mpfr_div_2exp (x, x, 1075, GMP_RNDN);
  test_sqrt (x, x, GMP_RNDN);
  mpfr_set_str (z, "1e55596835b5ef@-141", 16, GMP_RNDN);
  if (mpfr_cmp (x, z))
    {
      printf ("Error: square root of 8093416094703476*2^(-1075)\n");
      printf ("expected "); mpfr_dump (z);
      printf ("got      "); mpfr_dump (x);
      exit (1);
    }

  mpfr_set_prec (x, 33);
  mpfr_set_str_binary (x, "0.111011011011110001100111111001000e-10");
  mpfr_set_prec (z, 157);
  inexact = test_sqrt (z, x, GMP_RNDN);
  mpfr_set_prec (x, 157);
  mpfr_set_str_binary (x, "0.11110110101100101111001011100011100011100001101010111011010000100111011000111110100001001011110011111100101110010110010110011001011011010110010000011001101E-5");
  if (mpfr_cmp (x, z))
    {
      printf ("Error: square root (1)\n");
      exit (1);
    }
  if (inexact <= 0)
    {
      printf ("Error: wrong inexact flag (1)\n");
      exit (1);
    }

  /* case prec(result) << prec(input) */
  mpfr_set_prec (z, 2);
  for (p = 2; p < 1000; p++)
    {
      mpfr_set_prec (x, p);
      mpfr_set_ui (x, 1, GMP_RNDN);
      mpfr_nextabove (x);
      /* 1.0 < x <= 1.5 thus 1 < sqrt(x) <= 1.23 */
      inexact = test_sqrt (z, x, GMP_RNDN);
      MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0);
      inexact = test_sqrt (z, x, GMP_RNDZ);
      MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0);
      inexact = test_sqrt (z, x, GMP_RNDU);
      MPFR_ASSERTN(inexact > 0 && mpfr_cmp_ui_2exp (z, 3, -1) == 0);
      inexact = test_sqrt (z, x, GMP_RNDD);
      MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0);
    }

  /* corner case rw = 0 in rounding to nearest */
  mpfr_set_prec (z, BITS_PER_MP_LIMB - 1);
  mpfr_set_prec (y, BITS_PER_MP_LIMB - 1);
  mpfr_set_ui (y, 1, GMP_RNDN);
  mpfr_mul_2exp (y, y, BITS_PER_MP_LIMB - 1, GMP_RNDN);
  mpfr_nextabove (y);
  for (p = 2 * BITS_PER_MP_LIMB - 1; p <= 1000; p++)
    {
      mpfr_set_prec (x, p);
      mpfr_set_ui (x, 1, GMP_RNDN);
      mpfr_set_exp (x, BITS_PER_MP_LIMB);
      mpfr_add_ui (x, x, 1, GMP_RNDN);
      /* now x = 2^(BITS_PER_MP_LIMB - 1) + 1 (BITS_PER_MP_LIMB bits) */
      MPFR_ASSERTN (mpfr_mul (x, x, x, GMP_RNDN) == 0); /* exact */
      inexact = test_sqrt (z, x, GMP_RNDN);
      /* even rule: z should be 2^(BITS_PER_MP_LIMB - 1) */
      MPFR_ASSERTN (inexact < 0);
      MPFR_ASSERTN (mpfr_cmp_ui_2exp (z, 1, BITS_PER_MP_LIMB - 1) == 0);
      mpfr_nextbelow (x);
      /* now x is just below [2^(BITS_PER_MP_LIMB - 1) + 1]^2 */
      inexact = test_sqrt (z, x, GMP_RNDN);
      MPFR_ASSERTN(inexact < 0 &&
                   mpfr_cmp_ui_2exp (z, 1, BITS_PER_MP_LIMB - 1) == 0);
      mpfr_nextabove (x);
      mpfr_nextabove (x);
      /* now x is just above [2^(BITS_PER_MP_LIMB - 1) + 1]^2 */
      inexact = test_sqrt (z, x, GMP_RNDN);
      if (mpfr_cmp (z, y))
        {
          printf ("Error for sqrt(x) in rounding to nearest\n");
          printf ("x="); mpfr_dump (x);
          printf ("Expected "); mpfr_dump (y);
          printf ("Got      "); mpfr_dump (z);
          exit (1);
        }
      if (inexact <= 0)
        {
          printf ("Wrong inexact flag in corner case for p = %lu\n", (unsigned long) p);
          exit (1);
        }
    }

  mpfr_set_prec (x, 1000);
  mpfr_set_ui (x, 9, GMP_RNDN);
  mpfr_set_prec (y, 10);
  inexact = test_sqrt (y, x, GMP_RNDN);
  if (mpfr_cmp_ui (y, 3) || inexact != 0)
    {
      printf ("Error in sqrt(9:1000) for prec=10\n");
      exit (1);
    }
  mpfr_set_prec (y, BITS_PER_MP_LIMB);
  mpfr_nextabove (x);
  inexact = test_sqrt (y, x, GMP_RNDN);
  if (mpfr_cmp_ui (y, 3) || inexact >= 0)
    {
      printf ("Error in sqrt(9:1000) for prec=%u\n", BITS_PER_MP_LIMB);
      exit (1);
    }
  mpfr_set_prec (x, 2 * BITS_PER_MP_LIMB);
  mpfr_set_prec (y, BITS_PER_MP_LIMB);
  mpfr_set_ui (y, 1, GMP_RNDN);
  mpfr_nextabove (y);
  mpfr_set (x, y, GMP_RNDN);
  inexact = test_sqrt (y, x, GMP_RNDN);
  if (mpfr_cmp_ui (y, 1) || inexact >= 0)
    {
      printf ("Error in sqrt(1) for prec=%u\n", BITS_PER_MP_LIMB);
      mpfr_dump (y);
      exit (1);
    }

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

static void
check_inexact (mp_prec_t p)
{
  mpfr_t x, y, z;
  mp_rnd_t rnd;
  int inexact, sign;

  mpfr_init2 (x, p);
  mpfr_init2 (y, p);
  mpfr_init2 (z, 2*p);
  mpfr_random (x);
  rnd = (mp_rnd_t) RND_RAND();
  inexact = test_sqrt (y, x, rnd);
  if (mpfr_mul (z, y, y, rnd)) /* exact since prec(z) = 2*prec(y) */
    {
      printf ("Error: multiplication should be exact\n");
      exit (1);
    }
  mpfr_sub (z, z, x, rnd); /* exact also */
  sign = mpfr_cmp_ui (z, 0);
  if (((inexact == 0) && (sign)) ||
      ((inexact > 0) && (sign <= 0)) ||
      ((inexact < 0) && (sign >= 0)))
    {
      printf ("Error: wrong inexact flag, expected %d, got %d\n",
              sign, inexact);
      printf ("x=");
      mpfr_print_binary (x);
      printf (" rnd=%s\n", mpfr_print_rnd_mode (rnd));
      printf ("y="); mpfr_print_binary (y); puts ("");
      exit (1);
    }
  mpfr_clear (x);
  mpfr_clear (y);
  mpfr_clear (z);
}

static void
check_nan (void)
{
  mpfr_t  x, got;

  mpfr_init2 (x, 100L);
  mpfr_init2 (got, 100L);

  /* sqrt(NaN) == NaN */
  MPFR_CLEAR_FLAGS (x);
  MPFR_SET_NAN (x);
  MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_nan_p (got));

  /* sqrt(-1) == NaN */
  mpfr_set_si (x, -1L, GMP_RNDZ);
  MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_nan_p (got));

  /* sqrt(+inf) == +inf */
  MPFR_CLEAR_FLAGS (x);
  MPFR_SET_INF (x);
  MPFR_SET_POS (x);
  MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_inf_p (got));

  /* sqrt(-inf) == NaN */
  MPFR_CLEAR_FLAGS (x);
  MPFR_SET_INF (x);
  MPFR_SET_NEG (x);
  MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_nan_p (got));

  /* sqrt(-0) == 0 */
  mpfr_set_si (x, 0L, GMP_RNDZ);
  MPFR_SET_NEG (x);
  MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */
  MPFR_ASSERTN (mpfr_number_p (got));
  MPFR_ASSERTN (mpfr_cmp_ui (got, 0L) == 0);

  mpfr_clear (x);
  mpfr_clear (got);
}

/* check that -1 <= x/sqrt(x^2+y^2) <= 1 for rounding to nearest or up */
static void
test_property1 (mp_prec_t p, mp_rnd_t r)
{
  mpfr_t x, y, z, t;

  mpfr_init2 (x, p);
  mpfr_init2 (y, p);
  mpfr_init2 (z, p);
  mpfr_init2 (t, p);

  mpfr_random (x);
  mpfr_random (y);
  mpfr_mul (z, x, x, r);
  mpfr_mul (t, x, x, r);
  mpfr_add (z, z, t, r);
  mpfr_sqrt (z, z, r);
  mpfr_div (z, x, z, r);
  if (mpfr_cmp_si (z, -1) < 0 || mpfr_cmp_ui (z, 1) > 0)
    {
      printf ("Error, -1 <= x/sqrt(x^2+y^2) <= 1 does not hold for r=%s\n",
              mpfr_print_rnd_mode (r));
      printf ("x="); mpfr_dump (x);
      printf ("y="); mpfr_dump (y);
      printf ("got "); mpfr_dump (z);
      exit (1);
    }

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

/* check sqrt(x^2) = x */
static void
test_property2 (mp_prec_t p, mp_rnd_t r)
{
  mpfr_t x, y;

  mpfr_init2 (x, p);
  mpfr_init2 (y, p);

  mpfr_random (x);
  mpfr_mul (y, x, x, r);
  mpfr_sqrt (y, y, r);
  if (mpfr_cmp (y, x))
    {
      printf ("Error, sqrt(x^2) = x does not hold for r=%s\n",
              mpfr_print_rnd_mode (r));
      printf ("x="); mpfr_dump (x);
      printf ("got "); mpfr_dump (y);
      exit (1);
    }

  mpfr_clear (x);
  mpfr_clear (y);
}

int
main (void)
{
  mp_prec_t p;
  int k;

  tests_start_mpfr ();

  for (p = MPFR_PREC_MIN; p <= 128; p++)
    {
      test_property1 (p, GMP_RNDN);
      test_property1 (p, GMP_RNDU);
      test_property2 (p, GMP_RNDN);
    }

  check_diverse ("635030154261163106768013773815762607450069292760790610550915652722277604820131530404842415587328", 160, "796887792767063979679855997149887366668464780637");
  special ();
  check_nan ();

  for (p=2; p<200; p++)
    for (k=0; k<200; k++)
      check_inexact (p);
  check_float();

  check3 ("-0.0", GMP_RNDN, "0.0");
  check4 ("6.37983013646045901440e+32", GMP_RNDN, "5.9bc5036d09e0c@13");
  check4 ("1.0", GMP_RNDN, "1");
  check4 ("1.0", GMP_RNDZ, "1");
  check4 ("3.725290298461914062500000e-9", GMP_RNDN, "4@-4");
  check4 ("3.725290298461914062500000e-9", GMP_RNDZ, "4@-4");
  check4 ("1190456976439861.0", GMP_RNDZ, "2.0e7957873529a@6");
  check4 ("1219027943874417664.0", GMP_RNDZ, "4.1cf2af0e6a534@7");
  /* the following examples are bugs in Cygnus compiler/system, found by
     Fabrice Rouillier while porting mpfr to Windows */
  check4 ("9.89438396044940256501e-134", GMP_RNDU, "8.7af7bf0ebbee@-56");
  check4 ("7.86528588050363751914e+31", GMP_RNDZ, "1.f81fc40f32062@13");
  check4 ("0.99999999999999988897", GMP_RNDN, "f.ffffffffffff8@-1");
  check4 ("1.00000000000000022204", GMP_RNDN, "1");
  /* the following examples come from the paper "Number-theoretic Test
   Generation for Directed Rounding" from Michael Parks, Table 4 */

  check4 ("78652858805036375191418371571712.0", GMP_RNDN,
          "1.f81fc40f32063@13");
  check4 ("38510074998589467860312736661504.0", GMP_RNDN,
          "1.60c012a92fc65@13");
  check4 ("35318779685413012908190921129984.0", GMP_RNDN,
          "1.51d17526c7161@13");
  check4 ("26729022595358440976973142425600.0", GMP_RNDN,
          "1.25e19302f7e51@13");
  check4 ("22696567866564242819241453027328.0", GMP_RNDN,
          "1.0ecea7dd2ec3d@13");
  check4 ("22696888073761729132924856434688.0", GMP_RNDN,
          "1.0ecf250e8e921@13");
  check4 ("36055652513981905145251657416704.0", GMP_RNDN,
          "1.5552f3eedcf33@13");
  check4 ("30189856268896404997497182748672.0", GMP_RNDN,
          "1.3853ee10c9c99@13");
  check4 ("36075288240584711210898775080960.0", GMP_RNDN,
          "1.556abe212b56f@13");
  check4 ("72154663483843080704304789585920.0", GMP_RNDN,
          "1.e2d9a51977e6e@13");

  check4 ("78652858805036375191418371571712.0", GMP_RNDZ,
          "1.f81fc40f32062@13");
  check4 ("38510074998589467860312736661504.0", GMP_RNDZ,
          "1.60c012a92fc64@13");
  check4 ("35318779685413012908190921129984.0", GMP_RNDZ, "1.51d17526c716@13");
  check4 ("26729022595358440976973142425600.0", GMP_RNDZ, "1.25e19302f7e5@13");
  check4 ("22696567866564242819241453027328.0", GMP_RNDZ,
          "1.0ecea7dd2ec3c@13");
  check4 ("22696888073761729132924856434688.0", GMP_RNDZ, "1.0ecf250e8e92@13");
  check4 ("36055652513981905145251657416704.0", GMP_RNDZ,
          "1.5552f3eedcf32@13");
  check4 ("30189856268896404997497182748672.0", GMP_RNDZ,
          "1.3853ee10c9c98@13");
  check4 ("36075288240584711210898775080960.0", GMP_RNDZ,
          "1.556abe212b56e@13");
  check4 ("72154663483843080704304789585920.0", GMP_RNDZ,
          "1.e2d9a51977e6d@13");

  check4 ("78652858805036375191418371571712.0", GMP_RNDU,
          "1.f81fc40f32063@13");
  check4 ("38510074998589467860312736661504.0", GMP_RNDU,
          "1.60c012a92fc65@13");
  check4 ("35318779685413012908190921129984.0", GMP_RNDU,
          "1.51d17526c7161@13");
  check4 ("26729022595358440976973142425600.0", GMP_RNDU,
          "1.25e19302f7e51@13");
  check4 ("22696567866564242819241453027328.0", GMP_RNDU,
          "1.0ecea7dd2ec3d@13");
  check4 ("22696888073761729132924856434688.0", GMP_RNDU,
          "1.0ecf250e8e921@13");
  check4 ("36055652513981905145251657416704.0", GMP_RNDU,
          "1.5552f3eedcf33@13");
  check4 ("30189856268896404997497182748672.0", GMP_RNDU,
          "1.3853ee10c9c99@13");
  check4 ("36075288240584711210898775080960.0", GMP_RNDU,
          "1.556abe212b56f@13");
  check4 ("72154663483843080704304789585920.0", GMP_RNDU,
          "1.e2d9a51977e6e@13");

  check4 ("78652858805036375191418371571712.0", GMP_RNDD,
          "1.f81fc40f32062@13");
  check4 ("38510074998589467860312736661504.0", GMP_RNDD,
          "1.60c012a92fc64@13");
  check4 ("35318779685413012908190921129984.0", GMP_RNDD, "1.51d17526c716@13");
  check4 ("26729022595358440976973142425600.0", GMP_RNDD, "1.25e19302f7e5@13");
  check4 ("22696567866564242819241453027328.0", GMP_RNDD,
          "1.0ecea7dd2ec3c@13");
  check4 ("22696888073761729132924856434688.0", GMP_RNDD, "1.0ecf250e8e92@13");
  check4 ("36055652513981905145251657416704.0", GMP_RNDD,
          "1.5552f3eedcf32@13");
  check4 ("30189856268896404997497182748672.0", GMP_RNDD,
          "1.3853ee10c9c98@13");
  check4 ("36075288240584711210898775080960.0", GMP_RNDD,
          "1.556abe212b56e@13");
  check4 ("72154663483843080704304789585920.0", GMP_RNDD,
          "1.e2d9a51977e6d@13");

  tests_end_mpfr ();
  return 0;
}