From b8705b7dd5d55104ada4c36e6128d8f3a1f69941 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Fri, 24 Oct 2014 19:02:06 +0200 Subject: [PATCH] use the sqrt trick in exp --- arb.h | 3 ++- arb/exp.c | 38 ++++++++++++++++++++++++++++------ arb/sin_cos.c | 4 ++-- arb/sin_cos_taylor_rs.c | 21 +++++++++---------- arb/test/t-sin_cos_taylor_rf.c | 2 +- 5 files changed, 47 insertions(+), 21 deletions(-) diff --git a/arb.h b/arb.h index ae40d7d7..16b93e86 100644 --- a/arb.h +++ b/arb.h @@ -983,7 +983,8 @@ void _arb_sin_cos_taylor_naive(mp_ptr ysin, mp_ptr ycos, mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N); void _arb_sin_cos_taylor_rs(mp_ptr ysin, mp_ptr ycos, - mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N, int sinonly); + mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N, + int sinonly, int alternating); int _arb_get_mpn_fixed_mod_pi4(mp_ptr w, fmpz_t q, int * octant, mp_limb_t * error, const arf_t x, mp_size_t wn); diff --git a/arb/exp.c b/arb/exp.c index 6c6c450f..0a130024 100644 --- a/arb/exp.c +++ b/arb/exp.c @@ -259,14 +259,40 @@ arb_exp_arf(arb_t z, const arf_t x, long prec, int minus_one) error is <= 2^-wp */ N = _arb_exp_taylor_bound(-r, wp); - /* Evaluate Taylor series */ - _arb_exp_taylor_rs(t, &error2, w, wn, N); + if (N < 60) + { + /* Evaluate Taylor series */ + _arb_exp_taylor_rs(t, &error2, w, wn, N); + /* Taylor series evaluation error */ + error += error2; + /* Taylor series truncation error */ + error += 1UL << (wprounded-wp); + } + else /* Compute cosh(a) from sinh(a) using a square root. */ + { + /* the summation for sinh is actually done to (2N-1)! */ + N = (N + 1) / 2; - /* Taylor series evaluation error */ - error += error2; + /* Evaluate Taylor series for sinh */ + _arb_sin_cos_taylor_rs(t, u, &error2, w, wn, N, 1, 0); + error += error2; + error += 1UL << (wprounded-wp); - /* Taylor series truncation error */ - error += 1UL << (wprounded-wp); + /* 1 + sinh^2, with wn + 1 limbs */ + mpn_sqr(u, t, wn); + u[2 * wn] = 1; + + /* cosh, with wn + 1 limbs */ + mpn_sqrtrem(w, u, u, 2 * wn + 1); + + /* exp = sinh + cosh */ + t[wn] = w[wn] + mpn_add_n(t, t, w, wn); + + /* When converting sinh to cosh, the error for cosh must be + smaller than the error for sinh; but we also get 1 ulp + extra error from the square root. */ + error += 1; + } if (wp <= ARB_EXP_TAB1_PREC) { diff --git a/arb/sin_cos.c b/arb/sin_cos.c index f73a9352..82a46373 100644 --- a/arb/sin_cos.c +++ b/arb/sin_cos.c @@ -303,7 +303,7 @@ arb_sin_cos_arf_new(arb_t zsin, arb_t zcos, const arf_t x, long prec) if (N < 14) { /* Evaluate Taylor series */ - _arb_sin_cos_taylor_rs(sina, cosa, &error2, w, wn, N, 0); + _arb_sin_cos_taylor_rs(sina, cosa, &error2, w, wn, N, 0, 1); /* Taylor series evaluation error */ error += error2; /* Taylor series truncation error */ @@ -312,7 +312,7 @@ arb_sin_cos_arf_new(arb_t zsin, arb_t zcos, const arf_t x, long prec) else /* Compute cos(a) from sin(a) using a square root. */ { /* Evaluate Taylor series */ - _arb_sin_cos_taylor_rs(sina, cosa, &error2, w, wn, N, 1); + _arb_sin_cos_taylor_rs(sina, cosa, &error2, w, wn, N, 1, 1); error += error2; error += 1UL << (wprounded-wp); diff --git a/arb/sin_cos_taylor_rs.c b/arb/sin_cos_taylor_rs.c index c8bedbc5..8be8c6d7 100644 --- a/arb/sin_cos_taylor_rs.c +++ b/arb/sin_cos_taylor_rs.c @@ -33,7 +33,8 @@ extern const mp_limb_t factorial_tab_numer[FACTORIAL_TAB_SIZE]; extern const mp_limb_t factorial_tab_denom[FACTORIAL_TAB_SIZE]; void _arb_sin_cos_taylor_rs(mp_ptr ysin, mp_ptr ycos, - mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N, int sinonly) + mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N, + int sinonly, int alternating) { mp_ptr s, t, xpow; mp_limb_t new_denom, old_denom, c; @@ -109,24 +110,22 @@ void _arb_sin_cos_taylor_rs(mp_ptr ysin, mp_ptr ycos, /* change denominators */ if (new_denom != old_denom && k < N - 1) { - if (k % 2 == 0) - /* mpn_neg_n(s, s, xn + 1); */ + if (alternating && (k % 2 == 0)) s[xn] += old_denom; mpn_divrem_1(s, 0, s, xn + 1, old_denom); - if (k % 2 == 0) + if (alternating && (k % 2 == 0)) s[xn] -= 1; - /* mpn_neg_n(s, s, xn + 1); */ } if (power == 0) { /* add c * x^0 -- only top limb is affected */ - if (k % 2 == 0) - s[xn] += c; - else + if (alternating & k) s[xn] -= c; + else + s[xn] += c; /* Outer polynomial evaluation: multiply by x^m */ if (k != 0) @@ -139,10 +138,10 @@ void _arb_sin_cos_taylor_rs(mp_ptr ysin, mp_ptr ycos, } else { - if (k % 2 == 0) - s[xn] += mpn_addmul_1(s, XPOW_READ(power), xn, c); - else + if (alternating & k) s[xn] -= mpn_submul_1(s, XPOW_READ(power), xn, c); + else + s[xn] += mpn_addmul_1(s, XPOW_READ(power), xn, c); power--; } diff --git a/arb/test/t-sin_cos_taylor_rf.c b/arb/test/t-sin_cos_taylor_rf.c index 16157e19..5b9ce489 100644 --- a/arb/test/t-sin_cos_taylor_rf.c +++ b/arb/test/t-sin_cos_taylor_rf.c @@ -63,7 +63,7 @@ int main() x[xn - 1] &= (LIMB_ONES >> 4); _arb_sin_cos_taylor_naive(y1s, y1c, &err1, x, xn, N); - _arb_sin_cos_taylor_rs(y2s, y2c, &err2, x, xn, N, 0); + _arb_sin_cos_taylor_rs(y2s, y2c, &err2, x, xn, N, 0, 1); cmp = mpn_cmp(y1s, y2s, xn);