/* Implements asymptotic expansion for jn or yn (formulae 9.2.5 and 9.2.6
from Abramowitz & Stegun).
Assumes |z| > p log(2)/2, where p is the target precision
(z can be negative only for jn).
Return 0 if the expansion does not converge enough (the value 0 as inexact
flag should not happen for normal input).
*/
static int
FUNCTION (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r)
{
mpfr_t s, c, P, Q, t, iz, err_t, err_s, err_u;
mpfr_prec_t w;
long k;
int inex, stop, diverge = 0;
mpfr_exp_t err2, err;
MPFR_ZIV_DECL (loop);
mpfr_init (c);
w = MPFR_PREC(res) + MPFR_INT_CEIL_LOG2(MPFR_PREC(res)) + 4;
MPFR_ZIV_INIT (loop, w);
for (;;)
{
mpfr_set_prec (c, w);
mpfr_init2 (s, w);
mpfr_init2 (P, w);
mpfr_init2 (Q, w);
mpfr_init2 (t, w);
mpfr_init2 (iz, w);
mpfr_init2 (err_t, 31);
mpfr_init2 (err_s, 31);
mpfr_init2 (err_u, 31);
/* Approximate sin(z) and cos(z). In the following, err <= k means that
the approximate value y and the true value x are related by
y = x * (1 + u)^k with |u| <= 2^(-w), following Higham's method. */
mpfr_sin_cos (s, c, z, MPFR_RNDN);
if (MPFR_IS_NEG(z))
mpfr_neg (s, s, MPFR_RNDN); /* compute jn/yn(|z|), fix sign later */
/* The absolute error on s/c is bounded by 1/2 ulp(1/2) <= 2^(-w-1). */
mpfr_add (t, s, c, MPFR_RNDN);
mpfr_sub (c, s, c, MPFR_RNDN);
mpfr_swap (s, t);
/* now s approximates sin(z)+cos(z), and c approximates sin(z)-cos(z),
with total absolute error bounded by 2^(1-w). */
/* precompute 1/(8|z|) */
mpfr_si_div (iz, MPFR_IS_POS(z) ? 1 : -1, z, MPFR_RNDN); /* err <= 1 */
mpfr_div_2ui (iz, iz, 3, MPFR_RNDN);
/* compute P and Q */
mpfr_set_ui (P, 1, MPFR_RNDN);
mpfr_set_ui (Q, 0, MPFR_RNDN);
mpfr_set_ui (t, 1, MPFR_RNDN); /* current term */
mpfr_set_ui (err_t, 0, MPFR_RNDN); /* error on t */
mpfr_set_ui (err_s, 0, MPFR_RNDN); /* error on P and Q (sum of errors) */
for (k = 1, stop = 0; stop < 4; k++)
{
/* compute next term: t(k)/t(k-1) = (2n+2k-1)(2n-2k+1)/(8kz) */
mpfr_mul_si (t, t, 2 * (n + k) - 1, MPFR_RNDN); /* err <= err_k + 1 */
mpfr_mul_si (t, t, 2 * (n - k) + 1, MPFR_RNDN); /* err <= err_k + 2 */
mpfr_div_ui (t, t, k, MPFR_RNDN); /* err <= err_k + 3 */
mpfr_mul (t, t, iz, MPFR_RNDN); /* err <= err_k + 5 */
/* the relative error on t is bounded by (1+u)^(5k)-1, which is
bounded by 6ku for 6ku <= 0.02: first |5 log(1+u)| <= |5.5u|
for |u| <= 0.15, then |exp(5.5u)-1| <= 6u for |u| <= 0.02. */
mpfr_mul_ui (err_t, t, 6 * k, MPFR_IS_POS(t) ? MPFR_RNDU : MPFR_RNDD);
mpfr_abs (err_t, err_t, MPFR_RNDN); /* exact */
/* the absolute error on t is bounded by err_t * 2^(-w) */
mpfr_abs (err_u, t, MPFR_RNDU);
mpfr_mul_2ui (err_u, err_u, w, MPFR_RNDU); /* t * 2^w */
mpfr_add (err_u, err_u, err_t, MPFR_RNDU); /* max|t| * 2^w */
if (stop >= 2)
{
/* take into account the neglected terms: t * 2^w */
mpfr_div_2ui (err_s, err_s, w, MPFR_RNDU);
if (MPFR_IS_POS(t))
mpfr_add (err_s, err_s, t, MPFR_RNDU);
else
mpfr_sub (err_s, err_s, t, MPFR_RNDU);
mpfr_mul_2ui (err_s, err_s, w, MPFR_RNDU);
stop ++;
}
/* if k is odd, add to Q, otherwise to P */
else if (k & 1)
{
/* if k = 1 mod 4, add, otherwise subtract */
if ((k & 2) == 0)
mpfr_add (Q, Q, t, MPFR_RNDN);
else
mpfr_sub (Q, Q, t, MPFR_RNDN);
/* check if the next term is smaller than ulp(Q): if EXP(err_u)
<= EXP(Q), since the current term is bounded by
err_u * 2^(-w), it is bounded by ulp(Q) */
if (MPFR_EXP(err_u) <= MPFR_EXP(Q))
stop ++;
else
stop = 0;
}
//.........这里部分代码省略.........
开发者ID:Kirija，项目名称:XPIR，代码行数:101，代码来源:jyn_asympt.c
示例3: 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;
//.........这里部分代码省略.........
/* Assumes that the exponent range has already been extended and if y is
an integer, then the result is not exact in unbounded exponent range. */
int
mpfr_pow_general (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
mpfr_rnd_t rnd_mode, int y_is_integer, mpfr_save_expo_t *expo)
{
mpfr_t t, u, k, absx;
int neg_result = 0;
int k_non_zero = 0;
int check_exact_case = 0;
int inexact;
/* Declaration of the size variable */
mpfr_prec_t Nz = MPFR_PREC(z); /* target precision */
mpfr_prec_t Nt; /* working precision */
mpfr_exp_t err; /* error */
MPFR_ZIV_DECL (ziv_loop);
MPFR_LOG_FUNC
(("x[%Pu]=%.*Rg y[%Pu]=%.*Rg rnd=%d",
mpfr_get_prec (x), mpfr_log_prec, x,
mpfr_get_prec (y), mpfr_log_prec, y, rnd_mode),
("z[%Pu]=%.*Rg inexact=%d",
mpfr_get_prec (z), mpfr_log_prec, z, inexact));
/* We put the absolute value of x in absx, pointing to the significand
of x to avoid allocating memory for the significand of absx. */
MPFR_ALIAS(absx, x, /*sign=*/ 1, /*EXP=*/ MPFR_EXP(x));
/* We will compute the absolute value of the result. So, let's
invert the rounding mode if the result is negative. */
if (MPFR_IS_NEG (x) && is_odd (y))
{
neg_result = 1;
rnd_mode = MPFR_INVERT_RND (rnd_mode);
}
/* compute the precision of intermediary variable */
/* the optimal number of bits : see algorithms.tex */
Nt = Nz + 5 + MPFR_INT_CEIL_LOG2 (Nz);
/* initialise of intermediary variable */
mpfr_init2 (t, Nt);
MPFR_ZIV_INIT (ziv_loop, Nt);
for (;;)
{
MPFR_BLOCK_DECL (flags1);
/* compute exp(y*ln|x|), using MPFR_RNDU to get an upper bound, so
that we can detect underflows. */
mpfr_log (t, absx, MPFR_IS_NEG (y) ? MPFR_RNDD : MPFR_RNDU); /* ln|x| */
mpfr_mul (t, y, t, MPFR_RNDU); /* y*ln|x| */
if (k_non_zero)
{
MPFR_LOG_MSG (("subtract k * ln(2)\n", 0));
mpfr_const_log2 (u, MPFR_RNDD);
mpfr_mul (u, u, k, MPFR_RNDD);
/* Error on u = k * log(2): < k * 2^(-Nt) < 1. */
mpfr_sub (t, t, u, MPFR_RNDU);
MPFR_LOG_MSG (("t = y * ln|x| - k * ln(2)\n", 0));
MPFR_LOG_VAR (t);
}
/* estimate of the error -- see pow function in algorithms.tex.
The error on t is at most 1/2 + 3*2^(EXP(t)+1) ulps, which is
<= 2^(EXP(t)+3) for EXP(t) >= -1, and <= 2 ulps for EXP(t) <= -2.
Additional error if k_no_zero: treal = t * errk, with
1 - |k| * 2^(-Nt) <= exp(-|k| * 2^(-Nt)) <= errk <= 1,
i.e., additional absolute error <= 2^(EXP(k)+EXP(t)-Nt).
Total error <= 2^err1 + 2^err2 <= 2^(max(err1,err2)+1). */
err = MPFR_NOTZERO (t) && MPFR_GET_EXP (t) >= -1 ?
MPFR_GET_EXP (t) + 3 : 1;
if (k_non_zero)
{
if (MPFR_GET_EXP (k) > err)
err = MPFR_GET_EXP (k);
err++;
}
MPFR_BLOCK (flags1, mpfr_exp (t, t, MPFR_RNDN)); /* exp(y*ln|x|)*/
/* We need to test */
if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (t) || MPFR_UNDERFLOW (flags1)))
{
mpfr_prec_t Ntmin;
MPFR_BLOCK_DECL (flags2);
MPFR_ASSERTN (!k_non_zero);
MPFR_ASSERTN (!MPFR_IS_NAN (t));
/* Real underflow? */
if (MPFR_IS_ZERO (t))
{
/* Underflow. We computed rndn(exp(t)), where t >= y*ln|x|.
Therefore rndn(|x|^y) = 0, and we have a real underflow on
|x|^y. */
inexact = mpfr_underflow (z, rnd_mode == MPFR_RNDN ? MPFR_RNDZ
: rnd_mode, MPFR_SIGN_POS);
if (expo != NULL)
MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, MPFR_FLAGS_INEXACT
| MPFR_FLAGS_UNDERFLOW);
break;
//.........这里部分代码省略.........
开发者ID:Distrotech，项目名称:mpfr，代码行数:101，代码来源:pow.c
示例7: exp
/* The computation of z = pow(x,y) is done by
z = exp(y * log(x)) = x^y
For the special cases, see Section F.9.4.4 of the C standard:
_ pow(±0, y) = ±inf for y an odd integer < 0.
_ pow(±0, y) = +inf for y < 0 and not an odd integer.
_ pow(±0, y) = ±0 for y an odd integer > 0.
_ pow(±0, y) = +0 for y > 0 and not an odd integer.
_ pow(-1, ±inf) = 1.
_ pow(+1, y) = 1 for any y, even a NaN.
_ pow(x, ±0) = 1 for any x, even a NaN.
_ pow(x, y) = NaN for finite x < 0 and finite non-integer y.
_ pow(x, -inf) = +inf for |x| < 1.
_ pow(x, -inf) = +0 for |x| > 1.
_ pow(x, +inf) = +0 for |x| < 1.
_ pow(x, +inf) = +inf for |x| > 1.
_ pow(-inf, y) = -0 for y an odd integer < 0.
_ pow(-inf, y) = +0 for y < 0 and not an odd integer.
_ pow(-inf, y) = -inf for y an odd integer > 0.
_ pow(-inf, y) = +inf for y > 0 and not an odd integer.
_ pow(+inf, y) = +0 for y < 0.
_ pow(+inf, y) = +inf for y > 0. */
int
mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
{
int inexact;
int cmp_x_1;
int y_is_integer;
MPFR_SAVE_EXPO_DECL (expo);
MPFR_LOG_FUNC
(("x[%Pu]=%.*Rg y[%Pu]=%.*Rg rnd=%d",
mpfr_get_prec (x), mpfr_log_prec, x,
mpfr_get_prec (y), mpfr_log_prec, y, rnd_mode),
("z[%Pu]=%.*Rg inexact=%d",
mpfr_get_prec (z), mpfr_log_prec, z, inexact));
if (MPFR_ARE_SINGULAR (x, y))
{
/* pow(x, 0) returns 1 for any x, even a NaN. */
if (MPFR_UNLIKELY (MPFR_IS_ZERO (y)))
return mpfr_set_ui (z, 1, rnd_mode);
else if (MPFR_IS_NAN (x))
{
MPFR_SET_NAN (z);
MPFR_RET_NAN;
}
else if (MPFR_IS_NAN (y))
{
/* pow(+1, NaN) returns 1. */
if (mpfr_cmp_ui (x, 1) == 0)
return mpfr_set_ui (z, 1, rnd_mode);
MPFR_SET_NAN (z);
MPFR_RET_NAN;
}
else if (MPFR_IS_INF (y))
{
if (MPFR_IS_INF (x))
{
if (MPFR_IS_POS (y))
MPFR_SET_INF (z);
else
MPFR_SET_ZERO (z);
MPFR_SET_POS (z);
MPFR_RET (0);
}
else
{
int cmp;
cmp = mpfr_cmpabs (x, __gmpfr_one) * MPFR_INT_SIGN (y);
MPFR_SET_POS (z);
if (cmp > 0)
{
/* Return +inf. */
MPFR_SET_INF (z);
MPFR_RET (0);
}
else if (cmp < 0)
{
/* Return +0. */
MPFR_SET_ZERO (z);
MPFR_RET (0);
}
else
{
/* Return 1. */
return mpfr_set_ui (z, 1, rnd_mode);
}
}
}
else if (MPFR_IS_INF (x))
{
int negative;
/* Determine the sign now, in case y and z are the same object */
negative = MPFR_IS_NEG (x) && is_odd (y);
if (MPFR_IS_POS (y))
MPFR_SET_INF (z);
else
MPFR_SET_ZERO (z);
if (negative)
MPFR_SET_NEG (z);
//.........这里部分代码省略.........
开发者ID:Distrotech，项目名称:mpfr，代码行数:101，代码来源:pow.c
示例8: test_2exp
static void
test_2exp (void)
{
mpfr_t x;
int res;
mpfr_init2 (x, 32);
mpfr_set_ui_2exp (x, 1, 0, GMP_RNDN);
if (mpfr_cmp_ui(x, 1))
ERROR("(1U,0)");
mpfr_set_ui_2exp (x, 1024, -10, GMP_RNDN);
if (mpfr_cmp_ui(x, 1))
ERROR("(1024U,-10)");
mpfr_set_ui_2exp (x, 1024, 10, GMP_RNDN);
if (mpfr_cmp_ui(x, 1024*1024))
ERROR("(1024U,+10)");
mpfr_set_si_2exp (x, -1024L * 1024L, -10, GMP_RNDN);
if (mpfr_cmp_si(x, -1024))
ERROR("(1M,-10)");
mpfr_set_ui_2exp (x, 0x92345678, 16, GMP_RNDN);
if (mpfr_cmp_str (x, "[email protected]", 16, GMP_RNDN))
ERROR("(x92345678U,+16)");
mpfr_set_si_2exp (x, -0x1ABCDEF0, -256, GMP_RNDN);
if (mpfr_cmp_str (x, "[email protected]", 16, GMP_RNDN))
ERROR("(-x1ABCDEF0,-256)");
mpfr_set_prec (x, 2);
res = mpfr_set_si_2exp (x, 7, 10, GMP_RNDU);
if (mpfr_cmp_ui (x, 1<<13) || res <= 0)
ERROR ("Prec 2 + si_2exp");
res = mpfr_set_ui_2exp (x, 7, 10, GMP_RNDU);
if (mpfr_cmp_ui (x, 1<<13) || res <= 0)
ERROR ("Prec 2 + ui_2exp");
mpfr_clear_flags ();
mpfr_set_ui_2exp (x, 17, MPFR_EMAX_MAX, GMP_RNDN);
if (!mpfr_inf_p (x) || MPFR_IS_NEG (x))
ERROR ("mpfr_set_ui_2exp and overflow (bad result)");
if (!mpfr_overflow_p ())
ERROR ("mpfr_set_ui_2exp and overflow (overflow flag not set)");
mpfr_clear_flags ();
mpfr_set_si_2exp (x, 17, MPFR_EMAX_MAX, GMP_RNDN);
if (!mpfr_inf_p (x) || MPFR_IS_NEG (x))
ERROR ("mpfr_set_si_2exp (pos) and overflow (bad result)");
if (!mpfr_overflow_p ())
ERROR ("mpfr_set_si_2exp (pos) and overflow (overflow flag not set)");
mpfr_clear_flags ();
mpfr_set_si_2exp (x, -17, MPFR_EMAX_MAX, GMP_RNDN);
if (!mpfr_inf_p (x) || MPFR_IS_POS (x))
ERROR ("mpfr_set_si_2exp (neg) and overflow (bad result)");
if (!mpfr_overflow_p ())
ERROR ("mpfr_set_si_2exp (neg) and overflow (overflow flag not set)");
mpfr_clear (x);
}
int
mpfr_digamma (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
int inex;
MPFR_SAVE_EXPO_DECL (expo);
MPFR_LOG_FUNC
(("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode),
("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y, inex));
if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)))
{
if (MPFR_IS_NAN(x))
{
MPFR_SET_NAN(y);
MPFR_RET_NAN;
}
else if (MPFR_IS_INF(x))
{
if (MPFR_IS_POS(x)) /* Digamma(+Inf) = +Inf */
{
MPFR_SET_SAME_SIGN(y, x);
MPFR_SET_INF(y);
MPFR_RET(0);
}
else /* Digamma(-Inf) = NaN */
{
MPFR_SET_NAN(y);
MPFR_RET_NAN;
}
}
else /* Zero case */
{
/* the following works also in case of overlap */
MPFR_SET_INF(y);
MPFR_SET_OPPOSITE_SIGN(y, x);
mpfr_set_divby0 ();
MPFR_RET(0);
}
}
/* Digamma is undefined for negative integers */
if (MPFR_IS_NEG(x) && mpfr_integer_p (x))
{
MPFR_SET_NAN(y);
MPFR_RET_NAN;
}
/* now x is a normal number */
MPFR_SAVE_EXPO_MARK (expo);
/* for x very small, we have Digamma(x) = -1/x - gamma + O(x), more precisely
-1 < Digamma(x) + 1/x < 0 for -0.2 < x < 0.2, thus:
(i) either x is a power of two, then 1/x is exactly representable, and
as long as 1/2*ulp(1/x) > 1, we can conclude;
(ii) otherwise assume x has <= n bits, and y has <= n+1 bits, then
|y + 1/x| >= 2^(-2n) ufp(y), where ufp means unit in first place.
Since |Digamma(x) + 1/x| <= 1, if 2^(-2n) ufp(y) >= 2, then
|y - Digamma(x)| >= 2^(-2n-1)ufp(y), and rounding -1/x gives the correct result.
If x < 2^E, then y > 2^(-E), thus ufp(y) > 2^(-E-1).
A sufficient condition is thus EXP(x) <= -2 MAX(PREC(x),PREC(Y)). */
if (MPFR_EXP(x) < -2)
{
if (MPFR_EXP(x) <= -2 * (mpfr_exp_t) MAX(MPFR_PREC(x), MPFR_PREC(y)))
{
int signx = MPFR_SIGN(x);
inex = mpfr_si_div (y, -1, x, rnd_mode);
if (inex == 0) /* x is a power of two */
{ /* result always -1/x, except when rounding down */
if (rnd_mode == MPFR_RNDA)
rnd_mode = (signx > 0) ? MPFR_RNDD : MPFR_RNDU;
if (rnd_mode == MPFR_RNDZ)
rnd_mode = (signx > 0) ? MPFR_RNDU : MPFR_RNDD;
if (rnd_mode == MPFR_RNDU)
inex = 1;
else if (rnd_mode == MPFR_RNDD)
{
mpfr_nextbelow (y);
inex = -1;
}
else /* nearest */
inex = 1;
}
MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
goto end;
}
}
if (MPFR_IS_NEG(x))
inex = mpfr_digamma_reflection (y, x, rnd_mode);
/* if x < 1/2 we use the reflection formula */
else if (MPFR_EXP(x) < 0)
inex = mpfr_digamma_reflection (y, x, rnd_mode);
else
inex = mpfr_digamma_positive (y, x, rnd_mode);
end:
MPFR_SAVE_EXPO_FREE (expo);
return mpfr_check_range (y, inex, rnd_mode);
//.........这里部分代码省略.........
开发者ID:epowers，项目名称:mpfr，代码行数:101，代码来源:digamma.c
示例10: special
static void
special (void)
{
mpfr_t x, y;
mpq_t q;
mpz_t z;
int res = 0;
mpfr_init (x);
mpfr_init (y);
mpq_init (q);
mpz_init (z);
/* cancellation in mpfr_add_q */
mpfr_set_prec (x, 60);
mpfr_set_prec (y, 20);
mpz_set_str (mpq_numref (q), "-187207494", 10);
mpz_set_str (mpq_denref (q), "5721", 10);
mpfr_set_str_binary (x, "11111111101001011011100101100011011110010011100010000100001E-44");
mpfr_add_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("cancelation in add_q", mpfr_cmp_ui_2exp (y, 256783, -64) == 0);
mpfr_set_prec (x, 19);
mpfr_set_str_binary (x, "0.1011110101110011100E0");
mpz_set_str (mpq_numref (q), "187207494", 10);
mpz_set_str (mpq_denref (q), "5721", 10);
mpfr_set_prec (y, 29);
mpfr_add_q (y, x, q, MPFR_RNDD);
mpfr_set_prec (x, 29);
mpfr_set_str_binary (x, "11111111101001110011010001001E-14");
CHECK_FOR ("cancelation in add_q", mpfr_cmp (x,y) == 0);
/* Inf */
mpfr_set_inf (x, 1);
mpz_set_str (mpq_numref (q), "395877315", 10);
mpz_set_str (mpq_denref (q), "3508975966", 10);
mpfr_set_prec (y, 118);
mpfr_add_q (y, x, q, MPFR_RNDU);
CHECK_FOR ("inf", mpfr_inf_p (y) && mpfr_sgn (y) > 0);
mpfr_sub_q (y, x, q, MPFR_RNDU);
CHECK_FOR ("inf", mpfr_inf_p (y) && mpfr_sgn (y) > 0);
/* Nan */
MPFR_SET_NAN (x);
mpfr_add_q (y, x, q, MPFR_RNDU);
CHECK_FOR ("nan", mpfr_nan_p (y));
mpfr_sub_q (y, x, q, MPFR_RNDU);
CHECK_FOR ("nan", mpfr_nan_p (y));
/* Exact value */
mpfr_set_prec (x, 60);
mpfr_set_prec (y, 60);
mpfr_set_str1 (x, "0.5");
mpz_set_str (mpq_numref (q), "3", 10);
mpz_set_str (mpq_denref (q), "2", 10);
res = mpfr_add_q (y, x, q, MPFR_RNDU);
CHECK_FOR ("0.5+3/2", mpfr_cmp_ui(y, 2)==0 && res==0);
res = mpfr_sub_q (y, x, q, MPFR_RNDU);
CHECK_FOR ("0.5-3/2", mpfr_cmp_si(y, -1)==0 && res==0);
/* Inf Rationnal */
mpq_set_ui (q, 1, 0);
mpfr_set_str1 (x, "0.5");
res = mpfr_add_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("0.5+1/0", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0);
res = mpfr_sub_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("0.5-1/0", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0);
mpq_set_si (q, -1, 0);
res = mpfr_add_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("0.5+ -1/0", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0);
res = mpfr_sub_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("0.5- -1/0", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0);
res = mpfr_div_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("0.5 / (-1/0)", mpfr_zero_p (y) && MPFR_IS_NEG (y) && res == 0);
mpq_set_ui (q, 1, 0);
mpfr_set_inf (x, 1);
res = mpfr_add_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("+Inf + +Inf", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0);
res = mpfr_sub_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("+Inf - +Inf", MPFR_IS_NAN (y) && res == 0);
mpfr_set_inf (x, -1);
res = mpfr_add_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("-Inf + +Inf", MPFR_IS_NAN (y) && res == 0);
res = mpfr_sub_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("-Inf - +Inf", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0);
mpq_set_si (q, -1, 0);
mpfr_set_inf (x, 1);
res = mpfr_add_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("+Inf + -Inf", MPFR_IS_NAN (y) && res == 0);
res = mpfr_sub_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("+Inf - -Inf", mpfr_inf_p (y) && MPFR_IS_POS (y) && res == 0);
mpfr_set_inf (x, -1);
res = mpfr_add_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("-Inf + -Inf", mpfr_inf_p (y) && MPFR_IS_NEG (y) && res == 0);
res = mpfr_sub_q (y, x, q, MPFR_RNDN);
CHECK_FOR ("-Inf - -Inf", MPFR_IS_NAN (y) && res == 0);
/* 0 */
mpq_set_ui (q, 0, 1);
mpfr_set_ui (x, 42, MPFR_RNDN);
//.........这里部分代码省略.........
客服电话
APP下载
官方微信
请发表评论