summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/atan.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/atan.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/atan.c196
1 files changed, 137 insertions, 59 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/atan.c b/Build/source/libs/mpfr/mpfr-src/src/atan.c
index c6df3c7af54..8fa4f6c178f 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/atan.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/atan.c
@@ -23,31 +23,114 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
#define MPFR_NEED_LONGLONG_H
#include "mpfr-impl.h"
+#if GMP_NUMB_BITS == 64
+/* for each pair (r,p), we store a 192-bit approximation of atan(x)/x for
+ x=p/2^r, with lowest limb first.
+ Sage code:
+ for p in range(1,2^ceil(r/2)):
+ x=p/2^r
+ l=floor(2^192*n(atan(x)/x, 300)).digits(2^64)
+ print ("{0x%x, 0x%x, 0x%x}, /"+"* (%d,%d) *"+"/") % (l[0],l[1],l[2],r,p)
+*/
+static const mp_limb_t atan_table[][3] = {
+ {0x6e141587261cdf00, 0x6fe445ecbc3a8d03, 0xed63382b0dda7b45}, /* (1,1) */
+ {0xaa7fa90388b3836b, 0x6dc79ef5f7a217e5, 0xfadbafc96406eb15}, /* (2,1) */
+ {0x319c12cf59d4b2dc, 0xcb2792dc0e2e0d51, 0xffaaddb967ef4e36}, /* (4,1) */
+ {0x8b3957d95d9ad922, 0xc897989f3e888ef7, 0xfeadd4d5617b6e32}, /* (4,2) */
+ {0xc4e6abc8af62e439, 0x4eb9bf602625f0b4, 0xfd0fcdd343cac19b}, /* (4,3) */
+ {0x7c18baeb9bc95789, 0xb12afb6b6d4f7e16, 0xffffaaaaddddb94b}, /* (8,1) */
+ {0x6856a0171a2f001a, 0x62351fbbe60af47, 0xfffeaaadddd4b968}, /* (8,2) */
+ {0x69164c094f49da06, 0xd517294f7373d07a, 0xfffd001032cb1179}, /* (8,3) */
+ {0x20ef65c10deef460, 0xe78c564015f76048, 0xfffaaadddb94d5bb}, /* (8,4) */
+ {0x3ce233aa002f0344, 0x9dd8ea342a65d4cc, 0xfff7ab27a1f32f95}, /* (8,5) */
+ {0xa37f403c7279c5cb, 0x13ab53a1c8db8497, 0xfff40103192ce74d}, /* (8,6) */
+ {0xe5a85657103c1aa8, 0xb8409e6c914191d3, 0xffefac8a9c40a26b}, /* (8,7) */
+ {0x806d0294c0db8816, 0x779d776dda8c6213, 0xffeaaddd4bb12542}, /* (8,8) */
+ {0x5545d1914ef21478, 0x3aea58d6660f5a12, 0xffe5051f0aebf73a}, /* (8,9) */
+ {0x6e47a91d015f4133, 0xc085ab6b490b7f02, 0xffdeb2787d4adac1}, /* (8,10) */
+ {0x4efc1f931f7ec9b3, 0xb7f43cd16195ef4b, 0xffd7b61702b09aad}, /* (8,11) */
+ {0xd27d1dbf55fed60d, 0xd812c11d7d473e5e, 0xffd0102cb3c1bfbe}, /* (8,12) */
+ {0xca629e927383fe97, 0x8c61aedf58e42206, 0xffc7c0f05db9d1b6}, /* (8,13) */
+ {0x4eff0b53d4e905b7, 0x28ac1e800ca31e9d, 0xffbec89d7dddd7e9}, /* (8,14) */
+ {0xb0a7931deec6fe60, 0xb46feea78588554b, 0xffb527743c8cdd8f} /* (8,15) */
+ };
+
+static void
+set_table (mpfr_t y, const mp_limb_t x[3])
+{
+ mpfr_prec_t p = MPFR_PREC(y);
+ mp_size_t n = MPFR_PREC2LIMBS(p);
+ mpfr_prec_t sh;
+ mp_limb_t *yp = MPFR_MANT(y);
+
+ MPFR_UNSIGNED_MINUS_MODULO (sh, p);
+ mpn_copyi (yp, x + 3 - n, n);
+ yp[0] &= ~MPFR_LIMB_MASK(sh);
+ MPFR_SET_EXP(y, 0);
+}
+#endif
+
/* If x = p/2^r, put in y an approximation of atan(x)/x using 2^m terms
for the series expansion, with an error of at most 1 ulp.
- Assumes |x| < 1.
+ Assumes 0 < x < 1, thus 1 <= p < 2^r.
+ More precisely, p consists of the floor(r/2) bits of the binary expansion
+ of a number 0 < s < 1:
+ * the bit of weight 2^-1 is for r=1, thus p <= 1
+ * the bit of weight 2^-2 is for r=2, thus p <= 1
+ * the two bits of weight 2^-3 and 2^-4 are for r=4, thus p <= 3
+ * more generally p < 2^(r/2).
If X=x^2, we want 1 - X/3 + X^2/5 - ... + (-1)^k*X^k/(2k+1) + ...
- Assume p is non-zero.
-
When we sum terms up to x^k/(2k+1), the denominator Q[0] is
3*5*7*...*(2k+1) ~ (2k/e)^k.
+
+ The tab[] array should have at least 3*(m+1) entries.
*/
static void
-mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, long r, int m, mpz_t *tab)
+mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, unsigned long r, int m, mpz_t *tab)
{
mpz_t *S, *Q, *ptoj;
- unsigned long n, i, k, j, l;
+ mp_bitcnt_t n, h, j; /* unsigned type, which is >= unsigned long */
mpfr_exp_t diff, expo;
- int im, done;
- mpfr_prec_t mult, *accu, *log2_nb_terms;
+ int im, i, k, l, done;
+ mpfr_prec_t mult;
+ mpfr_prec_t accu[MPFR_PREC_BITS], log2_nb_terms[MPFR_PREC_BITS];
mpfr_prec_t precy = MPFR_PREC(y);
MPFR_ASSERTD(mpz_cmp_ui (p, 0) != 0);
+ MPFR_ASSERTD (m+1 <= MPFR_PREC_BITS);
- accu = (mpfr_prec_t*) (*__gmp_allocate_func) ((2 * m + 2) * sizeof (mpfr_prec_t));
- log2_nb_terms = accu + m + 1;
+#if GMP_NUMB_BITS == 64
+ /* tabulate values for small precision and small value of r (which are the
+ most expensive to compute) */
+ if (precy <= 192)
+ {
+ switch (r)
+ {
+ case 1:
+ /* p has 1 bit: necessarily p=1 */
+ MPFR_ASSERTD(mpz_cmp_ui (p, 1) == 0);
+ set_table (y, atan_table[0]);
+ return;
+ case 2:
+ /* p has 1 bit: necessarily p=1 too */
+ MPFR_ASSERTD(mpz_cmp_ui (p, 1) == 0);
+ set_table (y, atan_table[1]);
+ return;
+ case 4:
+ /* p has at most 2 bits: 1 <= p <= 3 */
+ MPFR_ASSERTD(1 <= mpz_get_ui (p) && mpz_get_ui (p) <= 3);
+ set_table (y, atan_table[1 + mpz_get_ui (p)]);
+ return;
+ case 8:
+ /* p has at most 4 bits: 1 <= p <= 15 */
+ MPFR_ASSERTD(1 <= mpz_get_ui (p) && mpz_get_ui (p) <= 15);
+ set_table (y, atan_table[4 + mpz_get_ui (p)]);
+ return;
+ }
+ }
+#endif
/* Set Tables */
S = tab; /* S */
@@ -61,9 +144,12 @@ mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, long r, int m, mpz_t *tab)
/* Normalize p */
n = mpz_scan1 (p, 0);
- mpz_tdiv_q_2exp (p, p, n); /* exact */
- MPFR_ASSERTD (r > n);
- r -= n;
+ if (n > 0)
+ {
+ mpz_tdiv_q_2exp (p, p, n); /* exact */
+ MPFR_ASSERTD (r > n);
+ r -= n;
+ }
/* since |p/2^r| < 1, and p is a non-zero integer, necessarily r > 0 */
MPFR_ASSERTD (mpz_sgn (p) > 0);
@@ -90,6 +176,7 @@ mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, long r, int m, mpz_t *tab)
mpz_mul (ptoj[im], ptoj[im - 1], ptoj[im - 1]);
/* main loop */
n = 1UL << m;
+ MPFR_ASSERTN (n != 0); /* no overflow */
/* the ith term being X^i/(2i+1) with X=p/2^r, we can stop when
p^i/2^(r*i) < 2^(-precy), i.e. r*i > precy + log2(p^i) */
for (i = k = done = 0; (i < n) && (done == 0); i += 2, k ++)
@@ -115,7 +202,6 @@ mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, long r, int m, mpz_t *tab)
log2_nb_terms[k-1] = l + 1;
/* now S[k-1]/Q[k-1] corresponds to 2^(l+1) terms */
MPFR_MPZ_SIZEINBASE2(mult, ptoj[l+1]);
- /* FIXME: precompute bits(ptoj[l+1]) outside the loop? */
mult = (r << (l + 1)) - mult - 1;
accu[k-1] = (k == 1) ? mult : accu[k-2] + mult;
if (accu[k-1] > precy)
@@ -127,7 +213,10 @@ mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, long r, int m, mpz_t *tab)
we can stop when r*i > precy i.e. i > precy/r */
{
n = 1UL << m;
- for (i = k = 0; (i < n) && (i <= precy / r); i += 2, k ++)
+ if (precy / r <= n)
+ n = (precy / r) + 1;
+ MPFR_ASSERTN (n != 0); /* no overflow */
+ for (i = k = 0; i < n; i += 2, k ++)
{
mpz_set_ui (Q[k + 1], 2 * i + 3);
mpz_mul_2exp (S[k], Q[k+1], r);
@@ -148,22 +237,20 @@ mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, long r, int m, mpz_t *tab)
}
/* we need to combine S[0]/Q[0]...S[k-1]/Q[k-1] */
- l = 0; /* number of terms accumulated in S[k]/Q[k] */
+ h = 0; /* number of terms accumulated in S[k]/Q[k] */
while (k > 1)
{
k --;
/* combine S[k-1]/Q[k-1] and S[k]/Q[k] */
- j = log2_nb_terms[k-1];
mpz_mul (S[k], S[k], Q[k-1]);
if (mpz_cmp_ui (p, 1) != 0)
- mpz_mul (S[k], S[k], ptoj[j]);
+ mpz_mul (S[k], S[k], ptoj[log2_nb_terms[k-1]]);
mpz_mul (S[k-1], S[k-1], Q[k]);
- l += 1 << log2_nb_terms[k];
- mpz_mul_2exp (S[k-1], S[k-1], r * l);
+ h += (mp_bitcnt_t) 1 << log2_nb_terms[k];
+ mpz_mul_2exp (S[k-1], S[k-1], r * h);
mpz_add (S[k-1], S[k-1], S[k]);
mpz_mul (Q[k-1], Q[k-1], Q[k]);
}
- (*__gmp_free_func) (accu, (2 * m + 2) * sizeof (mpfr_prec_t));
MPFR_MPZ_SIZEINBASE2 (diff, S[0]);
diff -= 2 * precy;
@@ -183,7 +270,9 @@ mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, long r, int m, mpz_t *tab)
mpz_tdiv_q (S[0], S[0], Q[0]);
mpfr_set_z (y, S[0], MPFR_RNDD);
- MPFR_SET_EXP (y, MPFR_EXP(y) + expo - r * (i - 1));
+ /* TODO: Check/prove that the following expression doesn't overflow. */
+ expo = MPFR_GET_EXP (y) + expo - r * (i - 1);
+ MPFR_SET_EXP (y, expo);
}
int
@@ -191,7 +280,7 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
mpfr_t xp, arctgt, sk, tmp, tmp2;
mpz_t ukz;
- mpz_t *tabz;
+ mpz_t tabz[3*(MPFR_PREC_BITS+1)];
mpfr_exp_t exptol;
mpfr_prec_t prec, realprec, est_lost, lost;
unsigned long twopoweri, log2p, red;
@@ -270,10 +359,9 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
prec = realprec + GMP_NUMB_BITS;
/* Initialisation */
- mpz_init (ukz);
+ mpz_init2 (ukz, prec); /* ukz will need 'prec' bits below */
MPFR_GROUP_INIT_4 (group, prec, sk, tmp, tmp2, arctgt);
oldn0 = 0;
- tabz = (mpz_t *) 0;
MPFR_ZIV_INIT (loop, prec);
for (;;)
@@ -301,21 +389,13 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
/* Initialisation */
MPFR_GROUP_REPREC_4 (group, prec, sk, tmp, tmp2, arctgt);
- if (MPFR_LIKELY (oldn0 == 0))
- {
- oldn0 = 3 * (n0 + 1);
- tabz = (mpz_t *) (*__gmp_allocate_func) (oldn0 * sizeof (mpz_t));
- for (i = 0; i < oldn0; i++)
- mpz_init (tabz[i]);
- }
- else if (MPFR_UNLIKELY (oldn0 < 3 * (n0 + 1)))
- {
- tabz = (mpz_t *) (*__gmp_reallocate_func)
- (tabz, oldn0 * sizeof (mpz_t), 3 * (n0 + 1)*sizeof (mpz_t));
- for (i = oldn0; i < 3 * (n0 + 1); i++)
- mpz_init (tabz[i]);
- oldn0 = 3 * (n0 + 1);
- }
+ MPFR_ASSERTD (n0 <= MPFR_PREC_BITS);
+ /* Note: the tabz[] entries are used to get a rational approximation
+ of atan(x) to precision 'prec', thus allocating them to 'prec' bits
+ should be good enough. */
+ for (i = oldn0; i < 3 * (n0 + 1); i++)
+ mpz_init2 (tabz[i], prec);
+ oldn0 = 3 * (n0 + 1);
/* The mpfr_ui_div below mustn't underflow. This is guaranteed by
MPFR_SAVE_EXPO_MARK, but let's check that for maintainability. */
@@ -334,7 +414,7 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
for (red = 0; MPFR_GET_EXP(sk) > - (mpfr_exp_t) log2p; red ++)
{
lost = 9 - 2 * MPFR_EXP(sk);
- mpfr_mul (tmp, sk, sk, MPFR_RNDN);
+ mpfr_sqr (tmp, sk, MPFR_RNDN);
mpfr_add_ui (tmp, tmp, 1, MPFR_RNDN);
mpfr_sqrt (tmp, tmp, MPFR_RNDN);
mpfr_sub_ui (tmp, tmp, 1, MPFR_RNDN);
@@ -345,25 +425,25 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
mpfr_div (sk, tmp, sk, MPFR_RNDN);
}
- /* we started from x0 = 1/|x| if |x| > 1, and |x| otherwise, thus
+ /* We started from x0 = 1/|x| if |x| > 1, and |x| otherwise, thus
we had x0 = min(|x|, 1/|x|) <= 1, and applied 'red' times the
- argument reduction x -> (sqrt(1+x^2)-1)/x, which keeps 0 < x < 1,
- thus 0 < sk <= 1, and sk=1 can occur only if red=0 */
-
- /* If sk=1, then if |x| < 1, we have 1 - 2^(-prec-1) <= |x| < 1,
- or if |x| > 1, we have 1 - 2^(-prec-1) <= 1/|x| < 1, thus in all
- cases ||x| - 1| <= 2^(-prec), from which it follows
- |atan|x| - Pi/4| <= 2^(-prec), given the Taylor expansion
- atan(1+x) = Pi/4 + x/2 - x^2/4 + ...
- Since Pi/4 = 0.785..., the error is at most one ulp.
+ argument reduction x -> (sqrt(1+x^2)-1)/x, which keeps 0 < x <= 1 */
+
+ /* We first show that if the for-loop is executed at least once, then
+ sk < 1 after the loop. Indeed for 1/2 <= x <= 1, interval
+ arithmetic with precision 5 shows that (sqrt(1+x^2)-1)/x,
+ when evaluated with rounding to nearest, gives a value <= 0.875.
+ Now assume 2^(-k-1) <= x <= 2^(-k) for k >= 1.
+ Then o(x^2) <= 2^(-2k), o(1+x^2) <= 1+2^(-2k),
+ o(sqrt(1+x^2)) <= 1+2^(-2k-1), o(sqrt(1+x^2)-1) <= 2^(-2k-1),
+ and o((sqrt(1+x^2)-1)/x) <= 2^(-k) <= 1/2.
+
+ Now if sk=1 before the loop, then EXP(sk)=1 and since log2p >= 0,
+ the loop is performed at least once, thus the case sk=1 cannot
+ happen below.
*/
- if (MPFR_UNLIKELY(mpfr_cmp_ui (sk, 1) == 0))
- {
- mpfr_const_pi (arctgt, MPFR_RNDN); /* 1/2 ulp extra error */
- mpfr_div_2ui (arctgt, arctgt, 2, MPFR_RNDN); /* exact */
- realprec = prec - 2;
- goto can_round;
- }
+
+ MPFR_ASSERTD(mpfr_cmp_ui (sk, 1) < 0);
/* Assignation */
MPFR_SET_ZERO (arctgt);
@@ -416,7 +496,6 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
}
MPFR_SET_POS (arctgt);
- can_round:
if (MPFR_LIKELY (MPFR_CAN_ROUND (arctgt, realprec + est_lost - lost,
MPFR_PREC (atan), rnd_mode)))
break;
@@ -429,7 +508,6 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
for (i = 0 ; i < oldn0 ; i++)
mpz_clear (tabz[i]);
mpz_clear (ukz);
- (*__gmp_free_func) (tabz, oldn0 * sizeof (mpz_t));
MPFR_GROUP_CLEAR (group);
MPFR_SAVE_EXPO_FREE (expo);