diff options
author | Akira Kakuto <kakuto@fuk.kindai.ac.jp> | 2017-12-26 11:50:03 +0000 |
---|---|---|
committer | Akira Kakuto <kakuto@fuk.kindai.ac.jp> | 2017-12-26 11:50:03 +0000 |
commit | 6ddc203e5368cf1d4108e64c18b0bf6c7d53d176 (patch) | |
tree | dda1d9b65325699c60418cfcd02fa5d800cde1f8 /Build/source/libs/mpfr/mpfr-src/src/sin_cos.c | |
parent | bbfe0bab9260ac7a06a69f710b053dcc73cecad2 (diff) |
mpfr-4.0.0
git-svn-id: svn://tug.org/texlive/trunk@46140 c570f23f-e606-0410-a88d-b1316a301751
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.c | 103 |
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); } |