当前位置: 首页>>代码示例>>C++>>正文


C++ MPFR_CAN_ROUND函数代码示例

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


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

示例1: tgeneric_cci

static void
tgeneric_cci (mpc_function *function, mpc_ptr op1, int op2,
              mpc_ptr rop, mpc_ptr rop4, mpc_ptr rop4rnd, mpc_rnd_t rnd)
{
  known_signs_t ks = {1, 1};

  function->pointer.CCI (rop4, op1, op2, rnd);
  function->pointer.CCI (rop, op1, op2, rnd);
  if (MPFR_CAN_ROUND (mpc_realref (rop4), 1, MPC_PREC_RE (rop),
                      MPC_RND_RE (rnd))
      && MPFR_CAN_ROUND (mpc_imagref (rop4), 1, MPC_PREC_IM (rop),
                         MPC_RND_IM (rnd)))
    mpc_set (rop4rnd, rop4, rnd);
  else
    return;

  if (same_mpc_value (rop, rop4rnd, ks))
    return;

  printf ("Rounding in %s might be incorrect for\n", function->name);
  MPC_OUT (op1);
  printf ("op2=%d\n", op2);
  printf ("with rounding mode (%s, %s)",
          mpfr_print_rnd_mode (MPC_RND_RE (rnd)),
          mpfr_print_rnd_mode (MPC_RND_IM (rnd)));

  printf ("\n%s                     gives ", function->name);
  MPC_OUT (rop);
  printf ("%s quadruple precision gives ", function->name);
  MPC_OUT (rop4);
  printf ("and is rounded to                  ");
  MPC_OUT (rop4rnd);

  exit (1);
}
开发者ID:Distrotech,项目名称:mpc,代码行数:35,代码来源:tgeneric.c

示例2: tgeneric_cc_c

static void
tgeneric_cc_c (mpc_function *function, mpc_ptr op, mpc_ptr rop1, mpc_ptr rop2,
   mpc_ptr rop14, mpc_ptr rop24, mpc_ptr rop14rnd, mpc_ptr rop24rnd,
   mpc_rnd_t rnd1, mpc_rnd_t rnd2)
{
   /* same as the previous function, but for mpc functions computing two
      results from one argument                                          */
   known_signs_t ks = {1, 1};

   function->pointer.CC_C (rop14, rop24, op, rnd1, rnd2);
   function->pointer.CC_C (rop1,  rop2,  op, rnd1, rnd2);

   if (   MPFR_CAN_ROUND (mpc_realref (rop14), 1, MPC_PREC_RE (rop1),
                          MPC_RND_RE (rnd1))
       && MPFR_CAN_ROUND (mpc_imagref (rop14), 1, MPC_PREC_IM (rop1),
                          MPC_RND_IM (rnd1))
       && MPFR_CAN_ROUND (mpc_realref (rop24), 1, MPC_PREC_RE (rop2),
                          MPC_RND_RE (rnd2))
       && MPFR_CAN_ROUND (mpc_imagref (rop24), 1, MPC_PREC_IM (rop2),
                          MPC_RND_IM (rnd2))) {
     mpc_set (rop14rnd, rop14, rnd1);
     mpc_set (rop24rnd, rop24, rnd2);
   }
   else
     return;

   if (!same_mpc_value (rop1, rop14rnd, ks)) {
      /* rounding failed for first result */
      printf ("Rounding might be incorrect for the first result of %s at\n", function->name);
      MPC_OUT (op);
      printf ("with rounding mode (%s, %s)",
          mpfr_print_rnd_mode (MPC_RND_RE (rnd1)),
          mpfr_print_rnd_mode (MPC_RND_IM (rnd1)));
      printf ("\n%s                     gives ", function->name);
      MPC_OUT (rop1);
      printf ("%s quadruple precision gives ", function->name);
      MPC_OUT (rop14);
      printf ("and is rounded to                  ");
      MPC_OUT (rop14rnd);
      exit (1);
   }
   else if (!same_mpc_value (rop2, rop24rnd, ks)) {
      /* rounding failed for second result */
      printf ("Rounding might be incorrect for the second result of %s at\n", function->name);
      MPC_OUT (op);
      printf ("with rounding mode (%s, %s)",
          mpfr_print_rnd_mode (MPC_RND_RE (rnd2)),
          mpfr_print_rnd_mode (MPC_RND_IM (rnd2)));
      printf ("\n%s                     gives ", function->name);
      MPC_OUT (rop2);
      printf ("%s quadruple precision gives ", function->name);
      MPC_OUT (rop24);
      printf ("and is rounded to                  ");
      MPC_OUT (rop24rnd);
      exit (1);
   }
}
开发者ID:Distrotech,项目名称:mpc,代码行数:57,代码来源:tgeneric.c

示例3: Li2

/* try asymptotic expansion when x is large and negative:
   Li2(x) = -log(-x)^2/2 - Pi^2/6 - 1/x + O(1/x^2).
   More precisely for x <= -2 we have for g(x) = -log(-x)^2/2 - Pi^2/6:
   |Li2(x) - g(x)| <= 1/|x|.
   Assumes x <= -7, which ensures |log(-x)^2/2| >= Pi^2/6, and g(x) <= -3.5.
   Return 0 if asymptotic expansion failed (unable to round), otherwise
   returns correct ternary value.
*/
static int
mpfr_li2_asympt_neg (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t g, h;
  mp_prec_t w = MPFR_PREC (y) + 20;
  int inex = 0;

  MPFR_ASSERTN (mpfr_cmp_si (x, -7) <= 0);

  mpfr_init2 (g, w);
  mpfr_init2 (h, w);
  mpfr_neg (g, x, GMP_RNDN);
  mpfr_log (g, g, GMP_RNDN);    /* rel. error <= |(1 + theta) - 1| */
  mpfr_sqr (g, g, GMP_RNDN);    /* rel. error <= |(1 + theta)^3 - 1| <= 2^(2-w) */
  mpfr_div_2ui (g, g, 1, GMP_RNDN);     /* rel. error <= 2^(2-w) */
  mpfr_const_pi (h, GMP_RNDN);  /* error <= 2^(1-w) */
  mpfr_sqr (h, h, GMP_RNDN);    /* rel. error <= 2^(2-w) */
  mpfr_div_ui (h, h, 6, GMP_RNDN);      /* rel. error <= |(1 + theta)^4 - 1|
                                           <= 5 * 2^(-w) */
  MPFR_ASSERTN (MPFR_EXP (g) >= MPFR_EXP (h));
  mpfr_add (g, g, h, GMP_RNDN); /* err <= ulp(g)/2 + g*2^(2-w) + g*5*2^(-w)
                                   <= ulp(g) * (1/2 + 4 + 5) < 10 ulp(g).

                                   If in addition |1/x| <= 4 ulp(g), then the
                                   total error is bounded by 16 ulp(g). */
  if ((MPFR_EXP (x) >= (mp_exp_t) (w - 2) - MPFR_EXP (g)) &&
      MPFR_CAN_ROUND (g, w - 4, MPFR_PREC (y), rnd_mode))
    inex = mpfr_neg (y, g, rnd_mode);

  mpfr_clear (g);
  mpfr_clear (h);

  return inex;
}
开发者ID:Scorpiion,项目名称:Renux_cross_gcc,代码行数:42,代码来源:li2.c

示例4: tgeneric_fc

static void
tgeneric_fc (mpc_function *function, mpc_ptr op, mpfr_ptr rop,
             mpfr_ptr rop4, mpfr_ptr rop4rnd, mpfr_rnd_t rnd)
{
  function->pointer.FC (rop4, op, rnd);
  function->pointer.FC (rop, op, rnd);
  if (MPFR_CAN_ROUND (rop4, 1, mpfr_get_prec (rop), rnd))
    mpfr_set (rop4rnd, rop4, rnd);
  else
    return;

  if (same_mpfr_value (rop, rop4rnd, 1))
    return;

  printf ("Rounding in %s might be incorrect for\n", function->name);
  MPC_OUT (op);
  printf ("with rounding mode %s", mpfr_print_rnd_mode (rnd));

  printf ("\n%s                     gives ", function->name);
  MPFR_OUT (rop);
  printf ("%s quadruple precision gives ", function->name);
  MPFR_OUT (rop4);
  printf ("and is rounded to                  ");
  MPFR_OUT (rop4rnd);

  exit (1);
}
开发者ID:Distrotech,项目名称:mpc,代码行数:27,代码来源:tgeneric.c

示例5: tgeneric_cc

/* functions with one input, one output */
static void
tgeneric_cc (mpc_function *function, mpc_ptr op, mpc_ptr rop,
             mpc_ptr rop4, mpc_ptr rop4rnd, mpc_rnd_t rnd)
{
  known_signs_t ks = {1, 1};

  /* We compute the result with four times the precision and check whether the
     rounding is correct. Error reports in this part of the algorithm might
     still be wrong, though, since there are two consecutive roundings (but we
     try to avoid them).  */
  function->pointer.CC (rop4, op, rnd);
  function->pointer.CC (rop, op, rnd);

  /* can't use the mpfr_can_round function when argument is singular,
     use a custom macro instead. */
  if (MPFR_CAN_ROUND (mpc_realref (rop4), 1, MPC_PREC_RE (rop),
                      MPC_RND_RE (rnd))
      && MPFR_CAN_ROUND (mpc_imagref (rop4), 1, MPC_PREC_IM (rop),
                         MPC_RND_IM (rnd)))
    mpc_set (rop4rnd, rop4, rnd);
  else
    /* avoid double rounding error */
    return;

  if (same_mpc_value (rop, rop4rnd, ks))
    return;

  /* rounding failed */
  printf ("Rounding in %s might be incorrect for\n", function->name);
  MPC_OUT (op);

  printf ("with rounding mode (%s, %s)",
          mpfr_print_rnd_mode (MPC_RND_RE (rnd)),
          mpfr_print_rnd_mode (MPC_RND_IM (rnd)));

  printf ("\n%s                     gives ", function->name);
  MPC_OUT (rop);
  printf ("%s quadruple precision gives ", function->name);
  MPC_OUT (rop4);
  printf ("and is rounded to                  ");
  MPC_OUT (rop4rnd);

  exit (1);
}
开发者ID:Distrotech,项目名称:mpc,代码行数:45,代码来源:tgeneric.c

示例6: mpfr_const_euler_internal

int
mpfr_const_euler_internal (mpfr_t x, mpfr_rnd_t rnd)
{
  mpfr_prec_t prec = MPFR_PREC(x), m, log2m;
  mpfr_t y, z;
  unsigned long n;
  int inexact;
  MPFR_ZIV_DECL (loop);

  log2m = MPFR_INT_CEIL_LOG2 (prec);
  m = prec + 2 * log2m + 23;

  mpfr_init2 (y, m);
  mpfr_init2 (z, m);

  MPFR_ZIV_INIT (loop, m);
  for (;;)
    {
      mpfr_exp_t exp_S, err;
      /* since prec >= 1, we have m >= 24 here, which ensures n >= 9 below */
      n = 1 + (unsigned long) ((double) m * LOG2 / 2.0);
      MPFR_ASSERTD (n >= 9);
      mpfr_const_euler_S2 (y, n); /* error <= 3 ulps */
      exp_S = MPFR_EXP(y);
      mpfr_set_ui (z, n, MPFR_RNDN);
      mpfr_log (z, z, MPFR_RNDD); /* error <= 1 ulp */
      mpfr_sub (y, y, z, MPFR_RNDN); /* S'(n) - log(n) */
      /* the error is less than 1/2 + 3*2^(exp_S-EXP(y)) + 2^(EXP(z)-EXP(y))
         <= 1/2 + 2^(exp_S+2-EXP(y)) + 2^(EXP(z)-EXP(y))
         <= 1/2 + 2^(1+MAX(exp_S+2,EXP(z))-EXP(y)) */
      err = 1 + MAX(exp_S + 2, MPFR_EXP(z)) - MPFR_EXP(y);
      err = (err >= -1) ? err + 1 : 0; /* error <= 2^err ulp(y) */
      exp_S = MPFR_EXP(y);
      mpfr_const_euler_R (z, n); /* err <= ulp(1/2) = 2^(-m) */
      mpfr_sub (y, y, z, MPFR_RNDN);
      /* err <= 1/2 ulp(y) + 2^(-m) + 2^(err + exp_S - EXP(y)) ulp(y).
         Since the result is between 0.5 and 1, ulp(y) = 2^(-m).
         So we get 3/2*ulp(y) + 2^(err + exp_S - EXP(y)) ulp(y).
         3/2 + 2^e <= 2^(e+1) for e>=1, and <= 2^2 otherwise */
      err = err + exp_S - MPFR_EXP(y);
      err = (err >= 1) ? err + 1 : 2;
      if (MPFR_LIKELY (MPFR_CAN_ROUND (y, m - err, prec, rnd)))
        break;
      MPFR_ZIV_NEXT (loop, m);
      mpfr_set_prec (y, m);
      mpfr_set_prec (z, m);
    }
  MPFR_ZIV_FREE (loop);

  inexact = mpfr_set (x, y, rnd);

  mpfr_clear (y);
  mpfr_clear (z);

  return inexact; /* always inexact */
}
开发者ID:119,项目名称:aircam-openwrt,代码行数:56,代码来源:const_euler.c

示例7: mpfr_atanh

int
mpfr_atanh (mpfr_ptr y, mpfr_srcptr xt , mpfr_rnd_t rnd_mode)
{
  int inexact;
  mpfr_t x, t, te;
  mpfr_prec_t Nx, Ny, Nt;
  mpfr_exp_t err;
  MPFR_ZIV_DECL (loop);
  MPFR_SAVE_EXPO_DECL (expo);

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

  /* Special cases */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
    {
      /* atanh(NaN) = NaN, and atanh(+/-Inf) = NaN since tanh gives a result
         between -1 and 1 */
      if (MPFR_IS_NAN (xt) || MPFR_IS_INF (xt))
        {
          MPFR_SET_NAN (y);
          MPFR_RET_NAN;
        }
      else /* necessarily xt is 0 */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (xt));
          MPFR_SET_ZERO (y);   /* atanh(0) = 0 */
          MPFR_SET_SAME_SIGN (y,xt);
          MPFR_RET (0);
        }
    }

  /* atanh (x) = NaN as soon as |x| > 1, and arctanh(+/-1) = +/-Inf */
  if (MPFR_UNLIKELY (MPFR_GET_EXP (xt) > 0))
    {
      if (MPFR_GET_EXP (xt) == 1 && mpfr_powerof2_raw (xt))
        {
          MPFR_SET_INF (y);
          MPFR_SET_SAME_SIGN (y, xt);
          mpfr_set_divby0 ();
          MPFR_RET (0);
        }
      MPFR_SET_NAN (y);
      MPFR_RET_NAN;
    }

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

  MPFR_SAVE_EXPO_MARK (expo);

  /* Compute initial precision */
  Nx = MPFR_PREC (xt);
  MPFR_TMP_INIT_ABS (x, xt);
  Ny = MPFR_PREC (y);
  Nt = MAX (Nx, Ny);
  /* the optimal number of bits : see algorithms.ps */
  Nt = Nt + MPFR_INT_CEIL_LOG2 (Nt) + 4;

  /* initialise of intermediary variable */
  mpfr_init2 (t, Nt);
  mpfr_init2 (te, Nt);

  /* First computation of cosh */
  MPFR_ZIV_INIT (loop, Nt);
  for (;;)
    {
      /* compute atanh */
      mpfr_ui_sub (te, 1, x, MPFR_RNDU);   /* (1-xt)*/
      mpfr_add_ui (t,  x, 1, MPFR_RNDD);   /* (xt+1)*/
      mpfr_div (t, t, te, MPFR_RNDN);      /* (1+xt)/(1-xt)*/
      mpfr_log (t, t, MPFR_RNDN);          /* ln((1+xt)/(1-xt))*/
      mpfr_div_2ui (t, t, 1, MPFR_RNDN);   /* (1/2)*ln((1+xt)/(1-xt))*/

      /* error estimate: see algorithms.tex */
      /* FIXME: this does not correspond to the value in algorithms.tex!!! */
      /* err=Nt-__gmpfr_ceil_log2(1+5*pow(2,1-MPFR_EXP(t)));*/
      err = Nt - (MAX (4 - MPFR_GET_EXP (t), 0) + 1);

      if (MPFR_LIKELY (MPFR_IS_ZERO (t)
                       || MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
        break;

      /* reactualisation of the precision */
      MPFR_ZIV_NEXT (loop, Nt);
      mpfr_set_prec (t, Nt);
      mpfr_set_prec (te, Nt);
    }
  MPFR_ZIV_FREE (loop);

  inexact = mpfr_set4 (y, t, rnd_mode, MPFR_SIGN (xt));

  mpfr_clear(t);
  mpfr_clear(te);

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

示例8: mpfr_exp2


//.........这里部分代码省略.........
            MPFR_SET_INF (y);
          else
            MPFR_SET_ZERO (y);
          MPFR_SET_POS (y);
          MPFR_RET (0);
        }
      else /* 2^0 = 1 */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO(x));
          return mpfr_set_ui (y, 1, rnd_mode);
        }
    }

  /* since the smallest representable non-zero float is 1/2*2^__gmpfr_emin,
     if x < __gmpfr_emin - 1, the result is either 1/2*2^__gmpfr_emin or 0 */
  MPFR_ASSERTN (MPFR_EMIN_MIN >= LONG_MIN + 2);
  if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emin - 1) < 0))
    {
      mpfr_rnd_t rnd2 = rnd_mode;
      /* in round to nearest mode, round to zero when x <= __gmpfr_emin-2 */
      if (rnd_mode == MPFR_RNDN &&
          mpfr_cmp_si_2exp (x, __gmpfr_emin - 2, 0) <= 0)
        rnd2 = MPFR_RNDZ;
      return mpfr_underflow (y, rnd2, 1);
    }

  MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
  if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emax) >= 0))
    return mpfr_overflow (y, rnd_mode, 1);

  /* We now know that emin - 1 <= x < emax. */

  MPFR_SAVE_EXPO_MARK (expo);

  /* 2^x = 1 + x*log(2) + O(x^2) for x near zero, and for |x| <= 1 we have
     |2^x - 1| <= x < 2^EXP(x). If x > 0 we must round away from 0 (dir=1);
     if x < 0 we must round toward 0 (dir=0). */
  MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, __gmpfr_one, - MPFR_GET_EXP (x), 0,
                                    MPFR_IS_POS (x), rnd_mode, expo, {});

  xint = mpfr_get_si (x, MPFR_RNDZ);
  mpfr_init2 (xfrac, MPFR_PREC (x));
  mpfr_sub_si (xfrac, x, xint, MPFR_RNDN); /* exact */

  if (MPFR_IS_ZERO (xfrac))
    {
      mpfr_set_ui (y, 1, MPFR_RNDN);
      inexact = 0;
    }
  else
    {
      /* Declaration of the intermediary variable */
      mpfr_t t;

      /* Declaration of the size variable */
      mpfr_prec_t Ny = MPFR_PREC(y);              /* target precision */
      mpfr_prec_t Nt;                             /* working precision */
      mpfr_exp_t err;                             /* error */
      MPFR_ZIV_DECL (loop);

      /* compute the precision of intermediary variable */
      /* the optimal number of bits : see algorithms.tex */
      Nt = Ny + 5 + MPFR_INT_CEIL_LOG2 (Ny);

      /* initialize of intermediary variable */
      mpfr_init2 (t, Nt);

      /* First computation */
      MPFR_ZIV_INIT (loop, Nt);
      for (;;)
        {
          /* compute exp(x*ln(2))*/
          mpfr_const_log2 (t, MPFR_RNDU);       /* ln(2) */
          mpfr_mul (t, xfrac, t, MPFR_RNDU);    /* xfrac * ln(2) */
          err = Nt - (MPFR_GET_EXP (t) + 2);   /* Estimate of the error */
          mpfr_exp (t, t, MPFR_RNDN);           /* exp(xfrac * ln(2)) */

          if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
            break;

          /* Actualisation of the precision */
          MPFR_ZIV_NEXT (loop, Nt);
          mpfr_set_prec (t, Nt);
        }
      MPFR_ZIV_FREE (loop);

      inexact = mpfr_set (y, t, rnd_mode);

      mpfr_clear (t);
    }

  mpfr_clear (xfrac);
  MPFR_CLEAR_FLAGS ();
  mpfr_mul_2si (y, y, xint, MPFR_RNDN); /* exact or overflow */
  /* Note: We can have an overflow only when t was rounded up to 2. */
  MPFR_ASSERTD (MPFR_IS_PURE_FP (y) || inexact > 0);
  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (y, inexact, rnd_mode);
}
开发者ID:Canar,项目名称:mpfr,代码行数:101,代码来源:exp2.c

示例9: 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:Kirija,项目名称:XPIR,代码行数:101,代码来源:jyn_asympt.c

示例10: mpfr_log1p


//.........这里部分代码省略.........
            {
              MPFR_SET_NAN (y);
              MPFR_RET_NAN;
            }
        }
      else /* x is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
          MPFR_SET_ZERO (y);   /* log1p(+/- 0) = +/- 0 */
          MPFR_SET_SAME_SIGN (y, x);
          MPFR_RET (0);
        }
    }

  ex = MPFR_GET_EXP (x);
  if (ex < 0)  /* -0.5 < x < 0.5 */
    {
      /* For x > 0,    abs(log(1+x)-x) < x^2/2.
         For x > -0.5, abs(log(1+x)-x) < x^2. */
      if (MPFR_IS_POS (x))
        MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex - 1, 0, 0, rnd_mode, {});
      else
        MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 0, 1, rnd_mode, {});
    }

  comp = mpfr_cmp_si (x, -1);
  /* log1p(x) is undefined for x < -1 */
  if (MPFR_UNLIKELY(comp <= 0))
    {
      if (comp == 0)
        /* x=0: log1p(-1)=-inf (division by zero) */
        {
          MPFR_SET_INF (y);
          MPFR_SET_NEG (y);
          MPFR_RET (0);
        }
      MPFR_SET_NAN (y);
      MPFR_RET_NAN;
    }

  MPFR_SAVE_EXPO_MARK (expo);

  /* General case */
  {
    /* Declaration of the intermediary variable */
    mpfr_t t;
    /* Declaration of the size variable */
    mp_prec_t Ny = MPFR_PREC(y);             /* target precision */
    mp_prec_t Nt;                            /* working precision */
    mp_exp_t err;                            /* error */
    MPFR_ZIV_DECL (loop);

    /* compute the precision of intermediary variable */
    /* the optimal number of bits : see algorithms.tex */
    Nt = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 6;

    /* if |x| is smaller than 2^(-e), we will loose about e bits
       in log(1+x) */
    if (MPFR_EXP(x) < 0)
      Nt += -MPFR_EXP(x);

    /* initialise of intermediary variable */
    mpfr_init2 (t, Nt);

    /* First computation of log1p */
    MPFR_ZIV_INIT (loop, Nt);
    for (;;)
      {
        /* compute log1p */
        inexact = mpfr_add_ui (t, x, 1, GMP_RNDN);      /* 1+x */
        /* if inexact = 0, then t = x+1, and the result is simply log(t) */
        if (inexact == 0)
          {
            inexact = mpfr_log (y, t, rnd_mode);
            goto end;
          }
        mpfr_log (t, t, GMP_RNDN);        /* log(1+x) */

        /* the error is bounded by (1/2+2^(1-EXP(t))*ulp(t) (cf algorithms.tex)
           if EXP(t)>=2, then error <= ulp(t)
           if EXP(t)<=1, then error <= 2^(2-EXP(t))*ulp(t) */
        err = Nt - MAX (0, 2 - MPFR_GET_EXP (t));

        if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
          break;

        /* increase the precision */
        MPFR_ZIV_NEXT (loop, Nt);
        mpfr_set_prec (t, Nt);
      }
    inexact = mpfr_set (y, t, rnd_mode);

  end:
    MPFR_ZIV_FREE (loop);
    mpfr_clear (t);
  }

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (y, inexact, rnd_mode);
}
开发者ID:Scorpiion,项目名称:Renux_cross_gcc,代码行数:101,代码来源:log1p.c

示例11: mpfr_pow_si


//.........这里部分代码省略.........
           *     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).
           */
          expy =
            n != -1 && expx > 0 && expx > (__gmpfr_emin - 1) / n ?
            MPFR_EMIN_MIN - 2 /* Underflow */ :
            n != -1 && expx < 0 && expx < (__gmpfr_emax - 1) / n ?
            MPFR_EMAX_MAX /* Overflow */ : n * expx;
          return mpfr_set_si_2exp (y, n % 2 ? MPFR_INT_SIGN (x) : 1,
                                   expy, rnd);
        }

      /* General case */
      {
        /* Declaration of the intermediary variable */
        mpfr_t t;
        /* Declaration of the size variable */
        mp_prec_t Ny = MPFR_PREC (y);               /* target precision */
        mp_prec_t Nt;                              /* working precision */
        mp_exp_t  err;                             /* error */
        int inexact;
        unsigned long abs_n;
        MPFR_SAVE_EXPO_DECL (expo);
        MPFR_ZIV_DECL (loop);

        abs_n = - (unsigned long) n;

        /* compute the precision of intermediary variable */
        /* the optimal number of bits : see algorithms.tex */
        Nt = Ny + 3 + MPFR_INT_CEIL_LOG2 (Ny);

        MPFR_SAVE_EXPO_MARK (expo);

        /* initialise of intermediary   variable */
        mpfr_init2 (t, Nt);

        MPFR_ZIV_INIT (loop, Nt);
        for (;;)
          {
            /* compute 1/(x^n), with n > 0 */
            mpfr_pow_ui (t, x, abs_n, GMP_RNDN);
            mpfr_ui_div (t, 1, t, GMP_RNDN);
            /* FIXME: old code improved, but I think this is still incorrect. */
            if (MPFR_UNLIKELY (MPFR_IS_ZERO (t)))
              {
                MPFR_ZIV_FREE (loop);
                mpfr_clear (t);
                MPFR_SAVE_EXPO_FREE (expo);
                return mpfr_underflow (y, rnd == GMP_RNDN ? GMP_RNDZ : rnd,
                                       abs_n & 1 ? MPFR_SIGN (x) :
                                       MPFR_SIGN_POS);
              }
            if (MPFR_UNLIKELY (MPFR_IS_INF (t)))
              {
                MPFR_ZIV_FREE (loop);
                mpfr_clear (t);
                MPFR_SAVE_EXPO_FREE (expo);
                return mpfr_overflow (y, rnd, abs_n & 1 ? MPFR_SIGN (x) :
                                      MPFR_SIGN_POS);
              }
            /* error estimate -- see pow function in algorithms.ps */
            err = Nt - 3;
            if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd)))
              break;

            /* actualisation of the precision */
            Nt += BITS_PER_MP_LIMB;
            mpfr_set_prec (t, Nt);
          }
        MPFR_ZIV_FREE (loop);

        inexact = mpfr_set (y, t, rnd);
        mpfr_clear (t);
        MPFR_SAVE_EXPO_FREE (expo);
        return mpfr_check_range (y, inexact, rnd);
      }
    }
}
开发者ID:mmanley,项目名称:Antares,代码行数:101,代码来源:pow_si.c

示例12: mpfr_pow_general


//.........这里部分代码省略.........
              mpfr_log (t, absx, MPFR_IS_NEG (y) ? MPFR_RNDU : MPFR_RNDD);
              mpfr_mul (t, y, t, MPFR_RNDD);            /* y * ln|x| */
              MPFR_BLOCK (flags2, mpfr_exp (t, t, MPFR_RNDD));
              /* t = lower bound on exp(y * ln|x|) */
              if (MPFR_OVERFLOW (flags2))
                {
                  /* We have computed a lower bound on |x|^y, and it
                     overflowed. Therefore we have a real overflow
                     on |x|^y. */
                  inexact = mpfr_overflow (z, rnd_mode, MPFR_SIGN_POS);
                  if (expo != NULL)
                    MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, MPFR_FLAGS_INEXACT
                                                 | MPFR_FLAGS_OVERFLOW);
                  break;
                }
            }

          k_non_zero = 1;
          Ntmin = sizeof(mpfr_exp_t) * CHAR_BIT;
          if (Ntmin > Nt)
            {
              Nt = Ntmin;
              mpfr_set_prec (t, Nt);
            }
          mpfr_init2 (u, Nt);
          mpfr_init2 (k, Ntmin);
          mpfr_log2 (k, absx, MPFR_RNDN);
          mpfr_mul (k, y, k, MPFR_RNDN);
          mpfr_round (k, k);
          MPFR_LOG_VAR (k);
          /* |y| < 2^Ntmin, therefore |k| < 2^Nt. */
          continue;
        }
      if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - err, Nz, rnd_mode)))
        {
          inexact = mpfr_set (z, t, rnd_mode);
          break;
        }

      /* check exact power, except when y is an integer (since the
         exact cases for y integer have already been filtered out) */
      if (check_exact_case == 0 && ! y_is_integer)
        {
          if (mpfr_pow_is_exact (z, absx, y, rnd_mode, &inexact))
            break;
          check_exact_case = 1;
        }

      /* reactualisation of the precision */
      MPFR_ZIV_NEXT (ziv_loop, Nt);
      mpfr_set_prec (t, Nt);
      if (k_non_zero)
        mpfr_set_prec (u, Nt);
    }
  MPFR_ZIV_FREE (ziv_loop);

  if (k_non_zero)
    {
      int inex2;
      long lk;

      /* The rounded result in an unbounded exponent range is z * 2^k. As
       * MPFR chooses underflow after rounding, the mpfr_mul_2si below will
       * correctly detect underflows and overflows. However, in rounding to
       * nearest, if z * 2^k = 2^(emin - 2), then the double rounding may
       * affect the result. We need to cope with that before overwriting z.
开发者ID:Distrotech,项目名称:mpfr,代码行数:67,代码来源:pow.c

示例13: mpfr_sin


//.........这里部分代码省略.........

  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);
      mpfr_sqrt (c, c, MPFR_RNDZ);
      if (MPFR_IS_NEG_SIGN(sign))
        MPFR_CHANGE_SIGN(c);

      /* Warning: c may be 0! */
      if (MPFR_UNLIKELY (MPFR_IS_ZERO (c)))
        {
          /* Huge cancellation: increase prec a lot! */
          m = MAX (m, MPFR_PREC (x));
          m = 2 * m;
        }
      else
        {
          /* the absolute error on c is at most 2^(3-m-EXP(c)),
             plus 2^(2-m) if there was an argument reduction.
             Since EXP(c) <= 1, 3-m-EXP(c) >= 2-m, thus the error
             is at most 2^(3-m-EXP(c)) in case of argument reduction. */
          err = 2 * MPFR_GET_EXP (c) + (mpfr_exp_t) m - 3 - (reduce != 0);
          if (MPFR_CAN_ROUND (c, err, precy, rnd_mode))
            break;

          /* check for huge cancellation (Near 0) */
          if (err < (mpfr_exp_t) MPFR_PREC (y))
            m += MPFR_PREC (y) - err;
          /* Check if near 1 */
          if (MPFR_GET_EXP (c) == 1)
            m += m;
        }

    ziv_next:
      /* Else generic increase */
      MPFR_ZIV_NEXT (loop, m);
    }
  MPFR_ZIV_FREE (loop);

  inexact = mpfr_set (y, c, rnd_mode);
  /* inexact cannot be 0, since this would mean that c was representable
     within the target precision, but in that case mpfr_can_round will fail */

  mpfr_clear (c);
  mpfr_clear (xr);

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (y, inexact, rnd_mode);
}
开发者ID:119,项目名称:aircam-openwrt,代码行数:101,代码来源:sin.c

示例14: mpfr_log


//.........这里部分代码省略.........
              MPFR_RET_NAN;
            }
          else /* log(+Inf) = +Inf */
            {
              MPFR_SET_INF (r);
              MPFR_SET_POS (r);
              MPFR_RET (0);
            }
        }
      else /* a is zero */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (a));
          MPFR_SET_INF (r);
          MPFR_SET_NEG (r);
          mpfr_set_divby0 ();
          MPFR_RET (0); /* log(0) is an exact -infinity */
        }
    }
  /* If a is negative, the result is NaN */
  else if (MPFR_UNLIKELY (MPFR_IS_NEG (a)))
    {
      MPFR_SET_NAN (r);
      MPFR_RET_NAN;
    }
  /* If a is 1, the result is 0 */
  else if (MPFR_UNLIKELY (MPFR_GET_EXP (a) == 1 && mpfr_cmp_ui (a, 1) == 0))
    {
      MPFR_SET_ZERO (r);
      MPFR_SET_POS (r);
      MPFR_RET (0); /* only "normal" case where the result is exact */
    }

  q = MPFR_PREC (r);

  /* use initial precision about q+lg(q)+5 */
  p = q + 5 + 2 * MPFR_INT_CEIL_LOG2 (q);
  /* % ~(mpfr_prec_t)GMP_NUMB_BITS  ;
     m=q; while (m) { p++; m >>= 1; }  */
  /* if (MPFR_LIKELY(p % GMP_NUMB_BITS != 0))
      p += GMP_NUMB_BITS - (p%GMP_NUMB_BITS); */

  MPFR_SAVE_EXPO_MARK (expo);
  MPFR_GROUP_INIT_2 (group, p, tmp1, tmp2);

  MPFR_ZIV_INIT (loop, p);
  for (;;)
    {
      long m;
      mpfr_exp_t cancel;

      /* Calculus of m (depends on p) */
      m = (p + 1) / 2 - MPFR_GET_EXP (a) + 1;

      mpfr_mul_2si (tmp2, a, m, MPFR_RNDN);    /* s=a*2^m,        err<=1 ulp  */
      mpfr_div (tmp1, __gmpfr_four, tmp2, MPFR_RNDN);/* 4/s,      err<=2 ulps */
      mpfr_agm (tmp2, __gmpfr_one, tmp1, MPFR_RNDN); /* AG(1,4/s),err<=3 ulps */
      mpfr_mul_2ui (tmp2, tmp2, 1, MPFR_RNDN); /* 2*AG(1,4/s),    err<=3 ulps */
      mpfr_const_pi (tmp1, MPFR_RNDN);         /* compute pi,     err<=1ulp   */
      mpfr_div (tmp2, tmp1, tmp2, MPFR_RNDN);  /* pi/2*AG(1,4/s), err<=5ulps  */
      mpfr_const_log2 (tmp1, MPFR_RNDN);      /* compute log(2),  err<=1ulp   */
      mpfr_mul_si (tmp1, tmp1, m, MPFR_RNDN); /* compute m*log(2),err<=2ulps  */
      mpfr_sub (tmp1, tmp2, tmp1, MPFR_RNDN); /* log(a),    err<=7ulps+cancel */

      if (MPFR_LIKELY (MPFR_IS_PURE_FP (tmp1) && MPFR_IS_PURE_FP (tmp2)))
        {
          cancel = MPFR_GET_EXP (tmp2) - MPFR_GET_EXP (tmp1);
          MPFR_LOG_MSG (("canceled bits=%ld\n", (long) cancel));
          MPFR_LOG_VAR (tmp1);
          if (MPFR_UNLIKELY (cancel < 0))
            cancel = 0;

          /* we have 7 ulps of error from the above roundings,
             4 ulps from the 4/s^2 second order term,
             plus the canceled bits */
          if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp1, p-cancel-4, q, rnd_mode)))
            break;

          /* VL: I think it is better to have an increment that it isn't
             too low; in particular, the increment must be positive even
             if cancel = 0 (can this occur?). */
          p += cancel >= 8 ? cancel : 8;
        }
      else
        {
          /* TODO: find why this case can occur and what is best to do
             with it. */
          p += 32;
        }

      MPFR_ZIV_NEXT (loop, p);
      MPFR_GROUP_REPREC_2 (group, p, tmp1, tmp2);
    }
  MPFR_ZIV_FREE (loop);
  inexact = mpfr_set (r, tmp1, rnd_mode);
  /* We clean */
  MPFR_GROUP_CLEAR (group);

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (r, inexact, rnd_mode);
}
开发者ID:pgundlach,项目名称:LuaTeX,代码行数:101,代码来源:log.c

示例15: mpfr_asin

int
mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
  mpfr_t xp;
  int compared, inexact;
  mpfr_prec_t prec;
  mpfr_exp_t xp_exp;
  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),
    ("asin[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (asin), mpfr_log_prec, asin,
     inexact));

  /* Special cases */
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    {
      if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
        {
          MPFR_SET_NAN (asin);
          MPFR_RET_NAN;
        }
      else /* x = 0 */
        {
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
          MPFR_SET_ZERO (asin);
          MPFR_SET_SAME_SIGN (asin, x);
          MPFR_RET (0); /* exact result */
        }
    }

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

  /* Set x_p=|x| (x is a normal number) */
  mpfr_init2 (xp, MPFR_PREC (x));
  inexact = mpfr_abs (xp, x, MPFR_RNDN);
  MPFR_ASSERTD (inexact == 0);

  compared = mpfr_cmp_ui (xp, 1);

  MPFR_SAVE_EXPO_MARK (expo);

  if (MPFR_UNLIKELY (compared >= 0))
    {
      mpfr_clear (xp);
      if (compared > 0)                  /* asin(x) = NaN for |x| > 1 */
        {
          MPFR_SAVE_EXPO_FREE (expo);
          MPFR_SET_NAN (asin);
          MPFR_RET_NAN;
        }
      else                              /* x = 1 or x = -1 */
        {
          if (MPFR_IS_POS (x)) /* asin(+1) = Pi/2 */
            inexact = mpfr_const_pi (asin, rnd_mode);
          else /* asin(-1) = -Pi/2 */
            {
              inexact = -mpfr_const_pi (asin, MPFR_INVERT_RND(rnd_mode));
              MPFR_CHANGE_SIGN (asin);
            }
          mpfr_div_2ui (asin, asin, 1, rnd_mode);
        }
    }
  else
    {
      /* Compute exponent of 1 - ABS(x) */
      mpfr_ui_sub (xp, 1, xp, MPFR_RNDD);
      MPFR_ASSERTD (MPFR_GET_EXP (xp) <= 0);
      MPFR_ASSERTD (MPFR_GET_EXP (x) <= 0);
      xp_exp = 2 - MPFR_GET_EXP (xp);

      /* Set up initial prec */
      prec = MPFR_PREC (asin) + 10 + xp_exp;

      /* use asin(x) = atan(x/sqrt(1-x^2)) */
      MPFR_ZIV_INIT (loop, prec);
      for (;;)
        {
          mpfr_set_prec (xp, prec);
          mpfr_sqr (xp, x, MPFR_RNDN);
          mpfr_ui_sub (xp, 1, xp, MPFR_RNDN);
          mpfr_sqrt (xp, xp, MPFR_RNDN);
          mpfr_div (xp, x, xp, MPFR_RNDN);
          mpfr_atan (xp, xp, MPFR_RNDN);
          if (MPFR_LIKELY (MPFR_CAN_ROUND (xp, prec - xp_exp,
                                           MPFR_PREC (asin), rnd_mode)))
            break;
          MPFR_ZIV_NEXT (loop, prec);
        }
      MPFR_ZIV_FREE (loop);
      inexact = mpfr_set (asin, xp, rnd_mode);

      mpfr_clear (xp);
    }

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


注:本文中的MPFR_CAN_ROUND函数示例由纯净天空整理自Github/MSDocs等开源代码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。