• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

C++ MPFR_SIGN函数代码示例

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

本文整理汇总了C++中MPFR_SIGN函数的典型用法代码示例。如果您正苦于以下问题:C++ MPFR_SIGN函数的具体用法?C++ MPFR_SIGN怎么用?C++ MPFR_SIGN使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。



在下文中一共展示了MPFR_SIGN函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。

示例1: mpfr_rem1

static int
mpfr_rem1 (mpfr_ptr rem, long *quo, mpfr_rnd_t rnd_q,
           mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
{
  mpfr_exp_t ex, ey;
  int compare, inex, q_is_odd, sign, signx = MPFR_SIGN (x);
  mpz_t mx, my, r;
  int tiny = 0;

  MPFR_ASSERTD (rnd_q == MPFR_RNDN || rnd_q == MPFR_RNDZ);

  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))
        {
          /* for remquo, quo is undefined */
          MPFR_SET_NAN (rem);
          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. */
        {
          if (quo)
            *quo = 0;
          return mpfr_set (rem, x, rnd);
        }
    }

  /* now neither x nor y is NaN, Inf or zero */

  mpz_init (mx);
  mpz_init (my);
  mpz_init (r);

  ex = mpfr_get_z_2exp (mx, x);  /* x = mx*2^ex */
  ey = mpfr_get_z_2exp (my, y);  /* y = my*2^ey */

  /* to get rid of sign problems, we compute it separately:
     quo(-x,-y) = quo(x,y), rem(-x,-y) = -rem(x,y)
     quo(-x,y) = -quo(x,y), rem(-x,y)  = -rem(x,y)
     thus quo = sign(x/y)*quo(|x|,|y|), rem = sign(x)*rem(|x|,|y|) */
  sign = (signx == MPFR_SIGN (y)) ? 1 : -1;
  mpz_abs (mx, mx);
  mpz_abs (my, my);
  q_is_odd = 0;

  /* divide my by 2^k if possible to make operations mod my easier */
  {
    unsigned long k = mpz_scan1 (my, 0);
    ey += k;
    mpz_fdiv_q_2exp (my, my, k);
  }

  if (ex <= ey)
    {
      /* q = x/y = mx/(my*2^(ey-ex)) */

      /* First detect cases where q=0, to avoid creating a huge number
         my*2^(ey-ex): if sx = mpz_sizeinbase (mx, 2) and sy =
         mpz_sizeinbase (my, 2), we have x < 2^(ex + sx) and
         y >= 2^(ey + sy - 1), thus if ex + sx <= ey + sy - 1
         the quotient is 0 */
      if (ex + (mpfr_exp_t) mpz_sizeinbase (mx, 2) <
          ey + (mpfr_exp_t) mpz_sizeinbase (my, 2))
        {
          tiny = 1;
          mpz_set (r, mx);
          mpz_set_ui (mx, 0);
        }
      else
        {
          mpz_mul_2exp (my, my, ey - ex);   /* divide mx by my*2^(ey-ex) */

          /* since mx > 0 and my > 0, we can use mpz_tdiv_qr in all cases */
          mpz_tdiv_qr (mx, r, mx, my);
          /* 0 <= |r| <= |my|, r has the same sign as mx */
        }

      if (rnd_q == MPFR_RNDN)
        q_is_odd = mpz_tstbit (mx, 0);
      if (quo)                  /* mx is the quotient */
        {
          mpz_tdiv_r_2exp (mx, mx, WANTED_BITS);
          *quo = mpz_get_si (mx);
        }
    }
  else                          /* ex > ey */
    {
      if (quo) /* remquo case */
        /* for remquo, to get the low WANTED_BITS more bits of the quotient,
           we first compute R =  X mod Y*2^WANTED_BITS, where X and Y are
           defined below. Then the low WANTED_BITS of the quotient are
           floor(R/Y). */
        mpz_mul_2exp (my, my, WANTED_BITS);     /* 2^WANTED_BITS*Y */

      else if (rnd_q == MPFR_RNDN) /* remainder case */
        /* Let X = mx*2^(ex-ey) and Y = my. Then both X and Y are integers.
           Assume X = R mod Y, then x = X*2^ey = R*2^ey mod (Y*2^ey=y).
//.........这里部分代码省略.........
开发者ID:BrianGladman,项目名称:mpfr,代码行数:101,代码来源:rem1.c


示例2: log


//.........这里部分代码省略.........
#ifdef MPFR_JN
          mpfr_mul (c, c, Q, MPFR_RNDN); /* Q * (sin - cos) */
          mpfr_mul (s, s, P, MPFR_RNDN); /* P * (sin + cos) */
#else
          mpfr_mul (c, c, P, MPFR_RNDN); /* P * (sin - cos) */
          mpfr_mul (s, s, Q, MPFR_RNDN); /* Q * (sin + cos) */
#endif
          err = MPFR_EXP(c);
          if (MPFR_EXP(s) > err)
            err = MPFR_EXP(s);
#ifdef MPFR_JN
          mpfr_sub (s, s, c, MPFR_RNDN);
#else
          mpfr_add (s, s, c, MPFR_RNDN);
#endif
        }
      else /* n odd: P * (sin - cos) + Q (cos + sin) for jn,
                     Q * (sin - cos) - P (cos + sin) for yn */
        {
#ifdef MPFR_JN
          mpfr_mul (c, c, P, MPFR_RNDN); /* P * (sin - cos) */
          mpfr_mul (s, s, Q, MPFR_RNDN); /* Q * (sin + cos) */
#else
          mpfr_mul (c, c, Q, MPFR_RNDN); /* Q * (sin - cos) */
          mpfr_mul (s, s, P, MPFR_RNDN); /* P * (sin + cos) */
#endif
          err = MPFR_EXP(c);
          if (MPFR_EXP(s) > err)
            err = MPFR_EXP(s);
#ifdef MPFR_JN
          mpfr_add (s, s, c, MPFR_RNDN);
#else
          mpfr_sub (s, c, s, MPFR_RNDN);
#endif
        }
      if ((n & 2) != 0)
        mpfr_neg (s, s, MPFR_RNDN);
      if (MPFR_EXP(s) > err)
        err = MPFR_EXP(s);
      /* the absolute error on s is bounded by P*err(s/c) + Q*err(s/c)
         + err(P)*(s/c) + err(Q)*(s/c) + 3 * 2^(err - w - 1)
         <= (|P|+|Q|) * 2^(1-w) + err_s * 2^(1-w) + 2^err * 2^(1-w),
         since |c|, |old_s| <= 2. */
      err2 = (MPFR_EXP(P) >= MPFR_EXP(Q)) ? MPFR_EXP(P) + 2 : MPFR_EXP(Q) + 2;
      /* (|P| + |Q|) * 2^(1 - w) <= 2^(err2 - w) */
      err = MPFR_EXP(err_s) >= err ? MPFR_EXP(err_s) + 2 : err + 2;
      /* err_s * 2^(1-w) + 2^old_err * 2^(1-w) <= 2^err * 2^(-w) */
      err2 = (err >= err2) ? err + 1 : err2 + 1;
      /* now the absolute error on s is bounded by 2^(err2 - w) */

      /* multiply by sqrt(1/(Pi*z)) */
      mpfr_const_pi (c, MPFR_RNDN);     /* Pi, err <= 1 */
      mpfr_mul (c, c, z, MPFR_RNDN);    /* err <= 2 */
      mpfr_si_div (c, MPFR_IS_POS(z) ? 1 : -1, c, MPFR_RNDN); /* err <= 3 */
      mpfr_sqrt (c, c, MPFR_RNDN);      /* err<=5/2, thus the absolute error is
                                          bounded by 3*u*|c| for |u| <= 0.25 */
      mpfr_mul (err_t, c, s, MPFR_SIGN(c)==MPFR_SIGN(s) ? MPFR_RNDU : MPFR_RNDD);
      mpfr_abs (err_t, err_t, MPFR_RNDU);
      mpfr_mul_ui (err_t, err_t, 3, MPFR_RNDU);
      /* 3*2^(-w)*|old_c|*|s| [see below] is bounded by err_t * 2^(-w) */
      err2 += MPFR_EXP(c);
      /* |old_c| * 2^(err2 - w) [see below] is bounded by 2^(err2-w) */
      mpfr_mul (c, c, s, MPFR_RNDN);    /* the absolute error on c is bounded by
                                          1/2 ulp(c) + 3*2^(-w)*|old_c|*|s|
                                          + |old_c| * 2^(err2 - w) */
      /* compute err_t * 2^(-w) + 1/2 ulp(c) = (err_t + 2^EXP(c)) * 2^(-w) */
      err = (MPFR_EXP(err_t) > MPFR_EXP(c)) ? MPFR_EXP(err_t) + 1 : MPFR_EXP(c) + 1;
      /* err_t * 2^(-w) + 1/2 ulp(c) <= 2^(err - w) */
      /* now err_t * 2^(-w) bounds 1/2 ulp(c) + 3*2^(-w)*|old_c|*|s| */
      err = (err >= err2) ? err + 1 : err2 + 1;
      /* the absolute error on c is bounded by 2^(err - w) */

      mpfr_clear (s);
      mpfr_clear (P);
      mpfr_clear (Q);
      mpfr_clear (t);
      mpfr_clear (iz);
      mpfr_clear (err_t);
      mpfr_clear (err_s);
      mpfr_clear (err_u);

      err -= MPFR_EXP(c);
      if (MPFR_LIKELY (MPFR_CAN_ROUND (c, w - err, MPFR_PREC(res), r)))
        break;
      if (diverge != 0)
        {
          mpfr_set (c, z, r); /* will force inex=0 below, which means the
                               asymptotic expansion failed */
          break;
        }
      MPFR_ZIV_NEXT (loop, w);
    }
  MPFR_ZIV_FREE (loop);

  inex = (MPFR_IS_POS(z) || ((n & 1) == 0)) ? mpfr_set (res, c, r)
    : mpfr_neg (res, c, r);
  mpfr_clear (c);

  return inex;
}
开发者ID:Akheon23,项目名称:chromecast-mirrored-source.toolchain,代码行数:101,代码来源:jyn_asympt.c


示例3: mpfr_digamma

int
mpfr_digamma (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  int inex;
  MPFR_SAVE_EXPO_DECL (expo);

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


  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))
        {
          if (MPFR_IS_POS(x)) /* Digamma(+Inf) = +Inf */
            {
              MPFR_SET_SAME_SIGN(y, x);
              MPFR_SET_INF(y);
              MPFR_RET(0);
            }
          else                /* Digamma(-Inf) = NaN */
            {
              MPFR_SET_NAN(y);
              MPFR_RET_NAN;
            }
        }
      else /* Zero case */
        {
          /* the following works also in case of overlap */
          MPFR_SET_INF(y);
          MPFR_SET_OPPOSITE_SIGN(y, x);
          mpfr_set_divby0 ();
          MPFR_RET(0);
        }
    }

  /* Digamma is undefined for negative integers */
  if (MPFR_IS_NEG(x) && mpfr_integer_p (x))
    {
      MPFR_SET_NAN(y);
      MPFR_RET_NAN;
    }

  /* now x is a normal number */

  MPFR_SAVE_EXPO_MARK (expo);
  /* for x very small, we have Digamma(x) = -1/x - gamma + O(x), more precisely
     -1 < Digamma(x) + 1/x < 0 for -0.2 < x < 0.2, thus:
     (i) either x is a power of two, then 1/x is exactly representable, and
         as long as 1/2*ulp(1/x) > 1, we can conclude;
     (ii) otherwise assume x has <= n bits, and y has <= n+1 bits, then
   |y + 1/x| >= 2^(-2n) ufp(y), where ufp means unit in first place.
   Since |Digamma(x) + 1/x| <= 1, if 2^(-2n) ufp(y) >= 2, then
   |y - Digamma(x)| >= 2^(-2n-1)ufp(y), and rounding -1/x gives the correct result.
   If x < 2^E, then y > 2^(-E), thus ufp(y) > 2^(-E-1).
   A sufficient condition is thus EXP(x) <= -2 MAX(PREC(x),PREC(Y)). */
  if (MPFR_EXP(x) < -2)
    {
      if (MPFR_EXP(x) <= -2 * (mpfr_exp_t) MAX(MPFR_PREC(x), MPFR_PREC(y)))
        {
          int signx = MPFR_SIGN(x);
          inex = mpfr_si_div (y, -1, x, rnd_mode);
          if (inex == 0) /* x is a power of two */
            { /* result always -1/x, except when rounding down */
              if (rnd_mode == MPFR_RNDA)
                rnd_mode = (signx > 0) ? MPFR_RNDD : MPFR_RNDU;
              if (rnd_mode == MPFR_RNDZ)
                rnd_mode = (signx > 0) ? MPFR_RNDU : MPFR_RNDD;
              if (rnd_mode == MPFR_RNDU)
                inex = 1;
              else if (rnd_mode == MPFR_RNDD)
                {
                  mpfr_nextbelow (y);
                  inex = -1;
                }
              else /* nearest */
                inex = 1;
            }
          MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
          goto end;
        }
    }

  if (MPFR_IS_NEG(x))
    inex = mpfr_digamma_reflection (y, x, rnd_mode);
  /* if x < 1/2 we use the reflection formula */
  else if (MPFR_EXP(x) < 0)
    inex = mpfr_digamma_reflection (y, x, rnd_mode);
  else
    inex = mpfr_digamma_positive (y, x, rnd_mode);

 end:
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (y, inex, rnd_mode);
//.........这里部分代码省略.........
开发者ID:epowers,项目名称:mpfr,代码行数:101,代码来源:digamma.c


示例4: mpfr_acos

int
mpfr_acos (mpfr_ptr acos, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t xp, arcc, tmp;
  mpfr_exp_t supplement;
  mpfr_prec_t prec;
  int sign, compared, inexact;
  MPFR_SAVE_EXPO_DECL (expo);
  MPFR_ZIV_DECL (loop);

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

  /* Singular cases */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
        {
          MPFR_SET_NAN (acos);
          MPFR_RET_NAN;
        }
      else /* necessarily x=0 */
        {
          MPFR_ASSERTD(MPFR_IS_ZERO(x));
          /* acos(0)=Pi/2 */
          MPFR_SAVE_EXPO_MARK (expo);
          inexact = mpfr_const_pi (acos, rnd_mode);
          mpfr_div_2ui (acos, acos, 1, rnd_mode); /* exact */
          MPFR_SAVE_EXPO_FREE (expo);
          return mpfr_check_range (acos, inexact, rnd_mode);
        }
    }

  /* Set x_p=|x| */
  sign = MPFR_SIGN (x);
  mpfr_init2 (xp, MPFR_PREC (x));
  mpfr_abs (xp, x, MPFR_RNDN); /* Exact */

  compared = mpfr_cmp_ui (xp, 1);

  if (MPFR_UNLIKELY (compared >= 0))
    {
      mpfr_clear (xp);
      if (compared > 0) /* acos(x) = NaN for x > 1 */
        {
          MPFR_SET_NAN(acos);
          MPFR_RET_NAN;
        }
      else
        {
          if (MPFR_IS_POS_SIGN (sign)) /* acos(+1) = +0 */
            return mpfr_set_ui (acos, 0, rnd_mode);
          else /* acos(-1) = Pi */
            return mpfr_const_pi (acos, rnd_mode);
        }
    }

  MPFR_SAVE_EXPO_MARK (expo);

  /* Compute the supplement */
  mpfr_ui_sub (xp, 1, xp, MPFR_RNDD);
  if (MPFR_IS_POS_SIGN (sign))
    supplement = 2 - 2 * MPFR_GET_EXP (xp);
  else
    supplement = 2 - MPFR_GET_EXP (xp);
  mpfr_clear (xp);

  prec = MPFR_PREC (acos);
  prec += MPFR_INT_CEIL_LOG2(prec) + 10 + supplement;

  /* VL: The following change concerning prec comes from r3145
     "Optimize mpfr_acos by choosing a better initial precision."
     but it doesn't seem to be correct and leads to problems (assertion
     failure or very important inefficiency) with tiny arguments.
     Therefore, I've disabled it. */
  /* If x ~ 2^-N, acos(x) ~ PI/2 - x - x^3/6
     If Prec < 2*N, we can't round since x^3/6 won't be counted. */
#if 0
  if (MPFR_PREC (acos) >= MPFR_PREC (x) && MPFR_GET_EXP (x) < 0)
    {
      mpfr_uexp_t pmin = (mpfr_uexp_t) (-2 * MPFR_GET_EXP (x)) + 5;
      MPFR_ASSERTN (pmin <= MPFR_PREC_MAX);
      if (prec < pmin)
        prec = pmin;
    }
#endif

  mpfr_init2 (tmp, prec);
  mpfr_init2 (arcc, prec);

  MPFR_ZIV_INIT (loop, prec);
  for (;;)
    {
      /* acos(x) = Pi/2 - asin(x) = Pi/2 - atan(x/sqrt(1-x^2)) */
      mpfr_sqr (tmp, x, MPFR_RNDN);
      mpfr_ui_sub (tmp, 1, tmp, MPFR_RNDN);
      mpfr_sqrt (tmp, tmp, MPFR_RNDN);
      mpfr_div (tmp, x, tmp, MPFR_RNDN);
//.........这里部分代码省略.........
开发者ID:MiKTeX,项目名称:miktex,代码行数:101,代码来源:acos.c


示例5: overfl_exp10_0

static void
overfl_exp10_0 (void)
{
  mpfr_t x, y;
  int emax, i, inex, rnd, err = 0;
  mpfr_exp_t old_emax;

  old_emax = mpfr_get_emax ();

  mpfr_init2 (x, 8);
  mpfr_init2 (y, 8);

  for (emax = -1; emax <= 0; emax++)
    {
      mpfr_set_ui_2exp (y, 1, emax, MPFR_RNDN);
      mpfr_nextbelow (y);
      set_emax (emax);  /* 1 is not representable. */
      /* and if emax < 0, 1 - eps is not representable either. */
      for (i = -1; i <= 1; i++)
        RND_LOOP (rnd)
          {
            mpfr_set_si_2exp (x, i, -512 * ABS (i), MPFR_RNDN);
            mpfr_clear_flags ();
            inex = mpfr_exp10 (x, x, (mpfr_rnd_t) rnd);
            if ((i >= 0 || emax < 0 || rnd == MPFR_RNDN || rnd == MPFR_RNDU) &&
                ! mpfr_overflow_p ())
              {
                printf ("Error in overfl_exp10_0 (i = %d, rnd = %s):\n"
                        "  The overflow flag is not set.\n",
                        i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
                err = 1;
              }
            if (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)
              {
                if (inex >= 0)
                  {
                    printf ("Error in overfl_exp10_0 (i = %d, rnd = %s):\n"
                            "  The inexact value must be negative.\n",
                            i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
                    err = 1;
                  }
                if (! mpfr_equal_p (x, y))
                  {
                    printf ("Error in overfl_exp10_0 (i = %d, rnd = %s):\n"
                            "  Got ", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
                    mpfr_print_binary (x);
                    printf (" instead of 0.11111111E%d.\n", emax);
                    err = 1;
                  }
              }
            else
              {
                if (inex <= 0)
                  {
                    printf ("Error in overfl_exp10_0 (i = %d, rnd = %s):\n"
                            "  The inexact value must be positive.\n",
                            i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
                    err = 1;
                  }
                if (! (mpfr_inf_p (x) && MPFR_SIGN (x) > 0))
                  {
                    printf ("Error in overfl_exp10_0 (i = %d, rnd = %s):\n"
                            "  Got ", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
                    mpfr_print_binary (x);
                    printf (" instead of +Inf.\n");
                    err = 1;
                  }
              }
          }
      set_emax (old_emax);
    }

  if (err)
    exit (1);
  mpfr_clear (x);
  mpfr_clear (y);
}
开发者ID:michalkonecny,项目名称:haskell-mpfr,代码行数:77,代码来源:texp10.c


示例6: mpfr_mul3

static int
mpfr_mul3 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
    /* Old implementation */
    int sign_product, cc, inexact;
    mpfr_exp_t ax;
    mp_limb_t *tmp;
    mp_limb_t b1;
    mpfr_prec_t bq, cq;
    mp_size_t bn, cn, tn, k;
    MPFR_TMP_DECL(marker);

    /* deal with special cases */
    if (MPFR_ARE_SINGULAR(b,c))
    {
        if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c))
        {
            MPFR_SET_NAN(a);
            MPFR_RET_NAN;
        }
        sign_product = MPFR_MULT_SIGN( MPFR_SIGN(b) , MPFR_SIGN(c) );
        if (MPFR_IS_INF(b))
        {
            if (MPFR_IS_INF(c) || MPFR_NOTZERO(c))
            {
                MPFR_SET_SIGN(a,sign_product);
                MPFR_SET_INF(a);
                MPFR_RET(0); /* exact */
            }
            else
            {
                MPFR_SET_NAN(a);
                MPFR_RET_NAN;
            }
        }
        else if (MPFR_IS_INF(c))
        {
            if (MPFR_NOTZERO(b))
            {
                MPFR_SET_SIGN(a, sign_product);
                MPFR_SET_INF(a);
                MPFR_RET(0); /* exact */
            }
            else
            {
                MPFR_SET_NAN(a);
                MPFR_RET_NAN;
            }
        }
        else
        {
            MPFR_ASSERTD(MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c));
            MPFR_SET_SIGN(a, sign_product);
            MPFR_SET_ZERO(a);
            MPFR_RET(0); /* 0 * 0 is exact */
        }
    }
    sign_product = MPFR_MULT_SIGN( MPFR_SIGN(b) , MPFR_SIGN(c) );

    ax = MPFR_GET_EXP (b) + MPFR_GET_EXP (c);

    bq = MPFR_PREC(b);
    cq = MPFR_PREC(c);

    MPFR_ASSERTD(bq+cq > bq); /* PREC_MAX is /2 so no integer overflow */

    bn = (bq+GMP_NUMB_BITS-1)/GMP_NUMB_BITS; /* number of limbs of b */
    cn = (cq+GMP_NUMB_BITS-1)/GMP_NUMB_BITS; /* number of limbs of c */
    k = bn + cn; /* effective nb of limbs used by b*c (= tn or tn+1) below */
    tn = (bq + cq + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
    /* <= k, thus no int overflow */
    MPFR_ASSERTD(tn <= k);

    /* Check for no size_t overflow*/
    MPFR_ASSERTD((size_t) k <= ((size_t) -1) / BYTES_PER_MP_LIMB);
    MPFR_TMP_MARK(marker);
    tmp = (mp_limb_t *) MPFR_TMP_ALLOC((size_t) k * BYTES_PER_MP_LIMB);

    /* multiplies two mantissa in temporary allocated space */
    b1 = (MPFR_LIKELY(bn >= cn)) ?
         mpn_mul (tmp, MPFR_MANT(b), bn, MPFR_MANT(c), cn)
         : mpn_mul (tmp, MPFR_MANT(c), cn, MPFR_MANT(b), bn);

    /* now tmp[0]..tmp[k-1] contains the product of both mantissa,
       with tmp[k-1]>=2^(GMP_NUMB_BITS-2) */
    b1 >>= GMP_NUMB_BITS - 1; /* msb from the product */

    /* if the mantissas of b and c are uniformly distributed in ]1/2, 1],
       then their product is in ]1/4, 1/2] with probability 2*ln(2)-1 ~ 0.386
       and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */
    tmp += k - tn;
    if (MPFR_UNLIKELY(b1 == 0))
        mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */
    cc = mpfr_round_raw (MPFR_MANT (a), tmp, bq + cq,
                         MPFR_IS_NEG_SIGN(sign_product),
                         MPFR_PREC (a), rnd_mode, &inexact);

    /* cc = 1 ==> result is a power of two */
    if (MPFR_UNLIKELY(cc))
        MPFR_MANT(a)[MPFR_LIMB_SIZE(a)-1] = MPFR_LIMB_HIGHBIT;
//.........这里部分代码省略.........
开发者ID:gnooth,项目名称:xcl,代码行数:101,代码来源:mul.c


示例7: test_overflow2

static void
test_overflow2 (void)
{
  mpfr_t x, y, z, r;
  int i, inex, rnd, err = 0;

  mpfr_inits2 (8, x, y, z, r, (mpfr_ptr) 0);

  MPFR_SET_POS (x);
  mpfr_setmin (x, mpfr_get_emax ());  /* x = [email protected] */
  mpfr_set_si (y, -2, MPFR_RNDN);      /* y = -2 */
  /* The intermediate multiplication x * y will overflow. */

  for (i = -9; i <= 9; i++)
    RND_LOOP (rnd)
      {
        int inf, overflow;

        inf = rnd == MPFR_RNDN || rnd == MPFR_RNDD || rnd == MPFR_RNDA;
        overflow = inf || i <= 0;

        inex = mpfr_set_si_2exp (z, i, mpfr_get_emin (), MPFR_RNDN);
        MPFR_ASSERTN (inex == 0);

        mpfr_clear_flags ();
        /* One has: x * y = [email protected] exactly (but not representable). */
        inex = mpfr_fma (r, x, y, z, (mpfr_rnd_t) rnd);
        if (overflow ^ (mpfr_overflow_p () != 0))
          {
            printf ("Error in test_overflow2 (i = %d, %s): wrong overflow"
                    " flag (should be %d)\n", i,
                    mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), overflow);
            err = 1;
          }
        if (mpfr_nanflag_p ())
          {
            printf ("Error in test_overflow2 (i = %d, %s): NaN flag should"
                    " not be set\n", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
            err = 1;
          }
        if (mpfr_nan_p (r))
          {
            printf ("Error in test_overflow2 (i = %d, %s): got NaN\n",
                    i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
            err = 1;
          }
        else if (MPFR_SIGN (r) >= 0)
          {
            printf ("Error in test_overflow2 (i = %d, %s): wrong sign "
                    "(+ instead of -)\n", i,
                    mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
            err = 1;
          }
        else if (inf && ! mpfr_inf_p (r))
          {
            printf ("Error in test_overflow2 (i = %d, %s): expected -Inf,"
                    " got\n", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
            mpfr_dump (r);
            err = 1;
          }
        else if (!inf && (mpfr_inf_p (r) ||
                          (mpfr_nextbelow (r), ! mpfr_inf_p (r))))
          {
            printf ("Error in test_overflow2 (i = %d, %s): expected -MAX,"
                    " got\n", i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
            mpfr_dump (r);
            err = 1;
          }
        if (inf ? inex >= 0 : inex <= 0)
          {
            printf ("Error in test_overflow2 (i = %d, %s): wrong inexact"
                    " flag (got %d)\n", i,
                    mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), inex);
            err = 1;
          }

      }

  if (err)
    exit (1);
  mpfr_clears (x, y, z, r, (mpfr_ptr) 0);
}
开发者ID:bsmr-common-lisp,项目名称:xcl,代码行数:82,代码来源:tfma.c


示例8: main

int
main (void)
{
  mpfr_t x, y;
  float f, g, infp;
  int i;

  infp = (float) DBL_POS_INF;
  if (infp * 0.5 != infp)
    {
      fprintf (stderr, "Error, FLT_MAX + FLT_MAX does not yield INFP\n");
      fprintf (stderr, "(this is probably a compiler bug, please report)\n");
      exit (1);
    }

  tests_start_mpfr ();

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

#if !defined(MPFR_ERRDIVZERO)
  mpfr_set_nan (x);
  f = mpfr_get_flt (x, MPFR_RNDN);
  if (f == f)
    {
      printf ("Error for mpfr_get_flt(NaN)\n");
      exit (1);
    }
  mpfr_set_flt (x, f, MPFR_RNDN);
  if (mpfr_nan_p (x) == 0)
    {
      printf ("Error for mpfr_set_flt(NaN)\n");
      exit (1);
    }

  mpfr_set_inf (x, 1);
  f = mpfr_get_flt (x, MPFR_RNDN);
  mpfr_set_flt (x, f, MPFR_RNDN);
  if (mpfr_inf_p (x) == 0 || mpfr_sgn (x) < 0)
    {
      printf ("Error for mpfr_set_flt(mpfr_get_flt(+Inf)):\n");
      printf ("f=%f, expected -Inf\n", f);
      printf ("got "); mpfr_dump (x);
      exit (1);
    }

  mpfr_set_inf (x, -1);
  f = mpfr_get_flt (x, MPFR_RNDN);
  mpfr_set_flt (x, f, MPFR_RNDN);
  if (mpfr_inf_p (x) == 0 || mpfr_sgn (x) > 0)
    {
      printf ("Error for mpfr_set_flt(mpfr_get_flt(-Inf)):\n");
      printf ("f=%f, expected -Inf\n", f);
      printf ("got "); mpfr_dump (x);
      exit (1);
    }
#endif

  mpfr_set_ui (x, 0, MPFR_RNDN);
  f = mpfr_get_flt (x, MPFR_RNDN);
  mpfr_set_flt (x, f, MPFR_RNDN);
  if (mpfr_zero_p (x) == 0 || MPFR_SIGN (x) < 0)
    {
      printf ("Error for mpfr_set_flt(mpfr_get_flt(+0))\n");
      exit (1);
    }

  mpfr_set_ui (x, 0, MPFR_RNDN);
  mpfr_neg (x, x, MPFR_RNDN);
  f = mpfr_get_flt (x, MPFR_RNDN);
  mpfr_set_flt (x, f, MPFR_RNDN);
  if (mpfr_zero_p (x) == 0 || MPFR_SIGN (x) > 0)
    {
      printf ("Error for mpfr_set_flt(mpfr_get_flt(-0))\n");
      exit (1);
    }

  mpfr_set_ui (x, 17, MPFR_RNDN);
  f = mpfr_get_flt (x, MPFR_RNDN);
  mpfr_set_flt (x, f, MPFR_RNDN);
  if (mpfr_cmp_ui (x, 17) != 0)
    {
      printf ("Error for mpfr_set_flt(mpfr_get_flt(17))\n");
      printf ("expected 17\n");
      printf ("got      ");
      mpfr_dump (x);
      exit (1);
    }

  mpfr_set_si (x, -42, MPFR_RNDN);
  f = mpfr_get_flt (x, MPFR_RNDN);
  mpfr_set_flt (x, f, MPFR_RNDN);
  if (mpfr_cmp_si (x, -42) != 0)
    {
      printf ("Error for mpfr_set_flt(mpfr_get_flt(-42))\n");
      printf ("expected -42\n");
      printf ("got      ");
      mpfr_dump (x);
      exit (1);
    }
//.........这里部分代码省略.........
开发者ID:Kirija,项目名称:XPIR,代码行数:101,代码来源:tget_flt.c


示例9: mpfr_pow_si

int
mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd)
{
  MPFR_LOG_FUNC
    (("x[%Pu]=%.*Rg n=%ld rnd=%d",
      mpfr_get_prec (x), mpfr_log_prec, x, n, rnd),
     ("y[%Pu]=%.*Rg", mpfr_get_prec (y), mpfr_log_prec, y));

  if (n >= 0)
    return mpfr_pow_ui (y, x, n, rnd);
  else
    {
      if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
        {
          if (MPFR_IS_NAN (x))
            {
              MPFR_SET_NAN (y);
              MPFR_RET_NAN;
            }
          else
            {
              int positive = MPFR_IS_POS (x) || ((unsigned long) n & 1) == 0;
              if (MPFR_IS_INF (x))
                MPFR_SET_ZERO (y);
              else /* x is zero */
                {
                  MPFR_ASSERTD (MPFR_IS_ZERO (x));
                  MPFR_SET_INF (y);
                  mpfr_set_divby0 ();
                }
              if (positive)
                MPFR_SET_POS (y);
              else
                MPFR_SET_NEG (y);
              MPFR_RET (0);
            }
        }

      /* detect exact powers: x^(-n) is exact iff x is a power of 2 */
      if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0)
        {
          mpfr_exp_t expx = MPFR_EXP (x) - 1, expy;
          MPFR_ASSERTD (n < 0);
          /* Warning: n * expx may overflow!
           *
           * Some systems (apparently alpha-freebsd) abort with
           * LONG_MIN / 1, and LONG_MIN / -1 is undefined.
           * http://www.freebsd.org/cgi/query-pr.cgi?pr=72024
           *
           * Proof of the overflow checking. The expressions below are
           * assumed to be on the rational numbers, but the word "overflow"
           * still has its own meaning in the C context. / still denotes
           * the integer (truncated) division, and // denotes the exact
           * division.
           * - First, (__gmpfr_emin - 1) / n and (__gmpfr_emax - 1) / n
           *   cannot overflow due to the constraints on the exponents of
           *   MPFR numbers.
           * - If n = -1, then n * expx = - expx, which is representable
           *   because of the constraints on the exponents of MPFR numbers.
           * - If expx = 0, then n * expx = 0, which is representable.
           * - If n < -1 and expx > 0:
           *   + If expx > (__gmpfr_emin - 1) / n, then
           *           expx >= (__gmpfr_emin - 1) / n + 1
           *                > (__gmpfr_emin - 1) // n,
           *     and
           *           n * expx < __gmpfr_emin - 1,
           *     i.e.
           *           n * expx <= __gmpfr_emin - 2.
           *     This corresponds to an underflow, with a null result in
           *     the rounding-to-nearest mode.
           *   + If expx <= (__gmpfr_emin - 1) / n, then n * expx cannot
           *     overflow since 0 < expx <= (__gmpfr_emin - 1) / n and
           *           0 > n * expx >= n * ((__gmpfr_emin - 1) / n)
           *                        >= __gmpfr_emin - 1.
           * - If n < -1 and expx < 0:
           *   + If expx < (__gmpfr_emax - 1) / n, then
           *           expx <= (__gmpfr_emax - 1) / n - 1
           *                < (__gmpfr_emax - 1) // n,
           *     and
           *           n * expx > __gmpfr_emax - 1,
           *     i.e.
           *           n * expx >= __gmpfr_emax.
           *     This corresponds to an overflow (2^(n * expx) has an
           *     exponent > __gmpfr_emax).
           *   + If expx >= (__gmpfr_emax - 1) / n, then n * expx cannot
           *     overflow since 0 > expx >= (__gmpfr_emax - 1) / n and
           *           0 < n * expx <= n * ((__gmpfr_emax - 1) / n)
           *                        <= __gmpfr_emax - 1.
           * Note: one could use expx bounds based on MPFR_EXP_MIN and
           * MPFR_EXP_MAX instead of __gmpfr_emin and __gmpfr_emax. The
           * current bounds do not lead to noticeably slower code and
           * allow us to avoid a bug in Sun's compiler for Solaris/x86
           * (when optimizations are enabled); known affected versions:
           *   cc: Sun C 5.8 2005/10/13
           *   cc: Sun C 5.8 Patch 121016-02 2006/03/31
           *   cc: Sun C 5.8 Patch 121016-04 2006/10/18
           */
          expy =
            n != -1 && expx > 0 && expx > (__gmpfr_emin - 1) / n ?
            MPFR_EMIN_MIN - 2 /* Underflow */ :
//.........这里部分代码省略.........
开发者ID:Kirija,项目名称:XPIR,代码行数:101,代码来源:pow_si.c


示例10: mpfr_sin

int
mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t c, xr;
  mpfr_srcptr xx;
  mpfr_exp_t expx, err;
  mpfr_prec_t precy, m;
  int inexact, sign, reduce;
  MPFR_ZIV_DECL (loop);
  MPFR_SAVE_EXPO_DECL (expo);

  MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd_mode),
                  ("y[%#R]=%R inexact=%d", y, y, inexact));

  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
        {
          MPFR_SET_NAN (y);
          MPFR_RET_NAN;

        }
      else /* x is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
          MPFR_SET_ZERO (y);
          MPFR_SET_SAME_SIGN (y, x);
          MPFR_RET (0);
        }
    }

  /* sin(x) = x - x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2 * MPFR_GET_EXP (x), 2, 0,
                                    rnd_mode, {});

  MPFR_SAVE_EXPO_MARK (expo);

  /* Compute initial precision */
  precy = MPFR_PREC (y);

  if (precy >= MPFR_SINCOS_THRESHOLD)
    return mpfr_sin_fast (y, x, rnd_mode);

  m = precy + MPFR_INT_CEIL_LOG2 (precy) + 13;
  expx = MPFR_GET_EXP (x);

  mpfr_init (c);
  mpfr_init (xr);

  MPFR_ZIV_INIT (loop, m);
  for (;;)
    {
      /* first perform argument reduction modulo 2*Pi (if needed),
         also helps to determine the sign of sin(x) */
      if (expx >= 2) /* If Pi < x < 4, we need to reduce too, to determine
                        the sign of sin(x). For 2 <= |x| < Pi, we could avoid
                        the reduction. */
        {
          reduce = 1;
          /* As expx + m - 1 will silently be converted into mpfr_prec_t
             in the mpfr_set_prec call, the assert below may be useful to
             avoid undefined behavior. */
          MPFR_ASSERTN (expx + m - 1 <= MPFR_PREC_MAX);
          mpfr_set_prec (c, expx + m - 1);
          mpfr_set_prec (xr, m);
          mpfr_const_pi (c, MPFR_RNDN);
          mpfr_mul_2ui (c, c, 1, MPFR_RNDN);
          mpfr_remainder (xr, x, c, MPFR_RNDN);
          /* The analysis is similar to that of cos.c:
             |xr - x - 2kPi| <= 2^(2-m). Thus we can decide the sign
             of sin(x) if xr is at distance at least 2^(2-m) of both
             0 and +/-Pi. */
          mpfr_div_2ui (c, c, 1, MPFR_RNDN);
          /* Since c approximates Pi with an error <= 2^(2-expx-m) <= 2^(-m),
             it suffices to check that c - |xr| >= 2^(2-m). */
          if (MPFR_SIGN (xr) > 0)
            mpfr_sub (c, c, xr, MPFR_RNDZ);
          else
            mpfr_add (c, c, xr, MPFR_RNDZ);
          if (MPFR_IS_ZERO(xr)
              || MPFR_EXP(xr) < (mpfr_exp_t) 3 - (mpfr_exp_t) m
              || MPFR_EXP(c) < (mpfr_exp_t) 3 - (mpfr_exp_t) m)
            goto ziv_next;

          /* |xr - x - 2kPi| <= 2^(2-m), thus |sin(xr) - sin(x)| <= 2^(2-m) */
          xx = xr;
        }
      else /* the input argument is already reduced */
        {
          reduce = 0;
          xx = x;
        }

      sign = MPFR_SIGN(xx);
      /* now that the argument is reduced, precision m is enough */
      mpfr_set_prec (c, m);
      mpfr_cos (c, xx, MPFR_RNDZ);    /* can't be exact */
      mpfr_nexttoinf (c);           /* now c = cos(x) rounded away */
      mpfr_mul (c, c, c, MPFR_RNDU); /* away */
      mpfr_ui_sub (c, 1, c, MPFR_RNDZ);
//.........这里部分代码省略.........
开发者ID:119,项目名称:aircam-openwrt,代码行数:101,代码来源:sin.c


示例11: main


//.........这里部分代码省略.........
        baseprec = 1 + (prec - 2 + logbase) / logbase;
      str = mpfr_get_str (NULL, &e, base, baseprec, x, rnd);
      mpfr_set_str (y, str, base, rnd);
      MPFR_EXP(y) += logbase * (e - strlen (str));
      if (mpfr_cmp (x, y))
        {
          printf ("mpfr_set_str o mpfr_get_str <> id for rnd_mode=%s\n",
                  mpfr_print_rnd_mode (rnd));
          printf ("x=");
          mpfr_print_binary (x);
          puts ("");
          printf ("s=%s, exp=%d, base=%d\n", str, (int) e, base);
          printf ("y=");
          mpfr_print_binary (y);
          puts ("");
          mpfr_clear (x);
          mpfr_clear (y);
          exit (1);
        }
      (*__gmp_free_func) (str, strlen (str) + 1);
    }

  for (i = 2; i <= 62; i++)
    {
      if (mpfr_set_str (x, "@[email protected](garbage)", i, MPFR_RNDN) != 0 ||
          !mpfr_nan_p(x))
        {
          printf ("mpfr_set_str failed on @[email protected](garbage)\n");
          exit (1);
        }

      /*
      if (mpfr_set_str (x, "@[email protected]", i, MPFR_RNDN) != 0 ||
          !mpfr_inf_p(x) || MPFR_SIGN(x) < 0)
        {
          printf ("mpfr_set_str failed on @[email protected]\n");
          exit (1);
        }

      if (mpfr_set_str (x, "[email protected]@garbage", i, MPFR_RNDN) != 0 ||
          !mpfr_inf_p(x) || MPFR_SIGN(x) > 0)
        {
          printf ("mpfr_set_str failed on [email protected]@garbage\n");
          exit (1);
        }

      if (mpfr_set_str (x, "[email protected]@garbage", i, MPFR_RNDN) != 0 ||
          !mpfr_inf_p(x) || MPFR_SIGN(x) < 0)
        {
          printf ("mpfr_set_str failed on [email protected]@garbage\n");
          exit (1);
        }
      */

      if (i > 16)
        continue;

      if (mpfr_set_str (x, "NaN", i, MPFR_RNDN) != 0 ||
          !mpfr_nan_p(x))
        {
          printf ("mpfr_set_str failed on NaN\n");
          exit (1);
        }

      if (mpfr_set_str (x, "Inf", i, MPFR_RNDN) != 0 ||
          !mpfr_inf_p(x) || MPFR_SIGN(x) < 0)
开发者ID:Kirija,项目名称:XPIR,代码行数:67,代码来源:tset_str.c


示例12: sign

 /// Returns sign of the number
 inline int sign() const { return MPFR_SIGN(val); }
开发者ID:dlevin256,项目名称:kfr,代码行数:2,代码来源:mpfrplus.hpp


示例13: main

int
main (void)
{
  mpfr_t x, y, z;
  int i, j, k;

  tests_start_mpfr ();

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

  for (i = 0; i <= 1; i++)
    for (j = 0; j <= 1; j++)
      for (k = 0; k <= 5; k++)
        {
          mpfr_set_nan (x);
          i ? MPFR_SET_NEG (x) : MPFR_SET_POS (x);
          mpfr_set_nan (y);
          j ? MPFR_SET_NEG (y) : MPFR_SET_POS (y);
          copysign_variant (z, x, y, GMP_RNDN, k);
          if (MPFR_SIGN (z) != MPFR_SIGN (y) || !mpfr_nanflag_p ())
            {
              printf ("Error in mpfr_copysign (%cNaN, %cNaN)\n",
                      i ? '-' : '+', j ? '-' : '+');
              exit (1);
            }

          mpfr_set_si (x, i ? -1250 : 1250, GMP_RNDN);
          mpfr_set_nan (y);
          j ? MPFR_SET_NEG (y) : MPFR_SET_POS (y);
          copysign_variant (z, x, y, GMP_RNDN, k);
          if (i != j)
            mpfr_neg (x, x, GMP_RNDN);
          if (! mpfr_equal_p (z, x) || mpfr_nanflag_p ())
            {
              printf ("Error in mpfr_copysign (%c1250, %cNaN)\n",
                      i ? '-' : '+', j ? '-' : '+');
              exit (1);
            }

          mpfr_set_si (x, i ? -1250 : 1250, GMP_RNDN);
          mpfr_set_si (y, j ? -1717 : 1717, GMP_RNDN);
          copysign_variant (z, x, y, GMP_RNDN, k);
          if (i != j)
            mpfr_neg (x, x, GMP_RNDN);
          if (! mpfr_equal_p (z, x) || mpfr_nanflag_p ())
            {
              printf ("Error in mpfr_copysign (%c1250, %c1717)\n",
                      i ? '-' : '+', j ? '-' : '+');
              exit (1);
            }
        }

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

  tests_end_mpfr ();
  return 0;
}
开发者ID:mmanley,项目名称:Antares,代码行数:61,代码来源:tcopysign.c


示例14: mpfr_sub1sp


//.........这里部分代码省略.........
                {
                  MPFR_ASSERTD( k > 0 );
                  limb = ap[--k];
                }
              while (limb == 0);
              MPFR_ASSERTD(limb != 0);
              count_leading_zeros(cnt, limb);
              k++;
              len = n - k; /* Number of last limb */
              MPFR_ASSERTD(k >= 0);
              if (MPFR_LIKELY(cnt))
                mpn_lshift(ap+len, ap, k, cnt); /* Normalize the High Limb*/
              else
                {
                  /* Must use DECR since src and dest may overlap & dest>=src*/
                  MPN_COPY_DECR(ap+len, ap, k);
                }
              MPN_ZERO(ap, len); /* Zeroing the last limbs */
              bx -= cnt + len*GMP_NUMB_BITS; /* Update Expo */
              /* Last limb should be ok */
              MPFR_ASSERTD(!(ap[len]&MPFR_LIMB_MASK((unsigned int) (-p)
                                                    % GMP_NUMB_BITS)));
            }
          /* Check expo underflow */
          if (MPFR_UNLIKELY(bx < __gmpfr_emin))
            {
              MPFR_TMP_FREE(marker);
              /* inexact=0 */
              DEBUG( printf("(D==0 Underflow)\n") );
              if (rnd_mode == MPFR_RNDN &&
                  (bx < __gmpfr_emin - 1 ||
                   (/*inexact >= 0 &&*/ mpfr_powerof2_raw (a))))
                rnd_mode = MPFR_RNDZ;
              return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
            }
          MPFR_SET_EXP (a, bx);
          /* No rounding is necessary since the result is exact */
          MPFR_ASSERTD(ap[n-1] > ~ap[n-1]);
          MPFR_TMP_FREE(marker);
          return 0;
        }
      else /* if (d == 1) */
        {
          /* | <-- b -->
             |  <-- c --> */
          mp_limb_t c0, mask;
          mp_size_t k;
          MPFR_UNSIGNED_MINUS_MODULO(sh, p);
          /* If we lose at least one bit, compute 2*b-c (Exact)
           * else compute b-c/2 */
          bp = MPFR_MANT(b);
          cp = MPFR_MANT(c);
          k = n-1;
          limb = bp[k] - cp[k]/2;
          if (limb > MPFR_LIMB_HIGHBIT)
            {
              /* We can't lose precision: compute b-c/2 */
              /* Shift c in the allocated temporary block */
            SubD1NoLose:
              c0 = cp[0] & (MPFR_LIMB_ONE<<sh);
              cp = MPFR_TMP_LIMBS_ALLOC (n);
              mpn_rshift(cp, MPFR_MANT(c), n, 1);
              if (MPFR_LIKELY(c0 == 0))
                {
                  /* Result is exact: no need of rounding! */
                  ap = MPFR_MANT(a);
开发者ID:texlive,项目名称:texlive-source,代码行数:67,代码来源:sub1sp.c


示例15: check_for_zero

static void
check_for_zero (void)
{
  /* Check that 0 is unsigned! */
  mpq_t q;
  mpz_t z;
  mpfr_t x;
  int r;
  mpfr_sign_t i;

  mpfr_init (x);
  mpz_init (z);
  mpq_init (q);

  mpz_set_ui (z, 0);
  mpq_set_ui (q, 0, 1);

  MPFR_SET_ZERO (x);
  RND_LOOP (r)
    {
      for (i = MPFR_SIGN_NEG ; i <= MPFR_SIGN_POS ;
           i+=MPFR_SIGN_POS-MPFR_SIGN_NEG)
        {
          MPFR_SET_SIGN(x, i);
          mpfr_add_z (x, x, z, (mpfr_rnd_t) r);
          if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
            {
              printf("GMP Zero errors for add_z & rnd=%s & s=%d\n",
                     mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
              mpfr_dump (x);
              exit (1);
            }
          mpfr_sub_z (x, x, z, (mpfr_rnd_t) r);
          if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
            {
              printf("GMP Zero errors for sub_z & rnd=%s & s=%d\n",
                     mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
              mpfr_dump (x);
              exit (1);
            }
          mpfr_mul_z (x, x, z, (mpfr_rnd_t) r);
          if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
            {
              printf("GMP Zero errors for mul_z & rnd=%s & s=%d\n",
                     mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
              mpfr_dump (x);
              exit (1);
            }
          mpfr_add_q (x, x, q, (mpfr_rnd_t) r);
          if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
            {
              printf("GMP Zero errors for add_q & rnd=%s & s=%d\n",
                     mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
              mpfr_dump (x);
              exit (1);
            }
          mpfr_sub_q (x, x, q, (mpfr_rnd_t) r);
          if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
            {
              printf("GMP Zero errors for sub_q & rnd=%s & s=%d\n",
                     mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
              mpfr_dump (x);
              exit (1);
             }
        }
    }

  mpq_clear (q);
  mpz_clear (z);
  mpfr_clear (x);
}
开发者ID:Canar,项目名称:mpfr,代码行数:71,代码来源:tgmpop.c


示例16: mpfr_add


鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
C++ MPFR_ZIV_DECL函数代码示例发布时间:2022-05-30
下一篇:
C++ MPFR_SET_ZERO函数代码示例发布时间:2022-05-30
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap