summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/sin_cos.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/sin_cos.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/sin_cos.c103
1 files changed, 52 insertions, 51 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/sin_cos.c b/Build/source/libs/mpfr/mpfr-src/src/sin_cos.c
index a7c44620b46..d2fd71e0d90 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/sin_cos.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/sin_cos.c
@@ -147,7 +147,7 @@ mpfr_sin_cos (mpfr_ptr y, mpfr_ptr z, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
mpfr_mul_2ui (c, c, 1, MPFR_RNDN);
mpfr_remainder (xr, x, c, MPFR_RNDN);
mpfr_div_2ui (c, c, 1, MPFR_RNDN);
- if (MPFR_SIGN (xr) > 0)
+ if (MPFR_IS_POS (xr))
mpfr_sub (c, c, xr, MPFR_RNDZ);
else
mpfr_add (c, c, xr, MPFR_RNDZ);
@@ -174,11 +174,10 @@ mpfr_sin_cos (mpfr_ptr y, mpfr_ptr z, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
err = m;
else
err = MPFR_GET_EXP (c) + (mpfr_exp_t) (m - 3);
- if (!mpfr_can_round (c, err, MPFR_RNDN, MPFR_RNDZ,
- MPFR_PREC (z) + (rnd_mode == MPFR_RNDN)))
+ if (!MPFR_CAN_ROUND (c, err, MPFR_PREC (z), rnd_mode))
goto next_step;
- /* we can't set z now, because in case z = x, and the mpfr_can_round()
+ /* we can't set z now, because in case z = x, and the MPFR_CAN_ROUND()
call below fails, we will have clobbered the input */
mpfr_set_prec (xr, MPFR_PREC(c));
mpfr_swap (xr, c); /* save the approximation of the cosine in xr */
@@ -197,8 +196,7 @@ mpfr_sin_cos (mpfr_ptr y, mpfr_ptr z, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
/* the absolute error on c is at most 2^(err-m), which we must put
in the form 2^(EXP(c)-err). */
err = MPFR_GET_EXP (c) + (mpfr_exp_t) m - err;
- if (mpfr_can_round (c, err, MPFR_RNDN, MPFR_RNDZ,
- MPFR_PREC (y) + (rnd_mode == MPFR_RNDN)))
+ if (MPFR_CAN_ROUND (c, err, MPFR_PREC (y), rnd_mode))
break;
/* check for huge cancellation */
if (err < (mpfr_exp_t) MPFR_PREC (y))
@@ -236,8 +234,9 @@ mpfr_sin_cos (mpfr_ptr y, mpfr_ptr z, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
static mpfr_prec_t
reduce (mpz_t Q, mpz_srcptr R, mpfr_prec_t prec)
{
- mpfr_prec_t l = mpz_sizeinbase (R, 2);
+ mpfr_prec_t l;
+ MPFR_MPZ_SIZEINBASE2(l, R);
l = (l > prec) ? l - prec : 0;
mpz_fdiv_q_2exp (Q, R, l);
return l;
@@ -249,10 +248,13 @@ reduce (mpz_t Q, mpz_srcptr R, mpfr_prec_t prec)
static unsigned long
reduce2 (mpz_t S, mpz_t C, mpfr_prec_t prec)
{
- unsigned long ls = mpz_sizeinbase (S, 2);
- unsigned long lc = mpz_sizeinbase (C, 2);
+ unsigned long ls;
+ unsigned long lc;
unsigned long l;
+ MPFR_MPZ_SIZEINBASE2(ls, S);
+ MPFR_MPZ_SIZEINBASE2(lc, C);
+
l = (ls < lc) ? ls : lc; /* smaller length */
l = (l > prec) ? l - prec : 0;
mpz_fdiv_q_2exp (S, S, l);
@@ -287,9 +289,9 @@ sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r,
mpz_t T[GMP_NUMB_BITS], Q[GMP_NUMB_BITS], ptoj[GMP_NUMB_BITS], pp;
mpfr_prec_t log2_nb_terms[GMP_NUMB_BITS], mult[GMP_NUMB_BITS];
mpfr_prec_t accu[GMP_NUMB_BITS], size_ptoj[GMP_NUMB_BITS];
- mpfr_prec_t prec_i_have, r0 = r;
- unsigned long alloc, i, j, k;
- mpfr_prec_t l;
+ mpfr_prec_t prec_i_have, h, r0 = r, pp_s, p_s;
+ unsigned long i, j, m;
+ int alloc, k, l;
if (MPFR_UNLIKELY(mpz_cmp_ui (p, 0) == 0)) /* sin(x)/x -> 1 */
{
@@ -305,10 +307,10 @@ sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r,
mpz_init (pp);
/* normalize p (non-zero here) */
- l = mpz_scan1 (p, 0);
- mpz_fdiv_q_2exp (pp, p, l); /* p = pp * 2^l */
+ h = mpz_scan1 (p, 0);
+ mpz_fdiv_q_2exp (pp, p, h); /* p = pp * 2^h */
mpz_mul (pp, pp, pp);
- r = 2 * (r - l); /* x^2 = (p/2^r0)^2 = pp / 2^r */
+ r = 2 * (r - h); /* x^2 = (p/2^r0)^2 = pp / 2^r */
/* now p is odd */
alloc = 2;
@@ -319,14 +321,16 @@ sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r,
mpz_init (Q[1]);
mpz_init (ptoj[1]);
mpz_mul (ptoj[1], pp, pp); /* ptoj[1] = pp^2 */
- size_ptoj[1] = mpz_sizeinbase (ptoj[1], 2);
+ MPFR_MPZ_SIZEINBASE2(size_ptoj[1], ptoj[1]);
mpz_mul_2exp (T[0], T[0], r);
mpz_sub (T[0], T[0], pp); /* 6*2^r - pp = 6*2^r*(1 - x^2/6) */
log2_nb_terms[0] = 1;
/* already take into account the factor x=p/2^r in sin(x) = x * (...) */
- mult[0] = r - mpz_sizeinbase (pp, 2) + r0 - mpz_sizeinbase (p, 2);
+ MPFR_MPZ_SIZEINBASE2(pp_s, pp);
+ MPFR_MPZ_SIZEINBASE2(p_s, p);
+ mult[0] = r - pp_s + r0 - p_s;
/* we have x^3 < 1/2^mult[0] */
for (i = 2, k = 0, prec_i_have = mult[0]; prec_i_have < prec; i += 2)
@@ -343,7 +347,7 @@ sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r,
mpz_init (Q[k+1]);
mpz_init (ptoj[k+1]);
mpz_mul (ptoj[k+1], ptoj[k], ptoj[k]); /* pp^(2^(k+1)) */
- size_ptoj[k+1] = mpz_sizeinbase (ptoj[k+1], 2);
+ MPFR_MPZ_SIZEINBASE2(size_ptoj[k+1], ptoj[k+1]);
}
/* for i even, we have Q[k] = (2*i)*(2*i+1), T[k] = 1,
then Q[k+1] = (2*i+2)*(2*i+3), T[k+1] = 1,
@@ -358,7 +362,8 @@ sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r,
mpz_mul_ui (Q[k], Q[k], 2 * i + 1);
/* the next term of the series is divided by Q[k] and multiplied
by pp^2/2^(2r), thus the mult. factor < 1/2^mult[k] */
- mult[k] = mpz_sizeinbase (Q[k], 2) + 2 * r - size_ptoj[1] - 1;
+ MPFR_MPZ_SIZEINBASE2(mult[k], Q[k]);
+ mult[k] += 2 * r - size_ptoj[1] - 1;
/* the absolute contribution of the next term is 1/2^accu[k] */
accu[k] = (k == 0) ? mult[k] : mult[k] + accu[k-1];
prec_i_have = accu[k]; /* the current term is < 1/2^accu[k] */
@@ -373,7 +378,7 @@ sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r,
mpz_mul (Q[k-1], Q[k-1], Q[k]);
log2_nb_terms[k-1] ++; /* number of terms in S[k-1]
is a power of 2 by construction */
- prec_i_have = mpz_sizeinbase (Q[k], 2);
+ MPFR_MPZ_SIZEINBASE2(prec_i_have, Q[k]);
mult[k-1] += prec_i_have + (r << l) - size_ptoj[l] - 1;
accu[k-1] = (k == 1) ? mult[k-1] : mult[k-1] + accu[k-2];
prec_i_have = accu[k-1];
@@ -385,46 +390,45 @@ sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r,
/* accumulate all products in T[0] and Q[0]. Warning: contrary to above,
here we do not have log2_nb_terms[k-1] = log2_nb_terms[k]+1. */
- l = 0; /* number of accumulated terms in the right part T[k]/Q[k] */
+ h = 0; /* number of accumulated terms in the right part T[k]/Q[k] */
while (k > 0)
{
- j = log2_nb_terms[k-1];
- mpz_mul (T[k], T[k], ptoj[j]);
+ mpz_mul (T[k], T[k], ptoj[log2_nb_terms[k-1]]);
mpz_mul (T[k-1], T[k-1], Q[k]);
- l += 1 << log2_nb_terms[k];
- mpz_mul_2exp (T[k-1], T[k-1], r * l);
+ h += (mpfr_prec_t) 1 << log2_nb_terms[k];
+ mpz_mul_2exp (T[k-1], T[k-1], r * h);
mpz_add (T[k-1], T[k-1], T[k]);
mpz_mul (Q[k-1], Q[k-1], Q[k]);
k--;
}
- l = r0 + r * (i - 1); /* implicit multiplier 2^r for Q0 */
- /* at this point T[0]/(2^l*Q[0]) is an approximation of sin(x) where the 1st
+ m = r0 + r * (i - 1); /* implicit multiplier 2^r for Q0 */
+ /* at this point T[0]/(2^m*Q[0]) is an approximation of sin(x) where the 1st
neglected term has contribution < 1/2^prec, thus since the series has
alternate signs, the error is < 1/2^prec */
/* we truncate Q0 to prec bits: the relative error is at most 2^(1-prec),
which means that Q0 = Q[0] * (1+theta) with |theta| <= 2^(1-prec)
[up to a power of two] */
- l += reduce (Q0, Q[0], prec);
- l -= reduce (T[0], T[0], prec);
- /* multiply by x = p/2^l */
+ m += reduce (Q0, Q[0], prec);
+ m -= reduce (T[0], T[0], prec);
+ /* multiply by x = p/2^m */
mpz_mul (S0, T[0], p);
- l -= reduce (S0, S0, prec); /* S0 = T[0] * (1 + theta)^2 up to power of 2 */
+ m -= reduce (S0, S0, prec); /* S0 = T[0] * (1 + theta)^2 up to power of 2 */
/* sin(X) ~ S0/Q0*(1 + theta)^3 + err with |theta| <= 2^(1-prec) and
|err| <= 2^(-prec), thus since |S0/Q0| <= 1:
|sin(X) - S0/Q0| <= 4*|theta*S0/Q0| + |err| <= 9*2^(-prec) */
mpz_clear (pp);
- for (j = 0; j < alloc; j ++)
+ for (k = 0; k < alloc; k ++)
{
- mpz_clear (T[j]);
- mpz_clear (Q[j]);
- mpz_clear (ptoj[j]);
+ mpz_clear (T[k]);
+ mpz_clear (Q[k]);
+ mpz_clear (ptoj[k]);
}
/* compute cos(X) from sin(X): sqrt(1-(S/Q)^2) = sqrt(Q^2-S^2)/Q
- = sqrt(Q0^2*2^(2l)-S0^2)/Q0.
+ = sqrt(Q0^2*2^(2m)-S0^2)/Q0.
Write S/Q = sin(X) + eps with |eps| <= 9*2^(-prec),
then sqrt(Q^2-S^2) = sqrt(Q^2-Q^2*(sin(X)+eps)^2)
= sqrt(Q^2*cos(X)^2-Q^2*(2*sin(X)*eps+eps^2))
@@ -438,14 +442,14 @@ sin_bs_aux (mpz_t Q0, mpz_t S0, mpz_t C0, mpz_srcptr p, mpfr_prec_t r,
= Q*cos(X)*(1+eps3+eps2/(Q*cos(X)))
= Q*cos(X)*(1+eps4) with |eps4| <= 9*2^(-prec)
since |Q| >= 2^(prec-1) */
- /* we assume that Q0*2^l >= 2^(prec-1) */
- MPFR_ASSERTN(l + mpz_sizeinbase (Q0, 2) >= prec);
+ /* we assume that Q0*2^m >= 2^(prec-1) */
+ MPFR_ASSERTN(m + mpz_sizeinbase (Q0, 2) >= prec);
mpz_mul (C0, Q0, Q0);
- mpz_mul_2exp (C0, C0, 2 * l);
+ mpz_mul_2exp (C0, C0, 2 * m);
mpz_submul (C0, S0, S0);
mpz_sqrt (C0, C0);
- return l;
+ return m;
}
/* Put in s and c approximations of sin(x) and cos(x) respectively.
@@ -499,7 +503,7 @@ sincos_aux (mpfr_t s, mpfr_t c, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
/* y <- trunc(x2 * 2^sh) = trunc(x * 2^(2*sh-1)) */
mpfr_mul_2exp (x2, x2, sh, MPFR_RNDN); /* exact */
- mpfr_get_z (y, x2, MPFR_RNDZ); /* round towards zero: now
+ mpfr_get_z (y, x2, MPFR_RNDZ); /* round toward zero: now
0 <= x2 < 2^sh, thus
0 <= x2/2^(sh-1) < 2^(1-sh) */
if (mpz_cmp_ui (y, 0) == 0)
@@ -573,6 +577,7 @@ mpfr_sincos_fast (mpfr_t s, mpfr_t c, mpfr_srcptr x, mpfr_rnd_t rnd)
mpfr_t x_red, ts, tc;
mpfr_prec_t w;
mpfr_exp_t err, errs, errc;
+ MPFR_GROUP_DECL (group);
MPFR_ZIV_DECL (loop);
MPFR_ASSERTN(s != c);
@@ -583,8 +588,8 @@ mpfr_sincos_fast (mpfr_t s, mpfr_t c, mpfr_srcptr x, mpfr_rnd_t rnd)
else
w = MPFR_PREC(s) >= MPFR_PREC(c) ? MPFR_PREC(s) : MPFR_PREC(c);
w += MPFR_INT_CEIL_LOG2(w) + 9; /* ensures w >= 10 (needed by sincos_aux) */
- mpfr_init2 (ts, w);
- mpfr_init2 (tc, w);
+
+ MPFR_GROUP_INIT_2(group, w, ts, tc);
MPFR_ZIV_INIT (loop, w);
for (;;)
@@ -597,11 +602,9 @@ mpfr_sincos_fast (mpfr_t s, mpfr_t c, mpfr_srcptr x, mpfr_rnd_t rnd)
/* if -Pi/4 <= x < 0, use sin(-x)=-sin(x) */
else if (MPFR_IS_NEG(x) && mpfr_cmp_si_2exp (x, -1686629713, -31) >= 0)
{
- mpfr_init2 (x_red, MPFR_PREC(x));
- mpfr_neg (x_red, x, rnd); /* exact */
+ MPFR_ALIAS(x_red, x, MPFR_SIGN_POS, MPFR_GET_EXP(x));
err = sincos_aux (ts, tc, x_red, MPFR_RNDN);
- mpfr_neg (ts, ts, MPFR_RNDN);
- mpfr_clear (x_red);
+ MPFR_CHANGE_SIGN(ts);
}
else /* argument reduction is needed */
{
@@ -650,15 +653,13 @@ mpfr_sincos_fast (mpfr_t s, mpfr_t c, mpfr_srcptr x, mpfr_rnd_t rnd)
(c == NULL || MPFR_CAN_ROUND (tc, w - errc, MPFR_PREC(c), rnd)))
break;
MPFR_ZIV_NEXT (loop, w);
- mpfr_set_prec (ts, w);
- mpfr_set_prec (tc, w);
+ MPFR_GROUP_REPREC_2(group, w, ts, tc);
}
MPFR_ZIV_FREE (loop);
inexs = (s == NULL) ? 0 : mpfr_set (s, ts, rnd);
inexc = (c == NULL) ? 0 : mpfr_set (c, tc, rnd);
- mpfr_clear (ts);
- mpfr_clear (tc);
+ MPFR_GROUP_CLEAR (group);
return INEX(inexs,inexc);
}