diff options
Diffstat (limited to 'src/sin_cos.c')
| -rw-r--r-- | src/sin_cos.c | 24 |
1 files changed, 21 insertions, 3 deletions
diff --git a/src/sin_cos.c b/src/sin_cos.c index 0bb00cf..7051708 100644 --- a/src/sin_cos.c +++ b/src/sin_cos.c @@ -299,11 +299,28 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_t s, c, sh, ch, sch, csh; mpfr_prec_t prec; int ok; - int inex_re, inex_im, inex_sin, inex_cos; + int inex_re, inex_im, inex_sin, inex_cos, loop = 0; prec = 2; if (rop_sin != NULL) - prec = MPC_MAX (prec, MPC_MAX_PREC (rop_sin)); + { + mp_prec_t er, ei; + prec = MPC_MAX (prec, MPC_MAX_PREC (rop_sin)); + /* since the Taylor expansion of sin(x) at x=0 is x - x^3/6 + O(x^5), + if x <= 2^(-p), then the second term/x is about 2^(-2p)/6, thus we + need at least 2p+3 bits of precision. This is true only when x is + exactly representable in the target precision. */ + if (MPC_MAX_PREC (op) <= prec) + { + er = mpfr_get_exp (mpc_realref (op)); + ei = mpfr_get_exp (mpc_imagref (op)); + /* consider the maximal exponent only */ + er = (er < ei) ? ei : er; + if (er < 0) + if (prec < 2 * (mp_prec_t) (-er) + 3) + prec = 2 * (mp_prec_t) (-er) + 3; + } + } if (rop_cos != NULL) prec = MPC_MAX (prec, MPC_MAX_PREC (rop_cos)); @@ -315,8 +332,9 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_init2 (csh, 2); do { + loop ++; ok = 1; - prec += mpc_ceil_log2 (prec) + 5; + prec += (loop <= 2) ? mpc_ceil_log2 (prec) + 5 : prec / 2; mpfr_set_prec (s, prec); mpfr_set_prec (c, prec); |