hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

use the sqrt trick in exp

Browse source
This commit is contained in:
Fredrik Johansson 2014-10-24 19:02:06 +02:00
parent 19c4aa3cb6
commit b8705b7dd5
5 changed files with 47 additions and 21 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
arb.h
View file

@ -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);

38
arb/exp.c
View file

@ -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)
{

4
arb/sin_cos.c
View file

@ -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);

21
arb/sin_cos_taylor_rs.c
View file

@ -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--;
}

2
arb/test/t-sin_cos_taylor_rf.c
View file

@ -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);