本文整理汇总了C++中MPFR_CAN_ROUND函数的典型用法代码示例。如果您正苦于以下问题：C++ MPFR_CAN_ROUND函数的具体用法？C++ MPFR_CAN_ROUND怎么用？C++ MPFR_CAN_ROUND使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了MPFR_CAN_ROUND函数的20个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于我们的系统推荐出更棒的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
示例16: mpfr_log2
int
mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mpfr_rnd_t rnd_mode)
{
int inexact;
MPFR_SAVE_EXPO_DECL (expo);
if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
{
/* If a is NaN, the result is NaN */
if (MPFR_IS_NAN (a))
{
MPFR_SET_NAN (r);
MPFR_RET_NAN;
}
/* check for infinity before zero */
else if (MPFR_IS_INF (a))
{
if (MPFR_IS_NEG (a))
/* log(-Inf) = NaN */
{
MPFR_SET_NAN (r);
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_RET (0); /* log2(0) is an exact -infinity */
}
}
/* If a is negative, the result is NaN */
if (MPFR_UNLIKELY (MPFR_IS_NEG (a)))
{
MPFR_SET_NAN (r);
MPFR_RET_NAN;
}
/* If a is 1, the result is 0 */
if (MPFR_UNLIKELY (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 */
}
/* If a is 2^N, log2(a) is exact*/
if (MPFR_UNLIKELY (mpfr_cmp_ui_2exp (a, 1, MPFR_GET_EXP (a) - 1) == 0))
return mpfr_set_si(r, MPFR_GET_EXP (a) - 1, rnd_mode);
MPFR_SAVE_EXPO_MARK (expo);
/* General case */
{
/* Declaration of the intermediary variable */
mpfr_t t, tt;
/* Declaration of the size variable */
mpfr_prec_t Ny = MPFR_PREC(r); /* 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 + 3 + MPFR_INT_CEIL_LOG2 (Ny);
/* initialise of intermediary variable */
mpfr_init2 (t, Nt);
mpfr_init2 (tt, Nt);
/* First computation of log2 */
MPFR_ZIV_INIT (loop, Nt);
for (;;)
{
/* compute log2 */
mpfr_const_log2(t,MPFR_RNDD); /* log(2) */
mpfr_log(tt,a,MPFR_RNDN);
|
客服电话
APP下载
官方微信
请发表评论