C++ MPFR_IS_NAN函数代码示例

int
mpfr_max (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
{
  if (MPFR_IS_NAN(x) && MPFR_IS_NAN(y) ) 
    {  
      MPFR_SET_NAN(z);
      MPFR_RET_NAN;
    }
  MPFR_CLEAR_NAN(z);

  if (MPFR_IS_NAN(x))
    return mpfr_set(z, y, rnd_mode);

  if (MPFR_IS_NAN(y))
    return mpfr_set(z, x, rnd_mode);

  if (MPFR_IS_FP(x) && MPFR_IS_ZERO(x) && MPFR_IS_FP(y) && MPFR_IS_ZERO(y))
    {
      if (MPFR_SIGN(x) < 0)
        return mpfr_set(z, y, rnd_mode);
      else
        return mpfr_set(z, x, rnd_mode);
    }

  if (mpfr_cmp(x,y) <= 0)
    return mpfr_set(z, y, rnd_mode);
  else
    return mpfr_set(z, x, rnd_mode);
}

/* compute remainder as in definition:
   r = x - n * y, where n = trunc(x/y).
   warning: may change flags. */
static int
slow_fmod (mpfr_ptr r, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
{
  mpfr_t q;
  int inexact;
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x) || MPFR_IS_SINGULAR (y)))
    {
      if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y) || MPFR_IS_INF (x)
          || MPFR_IS_ZERO (y))
        {
          MPFR_SET_NAN (r);
          MPFR_RET_NAN;
        }
      else                      /* either y is Inf and x is 0 or non-special,
                                   or x is 0 and y is non-special,
                                   in both cases the quotient is zero. */
        return mpfr_set (r, x, rnd);
    }
  /* regular cases */
  /* if 2^(ex-1) <= |x| < 2^ex, and 2^(ey-1) <= |y| < 2^ey,
     then |x/y| < 2^(ex-ey+1) */
  mpfr_init2 (q,
              MAX (MPFR_PREC_MIN, mpfr_get_exp (x) - mpfr_get_exp (y) + 1));
  mpfr_div (q, x, y, MPFR_RNDZ);
  mpfr_trunc (q, q);                            /* may change inexact flag */
  mpfr_prec_round (q, mpfr_get_prec (q) + mpfr_get_prec (y), MPFR_RNDZ);
  inexact = mpfr_mul (q, q, y, MPFR_RNDZ);       /* exact */
  inexact = mpfr_sub (r, x, q, rnd);
  mpfr_clear (q);
  return inexact;
}

int
mpfr_max (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
{
  if (MPFR_ARE_SINGULAR(x,y))
    {
      if (MPFR_IS_NAN(x) && MPFR_IS_NAN(y) )
        {
          MPFR_SET_NAN(z);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_NAN(x))
        return mpfr_set(z, y, rnd_mode);
      else if (MPFR_IS_NAN(y))
        return mpfr_set(z, x, rnd_mode);
      else if (MPFR_IS_ZERO(x) && MPFR_IS_ZERO(y))
        {
          if (MPFR_IS_NEG(x))
            return mpfr_set(z, y, rnd_mode);
          else
            return mpfr_set(z, x, rnd_mode);
        }
    }
  if (mpfr_cmp(x,y) <= 0)
    return mpfr_set(z, y, rnd_mode);
  else
    return mpfr_set(z, x, rnd_mode);
}

int
mpfr_cmpabs (mpfr_srcptr b, mpfr_srcptr c)
{
  mp_exp_t be, ce;
  mp_size_t bn, cn;
  mp_limb_t *bp, *cp;

  if (MPFR_ARE_SINGULAR (b, c))
    {
      if (MPFR_IS_NAN (b) || MPFR_IS_NAN (c))
        {
          MPFR_SET_ERANGE ();
          return 0;
        }
      else if (MPFR_IS_INF (b))
        return ! MPFR_IS_INF (c);
      else if (MPFR_IS_INF (c))
        return -1;
      else if (MPFR_IS_ZERO (c))
        return ! MPFR_IS_ZERO (b);
      else /* b == 0 */
        return -1;
    }

  be = MPFR_GET_EXP (b);
  ce = MPFR_GET_EXP (c);
  if (be > ce)
    return 1;
  if (be < ce)
    return -1;

  /* exponents are equal */

  bn = MPFR_LIMB_SIZE(b)-1;
  cn = MPFR_LIMB_SIZE(c)-1;

  bp = MPFR_MANT(b);
  cp = MPFR_MANT(c);

  for ( ; bn >= 0 && cn >= 0; bn--, cn--)
    {
      if (bp[bn] > cp[cn])
        return 1;
      if (bp[bn] < cp[cn])
        return -1;
    }

  for ( ; bn >= 0; bn--)
    if (bp[bn])
      return 1;

  for ( ; cn >= 0; cn--)
    if (cp[cn])
      return -1;

   return 0;
}

/* When a * b is exact, the FMA is equivalent to the separate operations. */
static void
test_exact (void)
{
  char *val[] =
    { "@[email protected]", "[email protected]@", "-2", "-1", "-0", "0", "1", "2", "@[email protected]" };
  int sv = sizeof (val) / sizeof (*val);
  int i, j, k;
  int rnd;
  mpfr_t a, b, c, r1, r2;

  mpfr_inits2 (8, a, b, c, r1, r2, (void *) 0);

  for (i = 0; i < sv; i++)
    for (j = 0; j < sv; j++)
      for (k = 0; k < sv; k++)
        RND_LOOP (rnd)
          {
            if (mpfr_set_str (a, val[i], 10, GMP_RNDN) ||
                mpfr_set_str (b, val[j], 10, GMP_RNDN) ||
                mpfr_set_str (c, val[k], 10, GMP_RNDN) ||
                mpfr_mul (r1, a, b, rnd) ||
                mpfr_add (r1, r1, c, rnd))
              {
                printf ("test_exact internal error for (%d,%d,%d,%d)\n",
                        i, j, k, rnd);
                exit (1);
              }
            if (mpfr_fma (r2, a, b, c, rnd))
              {
                printf ("test_exact(%d,%d,%d,%d): mpfr_fma should be exact\n",
                        i, j, k, rnd);
                exit (1);
              }
            if (MPFR_IS_NAN (r1))
              {
                if (MPFR_IS_NAN (r2))
                  continue;
                printf ("test_exact(%d,%d,%d,%d): mpfr_fma should be NaN\n",
                        i, j, k, rnd);
                exit (1);
              }
            if (mpfr_cmp (r1, r2) || MPFR_SIGN (r1) != MPFR_SIGN (r2))
              {
                printf ("test_exact(%d,%d,%d,%d):\nexpected ", i, j, k, rnd);
                mpfr_out_str (stdout, 10, 0, r1, GMP_RNDN);
                printf ("\n     got ");
                mpfr_out_str (stdout, 10, 0, r2, GMP_RNDN);
                printf ("\n");
                exit (1);
              }
          }

  mpfr_clears (a, b, c, r1, r2, (void *) 0);
}

/* Arguments:
 *   spx: non-zero if px is a stringm zero if px is a MPFR number.
 *   px: value of x (string or MPFR number).
 *   sy: value of y (string).
 *   rnd: rounding mode.
 *   z1: expected result (null pointer if unknown pure FP value).
 *   inex1: expected ternary value (if z1 is not a null pointer).
 *   z2: computed result.
 *   inex2: computed ternary value.
 *   flags1: expected flags (computed flags in __gmpfr_flags).
 *   s1, s2: strings about the context.
 */
static void
cmpres (int spx, const void *px, const char *sy, mpfr_rnd_t rnd,
        mpfr_srcptr z1, int inex1, mpfr_srcptr z2, int inex2,
        unsigned int flags1, const char *s1, const char *s2)
{
  unsigned int flags2 = __gmpfr_flags;

  if (flags1 == flags2)
    {
      /* Note: the test on the sign of z1 and z2 is needed
         in case they are both zeros. */
      if (z1 == NULL)
        {
          if (MPFR_IS_PURE_FP (z2))
            return;
        }
      else if (SAME_SIGN (inex1, inex2) &&
               ((MPFR_IS_NAN (z1) && MPFR_IS_NAN (z2)) ||
                ((MPFR_IS_NEG (z1) ^ MPFR_IS_NEG (z2)) == 0 &&
                 mpfr_equal_p (z1, z2))))
        return;
    }

  printf ("Error in %s\nwith %s%s\nx = ", s1, s2,
          ext ? ", extended exponent range" : "");
  if (spx)
    printf ("%s, ", (char *) px);
  else
    {
      mpfr_out_str (stdout, 16, 0, (mpfr_ptr) px, MPFR_RNDN);
      puts (",");
    }
  printf ("y = %s, %s\n", sy, mpfr_print_rnd_mode (rnd));
  printf ("Expected ");
  if (z1 == NULL)
    {
      printf ("pure FP value, flags = %u\n", flags1);
    }
  else
    {
      mpfr_out_str (stdout, 16, 0, z1, MPFR_RNDN);
      printf (", inex = %d, flags = %u\n", SIGN (inex1), flags1);
    }
  printf ("Got      ");
  mpfr_out_str (stdout, 16, 0, z2, MPFR_RNDN);
  printf (", inex = %d, flags = %u\n", SIGN (inex2), flags2);
  if (all_cmpres_errors != 0)
    all_cmpres_errors = -1;
  else
    exit (1);
}

/* reldiff(b, c) = abs(b-c)/b */
void
mpfr_reldiff (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
  mpfr_t b_copy;

  if (MPFR_ARE_SINGULAR (b, c))
    {
      if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c))
        {
          MPFR_SET_NAN(a);
          return;
        }
      else if (MPFR_IS_INF(b))
        {
          if (MPFR_IS_INF (c) && (MPFR_SIGN (c) == MPFR_SIGN (b)))
            MPFR_SET_ZERO(a);
          else
            MPFR_SET_NAN(a);
          return;
        }
      else if (MPFR_IS_INF(c))
        {
          MPFR_SET_SAME_SIGN (a, b);
          MPFR_SET_INF (a);
          return;
        }
      else if (MPFR_IS_ZERO(b)) /* reldiff = abs(c)/c = sign(c) */
        {
          mpfr_set_si (a, MPFR_INT_SIGN (c), rnd_mode);
          return;
        }
      /* Fall through */
    }

  if (a == b)
    {
      mpfr_init2 (b_copy, MPFR_PREC(b));
      mpfr_set (b_copy, b, MPFR_RNDN);
    }

  mpfr_sub (a, b, c, rnd_mode);
  mpfr_abs (a, a, rnd_mode); /* for compatibility with MPF */
  mpfr_div (a, a, (a == b) ? b_copy : b, rnd_mode);

  if (a == b)
    mpfr_clear (b_copy);

}

/* sets x to x+sign(x)*ulp(x) */
int
mpfr_add_one_ulp (mpfr_ptr x, mp_rnd_t rnd_mode)
{
  mp_size_t xn;
  int sh;
  mp_limb_t *xp;

  if (MPFR_IS_NAN(x))
    MPFR_RET_NAN;

  if (MPFR_IS_INF(x) || MPFR_IS_ZERO(x))
    return 0;

  xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB;
  sh = xn * BITS_PER_MP_LIMB - MPFR_PREC(x);
  xp = MPFR_MANT(x);
  if (mpn_add_1 (xp, xp, xn, MP_LIMB_T_ONE << sh)) /* got 1.0000... */
    {
      mp_exp_t exp = MPFR_EXP(x);
      if (exp == __mpfr_emax)
        return mpfr_set_overflow(x, rnd_mode, MPFR_SIGN(x));
      else
        {
          MPFR_EXP(x)++;
          xp[xn-1] = GMP_LIMB_HIGHBIT;
        }
    }
  return 0;
}

int
mpfr_ui_div (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t uu;
  mp_limb_t up[1];
  unsigned long cnt;

  MPFR_LOG_FUNC
    (("u=%lu x[%Pu]=%.*Rg rnd=%d",
      u, mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode),
     ("y[%Pu]=%.*Rg", mpfr_get_prec(y), mpfr_log_prec, y));

  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)))
    {
      if (MPFR_IS_NAN(x))
        {
          MPFR_SET_NAN(y);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF(x)) /* u/Inf = 0 */
        {
          MPFR_SET_ZERO(y);
          MPFR_SET_SAME_SIGN(y,x);
          MPFR_RET(0);
        }
      else /* u / 0 */
        {
          MPFR_ASSERTD(MPFR_IS_ZERO(x));
          if (u)
            {
              /* u > 0, so y = sign(x) * Inf */
              MPFR_SET_SAME_SIGN(y, x);
              MPFR_SET_INF(y);
              mpfr_set_divby0 ();
              MPFR_RET(0);
            }
          else
            {
              /* 0 / 0 */
              MPFR_SET_NAN(y);
              MPFR_RET_NAN;
            }
        }
    }
  else if (MPFR_LIKELY(u != 0))
    {
      MPFR_TMP_INIT1(up, uu, GMP_NUMB_BITS);
      MPFR_ASSERTN(u == (mp_limb_t) u);
      count_leading_zeros(cnt, (mp_limb_t) u);
      up[0] = (mp_limb_t) u << cnt;
      MPFR_SET_EXP (uu, GMP_NUMB_BITS - cnt);
      return mpfr_div (y, uu, x, rnd_mode);
    }
  else /* u = 0, and x != 0 */
    {
      MPFR_SET_ZERO(y);         /* if u=0, then set y to 0 */
      MPFR_SET_SAME_SIGN(y, x); /* u considered as +0: sign(+0/x) = sign(x) */
      MPFR_RET(0);
    }
}

int
mpfr_frexp (mpfr_exp_t *exp, mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
{
  int inex;

  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)))
    {
      if (MPFR_IS_NAN(x))
        {
          MPFR_SET_NAN(y);
          MPFR_RET_NAN; /* exp is unspecified */
        }
      else if (MPFR_IS_INF(x))
        {
          MPFR_SET_INF(y);
          MPFR_SET_SAME_SIGN(y,x);
          MPFR_RET(0); /* exp is unspecified */
        }
      else
        {
          MPFR_SET_ZERO(y);
          MPFR_SET_SAME_SIGN(y,x);
          *exp = 0;
          MPFR_RET(0);
        }
    }

  inex = mpfr_set (y, x, rnd);
  *exp = MPFR_GET_EXP (y);
  MPFR_SET_EXP (y, 0);
  return mpfr_check_range (y, inex, rnd);
}

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
check_equal (mpfr_srcptr a, mpfr_srcptr a2, char *s,
             mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t r)
{
  if ((MPFR_IS_NAN (a) && MPFR_IS_NAN (a2)) ||
      mpfr_equal_p (a, a2))
    return;
  printf ("Error in %s\n", mpfr_print_rnd_mode (r));
  printf ("b  = ");
  mpfr_dump (b);
  printf ("c  = ");
  mpfr_dump (c);
  printf ("mpfr_div    result: ");
  mpfr_dump (a);
  printf ("%s result: ", s);
  mpfr_dump (a2);
  exit (1);
}

int
mpfr_dim (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
{
  if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y))
    {
      MPFR_SET_NAN(z);
      MPFR_RET_NAN;
    }

  if (mpfr_cmp (x,y) > 0)
    return mpfr_sub (z, x, y, rnd_mode);
  else
    {
      MPFR_SET_ZERO(z);
      MPFR_SET_POS(z);
      MPFR_RET(0);
    }
}

int
mpfr_ui_div (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mp_rnd_t rnd_mode)
{
  mpfr_t uu;
  mp_limb_t up[1];
  unsigned long cnt;

  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)))
    {
      if (MPFR_IS_NAN(x))
        {
          MPFR_SET_NAN(y);
          MPFR_RET_NAN;
        }
      else if (MPFR_IS_INF(x)) /* u/Inf = 0 */
        {
          MPFR_SET_ZERO(y);
          MPFR_SET_SAME_SIGN(y,x);
          MPFR_RET(0);
        }
      else /* u / 0 */
        {
          MPFR_ASSERTD(MPFR_IS_ZERO(x));
          if (u)
            {
              /* u > 0, so y = sign(x) * Inf */
              MPFR_SET_SAME_SIGN(y, x);
              MPFR_SET_INF(y);
              MPFR_RET(0);
            }
          else
            {
              /* 0 / 0 */
              MPFR_SET_NAN(y);
              MPFR_RET_NAN;
            }
        }
    }
  else if (MPFR_LIKELY(u != 0))
    {
      MPFR_TMP_INIT1(up, uu, BITS_PER_MP_LIMB);
      MPFR_ASSERTN(u == (mp_limb_t) u);
      count_leading_zeros(cnt, (mp_limb_t) u);
      up[0] = (mp_limb_t) u << cnt;
      MPFR_SET_EXP (uu, BITS_PER_MP_LIMB - cnt);
      return mpfr_div (y, uu, x, rnd_mode);
    }
  else /* u = 0, and x != 0 */
    {
      MPFR_SET_ZERO(y);         /* if u=0, then set y to 0 */
      MPFR_SET_SAME_SIGN(y, x); /* u considered as +0: sign(+0/x) = sign(x) */
      MPFR_RET(0);
    }
}

int
mpfr_custom_get_kind (mpfr_srcptr x)
{
  if (MPFR_LIKELY (!MPFR_IS_SINGULAR (x)))
    return (int) MPFR_REGULAR_KIND * MPFR_INT_SIGN (x);
  if (MPFR_IS_INF (x))
    return (int) MPFR_INF_KIND * MPFR_INT_SIGN (x);
  if (MPFR_IS_NAN (x))
    return (int) MPFR_NAN_KIND;
  MPFR_ASSERTD (MPFR_IS_ZERO (x));
  return (int) MPFR_ZERO_KIND * MPFR_INT_SIGN (x);
}