summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/sin.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/sin.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/sin.c32
1 files changed, 22 insertions, 10 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/sin.c b/Build/source/libs/mpfr/mpfr-src/src/sin.c
index 860aa649c09..ea920c47788 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/sin.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/sin.c
@@ -38,7 +38,7 @@ 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_exp_t expx, err1, err;
mpfr_prec_t precy, m;
int inexact, sign, reduce;
MPFR_ZIV_DECL (loop);
@@ -55,7 +55,6 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
MPFR_SET_NAN (y);
MPFR_RET_NAN;
-
}
else /* x is zero */
{
@@ -66,9 +65,11 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
}
}
+ expx = MPFR_GET_EXP (x);
+ err1 = -2 * expx;
+
/* 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_FAST_COMPUTE_IF_SMALL_INPUT (y, x, err1, 2, 0, rnd_mode, {});
MPFR_SAVE_EXPO_MARK (expo);
@@ -81,8 +82,20 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
goto end;
}
- m = precy + MPFR_INT_CEIL_LOG2 (precy) + 13;
- expx = MPFR_GET_EXP (x);
+ m = precy + MPFR_INT_CEIL_LOG2 (precy) + 7;
+
+ /* since we compute sin(x) as sqrt(1-cos(x)^2), and for x small we have
+ cos(x)^2 ~ 1 - x^2, when subtracting cos(x)^2 from 1 we will lose
+ about -2*expx bits if expx < 0 */
+ if (expx < 0)
+ {
+ /* The following assertion includes a check for integer overflow.
+ At this point, precy < MPFR_SINCOS_THRESHOLD, so that both m and
+ err1 should be small enough. But the assertion makes the code
+ safer (a smart compiler might be able to remove it). */
+ MPFR_ASSERTN (err1 <= MPFR_PREC_MAX - m);
+ m += err1;
+ }
mpfr_init (c);
mpfr_init (xr);
@@ -113,7 +126,7 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
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)
+ if (MPFR_IS_POS (xr))
mpfr_sub (c, c, xr, MPFR_RNDZ);
else
mpfr_add (c, c, xr, MPFR_RNDZ);
@@ -135,9 +148,8 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
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_cos (c, xx, MPFR_RNDA); /* 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))