本文整理汇总了C++中MPFR_ZIV_INIT函数的典型用法代码示例。如果您正苦于以下问题：C++ MPFR_ZIV_INIT函数的具体用法？C++ MPFR_ZIV_INIT怎么用？C++ MPFR_ZIV_INIT使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了MPFR_ZIV_INIT函数的20个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: 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
示例2: mpfr_tan
/* computes tan(x) = sign(x)*sqrt(1/cos(x)^2-1) */
int
mpfr_tan (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
{
mp_prec_t precy, m;
int inexact;
mpfr_t s, c;
MPFR_ZIV_DECL (loop);
MPFR_SAVE_EXPO_DECL (expo);
MPFR_GROUP_DECL (group);
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);
}
}
/* tan(x) = x + x^3/3 + ... so the error is < 2^(3*EXP(x)-1) */
MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2 * MPFR_GET_EXP (x), 1, 1,
rnd_mode, {});
MPFR_SAVE_EXPO_MARK (expo);
/* Compute initial precision */
precy = MPFR_PREC (y);
m = precy + MPFR_INT_CEIL_LOG2 (precy) + 13;
MPFR_ASSERTD (m >= 2); /* needed for the error analysis in algorithms.tex */
MPFR_GROUP_INIT_2 (group, m, s, c);
MPFR_ZIV_INIT (loop, m);
for (;;)
{
/* The only way to get an overflow is to get ~ Pi/2
But the result will be ~ 2^Prec(y). */
mpfr_sin_cos (s, c, x, GMP_RNDN); /* err <= 1/2 ulp on s and c */
mpfr_div (c, s, c, GMP_RNDN); /* err <= 4 ulps */
MPFR_ASSERTD (!MPFR_IS_SINGULAR (c));
if (MPFR_LIKELY (MPFR_CAN_ROUND (c, m - 2, precy, rnd_mode)))
break;
MPFR_ZIV_NEXT (loop, m);
MPFR_GROUP_REPREC_2 (group, m, s, c);
}
MPFR_ZIV_FREE (loop);
inexact = mpfr_set (y, c, rnd_mode);
MPFR_GROUP_CLEAR (group);
MPFR_SAVE_EXPO_FREE (expo);
return mpfr_check_range (y, inexact, rnd_mode);
}
开发者ID:Scorpiion，项目名称:Renux_cross_gcc，代码行数:62，代码来源:tan.c
示例3: 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
示例4: mpfr_exp2
int
mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
int inexact;
long xint;
mpfr_t xfrac;
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,
inexact));
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))
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)) */
//.........这里部分代码省略.........
开发者ID:Canar，项目名称:mpfr，代码行数:101，代码来源:exp2.c
示例5: mpfr_const_log2_internal
/* Don't need to save / restore exponent range: the cache does it */
int
mpfr_const_log2_internal (mpfr_ptr x, mpfr_rnd_t rnd_mode)
{
unsigned long n = MPFR_PREC (x);
mpfr_prec_t w; /* working precision */
unsigned long N;
mpz_t *T, *P, *Q;
mpfr_t t, q;
int inexact;
int ok = 1; /* ensures that the 1st try will give correct rounding */
unsigned long lgN, i;
MPFR_GROUP_DECL(group);
MPFR_TMP_DECL(marker);
MPFR_ZIV_DECL(loop);
MPFR_LOG_FUNC (
("rnd_mode=%d", rnd_mode),
("x[%Pu]=%.*Rg inex=%d", mpfr_get_prec(x), mpfr_log_prec, x, inexact));
if (n < 1253)
w = n + 10; /* ensures correct rounding for the four rounding modes,
together with N = w / 3 + 1 (see below). */
else if (n < 2571)
w = n + 11; /* idem */
else if (n < 3983)
w = n + 12;
else if (n < 4854)
w = n + 13;
else if (n < 26248)
w = n + 14;
else
{
w = n + 15;
ok = 0;
}
MPFR_TMP_MARK(marker);
MPFR_GROUP_INIT_2(group, w, t, q);
MPFR_ZIV_INIT (loop, w);
for (;;)
{
N = w / 3 + 1; /* Warning: do not change that (even increasing N!)
without checking correct rounding in the above
ranges for n. */
/* the following are needed for error analysis (see algorithms.tex) */
MPFR_ASSERTD(w >= 3 && N >= 2);
lgN = MPFR_INT_CEIL_LOG2 (N) + 1;
T = (mpz_t *) MPFR_TMP_ALLOC (3 * lgN * sizeof (mpz_t));
P = T + lgN;
Q = T + 2*lgN;
for (i = 0; i < lgN; i++)
{
mpz_init (T[i]);
mpz_init (P[i]);
mpz_init (Q[i]);
}
S (T, P, Q, 0, N, 0);
mpfr_set_z (t, T[0], MPFR_RNDN);
mpfr_set_z (q, Q[0], MPFR_RNDN);
mpfr_div (t, t, q, MPFR_RNDN);
for (i = 0; i < lgN; i++)
{
mpz_clear (T[i]);
mpz_clear (P[i]);
mpz_clear (Q[i]);
}
if (MPFR_LIKELY (ok != 0
|| mpfr_can_round (t, w - 2, MPFR_RNDN, rnd_mode, n)))
break;
MPFR_ZIV_NEXT (loop, w);
MPFR_GROUP_REPREC_2(group, w, t, q);
}
MPFR_ZIV_FREE (loop);
inexact = mpfr_set (x, t, rnd_mode);
MPFR_GROUP_CLEAR(group);
MPFR_TMP_FREE(marker);
return inexact;
}
开发者ID:Canar，项目名称:mpfr，代码行数:90，代码来源:const_log2.c
示例6: mpfr_sin_cos
/* (y, z) <- (sin(x), cos(x)), return value is 0 iff both results are exact
ie, iff x = 0 */
int
mpfr_sin_cos (mpfr_ptr y, mpfr_ptr z, mpfr_srcptr x, mp_rnd_t rnd_mode)
{
mp_prec_t prec, m;
int neg, reduce;
mpfr_t c, xr;
mpfr_srcptr xx;
mp_exp_t err, expx;
MPFR_ZIV_DECL (loop);
if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
{
if (MPFR_IS_NAN(x) || MPFR_IS_INF(x))
{
MPFR_SET_NAN (y);
MPFR_SET_NAN (z);
MPFR_RET_NAN;
}
else /* x is zero */
{
MPFR_ASSERTD (MPFR_IS_ZERO (x));
MPFR_SET_ZERO (y);
MPFR_SET_SAME_SIGN (y, x);
/* y = 0, thus exact, but z is inexact in case of underflow
or overflow */
return mpfr_set_ui (z, 1, rnd_mode);
}
}
MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd_mode),
("sin[%#R]=%R cos[%#R]=%R", y, y, z, z));
prec = MAX (MPFR_PREC (y), MPFR_PREC (z));
m = prec + MPFR_INT_CEIL_LOG2 (prec) + 13;
expx = MPFR_GET_EXP (x);
mpfr_init (c);
mpfr_init (xr);
MPFR_ZIV_INIT (loop, m);
for (;;)
{
/* the following is copied from sin.c */
if (expx >= 2) /* reduce the argument */
{
reduce = 1;
mpfr_set_prec (c, expx + m - 1);
mpfr_set_prec (xr, m);
mpfr_const_pi (c, GMP_RNDN);
mpfr_mul_2ui (c, c, 1, GMP_RNDN);
mpfr_remainder (xr, x, c, GMP_RNDN);
mpfr_div_2ui (c, c, 1, GMP_RNDN);
if (MPFR_SIGN (xr) > 0)
mpfr_sub (c, c, xr, GMP_RNDZ);
else
mpfr_add (c, c, xr, GMP_RNDZ);
if (MPFR_IS_ZERO(xr) || MPFR_EXP(xr) < (mp_exp_t) 3 - (mp_exp_t) m
|| MPFR_EXP(c) < (mp_exp_t) 3 - (mp_exp_t) m)
goto next_step;
xx = xr;
}
else /* the input argument is already reduced */
{
reduce = 0;
xx = x;
}
neg = MPFR_IS_NEG (xx); /* gives sign of sin(x) */
mpfr_set_prec (c, m);
mpfr_cos (c, xx, GMP_RNDZ);
/* If no argument reduction was performed, the error is at most ulp(c),
otherwise it is at most ulp(c) + 2^(2-m). Since |c| < 1, we have
ulp(c) <= 2^(-m), thus the error is bounded by 2^(3-m) in that later
case. */
if (reduce == 0)
err = m;
else
err = MPFR_GET_EXP (c) + (mp_exp_t) (m - 3);
if (!mpfr_can_round (c, err, GMP_RNDN, rnd_mode,
MPFR_PREC (z) + (rnd_mode == GMP_RNDN)))
goto next_step;
mpfr_set (z, c, rnd_mode);
mpfr_sqr (c, c, GMP_RNDU);
mpfr_ui_sub (c, 1, c, GMP_RNDN);
err = 2 + (- MPFR_GET_EXP (c)) / 2;
mpfr_sqrt (c, c, GMP_RNDN);
if (neg)
MPFR_CHANGE_SIGN (c);
/* the absolute error on c is at most 2^(err-m), which we must put
in the form 2^(EXP(c)-err). If there was an argument reduction,
we need to add 2^(2-m); since err >= 2, the error is bounded by
2^(err+1-m) in that case. */
err = MPFR_GET_EXP (c) + (mp_exp_t) m - (err + reduce);
if (mpfr_can_round (c, err, GMP_RNDN, rnd_mode,
MPFR_PREC (y) + (rnd_mode == GMP_RNDN)))
break;
//.........这里部分代码省略.........
开发者ID:mmanley，项目名称:Antares，代码行数:101，代码来源:sin_cos.c
示例7: mpfr_log1p
int
mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
{
int comp, inexact;
mp_exp_t ex;
MPFR_SAVE_EXPO_DECL (expo);
if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
{
if (MPFR_IS_NAN (x))
{
MPFR_SET_NAN (y);
MPFR_RET_NAN;
}
/* check for inf or -inf (result is not defined) */
else if (MPFR_IS_INF (x))
{
if (MPFR_IS_POS (x))
{
MPFR_SET_INF (y);
MPFR_SET_POS (y);
MPFR_RET (0);
}
else
{
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;
}
//.........这里部分代码省略.........
开发者ID:Scorpiion，项目名称:Renux_cross_gcc，代码行数:101，代码来源:log1p.c
示例8: 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
示例9: mpfr_sinh_cosh
int
mpfr_sinh_cosh (mpfr_ptr sh, mpfr_ptr ch, mpfr_srcptr xt, mpfr_rnd_t rnd_mode)
{
mpfr_t x;
int inexact_sh, inexact_ch;
MPFR_ASSERTN (sh != ch);
MPFR_LOG_FUNC
(("x[%Pu]=%.*Rg rnd=%d",
mpfr_get_prec (xt), mpfr_log_prec, xt, rnd_mode),
("sh[%Pu]=%.*Rg ch[%Pu]=%.*Rg",
mpfr_get_prec (sh), mpfr_log_prec, sh,
mpfr_get_prec (ch), mpfr_log_prec, ch));
if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
{
if (MPFR_IS_NAN (xt))
{
MPFR_SET_NAN (ch);
MPFR_SET_NAN (sh);
MPFR_RET_NAN;
}
else if (MPFR_IS_INF (xt))
{
MPFR_SET_INF (sh);
MPFR_SET_SAME_SIGN (sh, xt);
MPFR_SET_INF (ch);
MPFR_SET_POS (ch);
MPFR_RET (0);
}
else /* xt is zero */
{
MPFR_ASSERTD (MPFR_IS_ZERO (xt));
MPFR_SET_ZERO (sh); /* sinh(0) = 0 */
MPFR_SET_SAME_SIGN (sh, xt);
inexact_sh = 0;
inexact_ch = mpfr_set_ui (ch, 1, rnd_mode); /* cosh(0) = 1 */
return INEX(inexact_sh,inexact_ch);
}
}
/* Warning: if we use MPFR_FAST_COMPUTE_IF_SMALL_INPUT here, make sure
that the code also works in case of overlap (see sin_cos.c) */
MPFR_TMP_INIT_ABS (x, xt);
{
mpfr_t s, c, ti;
mpfr_exp_t d;
mpfr_prec_t N; /* Precision of the intermediary variables */
long int err; /* Precision of error */
MPFR_ZIV_DECL (loop);
MPFR_SAVE_EXPO_DECL (expo);
MPFR_GROUP_DECL (group);
MPFR_SAVE_EXPO_MARK (expo);
/* compute the precision of intermediary variable */
N = MPFR_PREC (ch);
N = MAX (N, MPFR_PREC (sh));
/* the optimal number of bits : see algorithms.ps */
N = N + MPFR_INT_CEIL_LOG2 (N) + 4;
/* initialise of intermediary variables */
MPFR_GROUP_INIT_3 (group, N, s, c, ti);
/* First computation of sinh_cosh */
MPFR_ZIV_INIT (loop, N);
for (;;)
{
MPFR_BLOCK_DECL (flags);
/* compute sinh_cosh */
MPFR_BLOCK (flags, mpfr_exp (s, x, MPFR_RNDD));
if (MPFR_OVERFLOW (flags))
/* exp(x) does overflow */
{
/* since cosh(x) >= exp(x), cosh(x) overflows too */
inexact_ch = mpfr_overflow (ch, rnd_mode, MPFR_SIGN_POS);
/* sinh(x) may be representable */
inexact_sh = mpfr_sinh (sh, xt, rnd_mode);
MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
break;
}
d = MPFR_GET_EXP (s);
mpfr_ui_div (ti, 1, s, MPFR_RNDU); /* 1/exp(x) */
mpfr_add (c, s, ti, MPFR_RNDU); /* exp(x) + 1/exp(x) */
mpfr_sub (s, s, ti, MPFR_RNDN); /* exp(x) - 1/exp(x) */
mpfr_div_2ui (c, c, 1, MPFR_RNDN); /* 1/2(exp(x) + 1/exp(x)) */
mpfr_div_2ui (s, s, 1, MPFR_RNDN); /* 1/2(exp(x) - 1/exp(x)) */
/* it may be that s is zero (in fact, it can only occur when exp(x)=1,
and thus ti=1 too) */
if (MPFR_IS_ZERO (s))
err = N; /* double the precision */
else
{
/* calculation of the error */
d = d - MPFR_GET_EXP (s) + 2;
//.........这里部分代码省略.........
开发者ID:aimanqais，项目名称:gerardus，代码行数:101，代码来源:sinh_cosh.c
示例10: mpfr_cos
int
mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
mpfr_prec_t K0, K, precy, m, k, l;
int inexact, reduce = 0;
mpfr_t r, s, xr, c;
mpfr_exp_t exps, cancel = 0, expx;
MPFR_ZIV_DECL (loop);
MPFR_SAVE_EXPO_DECL (expo);
MPFR_GROUP_DECL (group);
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,
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
{
MPFR_ASSERTD (MPFR_IS_ZERO (x));
return mpfr_set_ui (y, 1, rnd_mode);
}
}
MPFR_SAVE_EXPO_MARK (expo);
/* cos(x) = 1-x^2/2 + ..., so error < 2^(2*EXP(x)-1) */
expx = MPFR_GET_EXP (x);
MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, __gmpfr_one, -2 * expx,
1, 0, rnd_mode, expo, {});
/* Compute initial precision */
precy = MPFR_PREC (y);
if (precy >= MPFR_SINCOS_THRESHOLD)
{
MPFR_SAVE_EXPO_FREE (expo);
return mpfr_cos_fast (y, x, rnd_mode);
}
K0 = __gmpfr_isqrt (precy / 3);
m = precy + 2 * MPFR_INT_CEIL_LOG2 (precy) + 2 * K0;
if (expx >= 3)
{
reduce = 1;
/* As expx + m - 1 will silently be converted into mpfr_prec_t
in the mpfr_init2 call, the assert below may be useful to
avoid undefined behavior. */
MPFR_ASSERTN (expx + m - 1 <= MPFR_PREC_MAX);
mpfr_init2 (c, expx + m - 1);
mpfr_init2 (xr, m);
}
MPFR_GROUP_INIT_2 (group, m, r, s);
MPFR_ZIV_INIT (loop, m);
for (;;)
{
/* If |x| >= 4, first reduce x cmod (2*Pi) into xr, using mpfr_remainder:
let e = EXP(x) >= 3, and m the target precision:
(1) c <- 2*Pi [precision e+m-1, nearest]
(2) xr <- remainder (x, c) [precision m, nearest]
We have |c - 2*Pi| <= 1/2ulp(c) = 2^(3-e-m)
|xr - x - k c| <= 1/2ulp(xr) <= 2^(1-m)
|k| <= |x|/(2*Pi) <= 2^(e-2)
Thus |xr - x - 2kPi| <= |k| |c - 2Pi| + 2^(1-m) <= 2^(2-m).
It follows |cos(xr) - cos(x)| <= 2^(2-m). */
if (reduce)
{
mpfr_const_pi (c, MPFR_RNDN);
mpfr_mul_2ui (c, c, 1, MPFR_RNDN); /* 2Pi */
mpfr_remainder (xr, x, c, MPFR_RNDN);
if (MPFR_IS_ZERO(xr))
goto ziv_next;
/* now |xr| <= 4, thus r <= 16 below */
mpfr_mul (r, xr, xr, MPFR_RNDU); /* err <= 1 ulp */
}
else
mpfr_mul (r, x, x, MPFR_RNDU); /* err <= 1 ulp */
/* now |x| < 4 (or xr if reduce = 1), thus |r| <= 16 */
/* we need |r| < 1/2 for mpfr_cos2_aux, i.e., EXP(r) - 2K <= -1 */
K = K0 + 1 + MAX(0, MPFR_GET_EXP(r)) / 2;
/* since K0 >= 0, if EXP(r) < 0, then K >= 1, thus EXP(r) - 2K <= -3;
otherwise if EXP(r) >= 0, then K >= 1/2 + EXP(r)/2, thus
EXP(r) - 2K <= -1 */
MPFR_SET_EXP (r, MPFR_GET_EXP (r) - 2 * K); /* Can't overflow! */
/* s <- 1 - r/2! + ... + (-1)^l r^l/(2l)! */
l = mpfr_cos2_aux (s, r);
/* l is the error bound in ulps on s */
MPFR_SET_ONE (r);
//.........这里部分代码省略.........
开发者ID:Distrotech，项目名称:mpfr，代码行数:101，代码来源:cos.c
示例11: mpfr_pow_si
//.........这里部分代码省略.........
/* General case */
{
/* Declaration of the intermediary variable */
mpfr_t t;
/* Declaration of the size variable */
mpfr_prec_t Ny; /* target precision */
mpfr_prec_t Nt; /* working precision */
mpfr_rnd_t rnd1;
int size_n;
int inexact;
unsigned long abs_n;
MPFR_SAVE_EXPO_DECL (expo);
MPFR_ZIV_DECL (loop);
abs_n = - (unsigned long) n;
count_leading_zeros (size_n, (mp_limb_t) abs_n);
size_n = GMP_NUMB_BITS - size_n;
/* initial working precision */
Ny = MPFR_PREC (y);
Nt = Ny + size_n + 3 + MPFR_INT_CEIL_LOG2 (Ny);
MPFR_SAVE_EXPO_MARK (expo);
/* initialise of intermediary variable */
mpfr_init2 (t, Nt);
/* We will compute rnd(rnd1(1/x) ^ |n|), where rnd1 is the rounding
toward sign(x), to avoid spurious overflow or underflow, as in
mpfr_pow_z. */
rnd1 = MPFR_EXP (x) < 1 ? MPFR_RNDZ :
(MPFR_SIGN (x) > 0 ? MPFR_RNDU : MPFR_RNDD);
MPFR_ZIV_INIT (loop, Nt);
for (;;)
{
MPFR_BLOCK_DECL (flags);
/* compute (1/x)^|n| */
MPFR_BLOCK (flags, mpfr_ui_div (t, 1, x, rnd1));
MPFR_ASSERTD (! MPFR_UNDERFLOW (flags));
/* t = (1/x)*(1+theta) where |theta| <= 2^(-Nt) */
if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
goto overflow;
MPFR_BLOCK (flags, mpfr_pow_ui (t, t, abs_n, rnd));
/* t = (1/x)^|n|*(1+theta')^(|n|+1) where |theta'| <= 2^(-Nt).
If (|n|+1)*2^(-Nt) <= 1/2, which is satisfied as soon as
Nt >= bits(n)+2, then we can use Lemma \ref{lemma_graillat}
from algorithms.tex, which yields x^n*(1+theta) with
|theta| <= 2(|n|+1)*2^(-Nt), thus the error is bounded by
2(|n|+1) ulps <= 2^(bits(n)+2) ulps. */
if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
{
overflow:
MPFR_ZIV_FREE (loop);
mpfr_clear (t);
MPFR_SAVE_EXPO_FREE (expo);
MPFR_LOG_MSG (("overflow\n", 0));
return mpfr_overflow (y, rnd, abs_n & 1 ?
MPFR_SIGN (x) : MPFR_SIGN_POS);
}
if (MPFR_UNLIKELY (MPFR_UNDERFLOW (flags)))
{
MPFR_ZIV_FREE (loop);
mpfr_clear (t);
MPFR_LOG_MSG (("underflow\n", 0));
开发者ID:Kirija，项目名称:XPIR，代码行数:67，代码来源:pow_si.c
示例12: 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
示例13: mpfr_log
int
mpfr_log (mpfr_ptr r, mpfr_srcptr a, mpfr_rnd_t rnd_mode)
{
int inexact;
mpfr_prec_t p, q;
mpfr_t tmp1, tmp2;
MPFR_SAVE_EXPO_DECL (expo);
MPFR_ZIV_DECL (loop);
MPFR_GROUP_DECL(group);
MPFR_LOG_FUNC
(("a[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (a), mpfr_log_prec, a, rnd_mode),
("r[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (r), mpfr_log_prec, r,
inexact));
/* Special cases */
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_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))
//.........这里部分代码省略.........
开发者ID:pgundlach，项目名称:LuaTeX，代码行数:101，代码来源:log.c
示例14: mpfr_ui_pow_ui
int
mpfr_ui_pow_ui (mpfr_ptr x, unsigned long int y, unsigned long int n,
mpfr_rnd_t rnd)
{
mpfr_exp_t err;
unsigned long m;
mpfr_t res;
mpfr_prec_t prec;
int size_n;
int inexact;
MPFR_ZIV_DECL (loop);
MPFR_SAVE_EXPO_DECL (expo);
if (MPFR_UNLIKELY (n <= 1))
{
if (n == 1)
return mpfr_set_ui (x, y, rnd); /* y^1 = y */
else
return mpfr_set_ui (x, 1, rnd); /* y^0 = 1 for any y */
}
else if (MPFR_UNLIKELY (y <= 1))
{
if (y == 1)
return mpfr_set_ui (x, 1, rnd); /* 1^n = 1 for any n > 0 */
else
return mpfr_set_ui (x, 0, rnd); /* 0^n = 0 for any n > 0 */
}
for (size_n = 0, m = n; m; size_n++, m >>= 1);
MPFR_SAVE_EXPO_MARK (expo);
prec = MPFR_PREC (x) + 3 + size_n;
mpfr_init2 (res, prec);
MPFR_ZIV_INIT (loop, prec);
for (;;)
{
int i = size_n;
inexact = mpfr_set_ui (res, y, MPFR_RNDU);
err = 1;
/* now 2^(i-1) <= n < 2^i: i=1+floor(log2(n)) */
for (i -= 2; i >= 0; i--)
{
inexact |= mpfr_mul (res, res, res, MPFR_RNDU);
err++;
if (n & (1UL << i))
inexact |= mpfr_mul_ui (res, res, y, MPFR_RNDU);
}
/* since the loop is executed floor(log2(n)) times,
we have err = 1+floor(log2(n)).
Since prec >= MPFR_PREC(x) + 4 + floor(log2(n)), prec > err */
err = prec - err;
if (MPFR_LIKELY (inexact == 0
|| MPFR_CAN_ROUND (res, err, MPFR_PREC (x), rnd)))
break;
/* Actualisation of the precision */
MPFR_ZIV_NEXT (loop, prec);
mpfr_set_prec (res, prec);
}
MPFR_ZIV_FREE (loop);
inexact = mpfr_set (x, res, rnd);
mpfr_clear (res);
MPFR_SAVE_EXPO_FREE (expo);
return mpfr_check_range (x, inexact, rnd);
}
开发者ID:Kirija，项目名称:XPIR，代码行数:71，代码来源:ui_pow_ui.c
示例15: mpfr_log10
int
mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mpfr_rnd_t rnd_mode)
{
int inexact;
MPFR_SAVE_EXPO_DECL (expo);
/* If a is NaN, the result is NaN */
if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
{
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))
/* log10(-Inf) = NaN */
{
MPFR_SET_NAN (r);
MPFR_RET_NAN;
}
else /* log10(+Inf) = +Inf */
{
MPFR_SET_INF (r);
MPFR_SET_POS (r);
MPFR_RET (0); /* exact */
}
}
else /* a = 0 */
{
MPFR_ASSERTD (MPFR_IS_ZERO (a));
MPFR_SET_INF (r);
MPFR_SET_NEG (r);
MPFR_RET (0); /* log10(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_cmp_ui (a, 1) == 0)
{
MPFR_SET_ZERO (r);
MPFR_SET_POS (r);
MPFR_RET (0); /* result is exact */
}
MPFR_SAVE_EXPO_MARK (expo);
/* General case */
{
/* Declaration of the intermediary variable */
mpfr_t t, tt;
MPFR_ZIV_DECL (loop);
/* Declaration of the size variable */
mpfr_prec_t Ny = MPFR_PREC(r); /* Precision of output variable */
mpfr_prec_t Nt; /* Precision of the intermediary variable */
mpfr_exp_t err; /* Precision of error */
/* compute the precision of intermediary variable */
/* the optimal number of bits : see algorithms.tex */
Nt = Ny + 4 + MPFR_INT_CEIL_LOG2 (Ny);
/* initialise of intermediary variables */
mpfr_init2 (t, Nt);
mpfr_init2 (tt, Nt);
/* First computation of log10 */
MPFR_ZIV_INIT (loop, Nt);
for (;;)
{
/* compute log10 */
mpfr_set_ui (t, 10, MPFR_RNDN); /* 10 */
mpfr_log (t, t, MPFR_RNDD); /* log(10) */
mpfr_log (tt, a, MPFR_RNDN); /* log(a) */
mpfr_div (t, tt, t, MPFR_RNDN); /* log(a)/log(10) */
/* estimation of the error */
err = Nt - 4;
if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
break;
/* log10(10^n) is exact:
FIXME: Can we have 10^n exactly representable as a mpfr_t
but n can't fit an unsigned long? */
if (MPFR_IS_POS (t)
&& mpfr_integer_p (t) && mpfr_fits_ulong_p (t, MPFR_RNDN)
&& !mpfr_ui_pow
|
客服电话
APP下载
官方微信
请发表评论