Aren/NetBSD
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
efee5258bc
NetBSD/external/lgpl3/mpfr/dist/atan.c
mrg efee5258bc initial import of MPRF 3.0.1. 
The MPFR library is a C library for multiple-precision floating-point
computations with exact rounding (also called correct rounding).  It is
based on the GMP multiple-precision library and should replace the MPF
class in further releases of GMP.

GCC >= 4.2 requires MPFR.
2011-06-20 05:53:01 +00:00

436 lines
15 KiB
C
Raw Blame History

/* mpfr_atan -- arc-tangent of a floating-point number
Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
Contributed by the Arenaire and Cacao projects, INRIA.
This file is part of the GNU MPFR Library, and was contributed by Mathieu Dutour.
The GNU MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.
The GNU MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
#define MPFR_NEED_LONGLONG_H
#include "mpfr-impl.h"
/* 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.
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.
*/
static void
mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, long r, int m, mpz_t *tab)
{
mpz_t *S, *Q, *ptoj;
unsigned long n, i, k, j, l;
mpfr_exp_t diff, expo;
int im, done;
mpfr_prec_t mult, *accu, *log2_nb_terms;
mpfr_prec_t precy = MPFR_PREC(y);
MPFR_ASSERTD(mpz_cmp_ui (p, 0) != 0);
accu = (mpfr_prec_t*) (*__gmp_allocate_func) ((2 * m + 2) * sizeof (mpfr_prec_t));
log2_nb_terms = accu + m + 1;
/* Set Tables */
S = tab; /* S */
ptoj = S + 1*(m+1); /* p^2^j Precomputed table */
Q = S + 2*(m+1); /* Product of Odd integer table */
/* From p to p^2, and r to 2r */
mpz_mul (p, p, p);
MPFR_ASSERTD (2 * r > r);
r = 2 * r;
/* Normalize p */
n = mpz_scan1 (p, 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);
MPFR_ASSERTD (m > 0);
/* check if p=1 (special case) */
l = 0;
/*
We compute by binary splitting, with X = x^2 = p/2^r:
P(a,b) = p if a+1=b, P(a,c)*P(c,b) otherwise
Q(a,b) = (2a+1)*2^r if a+1=b [except Q(0,1)=1], Q(a,c)*Q(c,b) otherwise
S(a,b) = p*(2a+1) if a+1=b, Q(c,b)*S(a,c)+Q(a,c)*P(a,c)*S(c,b) otherwise
Then atan(x)/x ~ S(0,i)/Q(0,i) for i so that (p/2^r)^i/i is small enough.
The factor 2^(r*(b-a)) in Q(a,b) is implicit, thus we have to take it
into account when we compute with Q.
*/
accu[0] = 0; /* accu[k] = Mult[0] + ... + Mult[k], where Mult[j] is the
number of bits of the corresponding term S[j]/Q[j] */
if (mpz_cmp_ui (p, 1) != 0)
{
/* p <> 1: precompute ptoj table */
mpz_set (ptoj[0], p);
for (im = 1 ; im <= m ; im ++)
mpz_mul (ptoj[im], ptoj[im - 1], ptoj[im - 1]);
/* main loop */
n = 1UL << m;
/* 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 ++)
{
/* initialize both S[k],Q[k] and S[k+1],Q[k+1] */
mpz_set_ui (Q[k+1], 2 * i + 3); /* Q(i+1,i+2) */
mpz_mul_ui (S[k+1], p, 2 * i + 1); /* S(i+1,i+2) */
mpz_mul_2exp (S[k], Q[k+1], r);
mpz_sub (S[k], S[k], S[k+1]); /* S(i,i+2) */
mpz_mul_ui (Q[k], Q[k+1], 2 * i + 1); /* Q(i,i+2) */
log2_nb_terms[k] = 1; /* S[k]/Q[k] corresponds to 2 terms */
for (j = (i + 2) >> 1, l = 1; (j & 1) == 0; l ++, j >>= 1, k --)
{
/* invariant: S[k-1]/Q[k-1] and S[k]/Q[k] correspond
to 2^l terms each. We combine them into S[k-1]/Q[k-1] */
MPFR_ASSERTD (k > 0);
mpz_mul (S[k], S[k], Q[k-1]);
mpz_mul (S[k], S[k], ptoj[l]);
mpz_mul (S[k-1], S[k-1], Q[k]);
mpz_mul_2exp (S[k-1], S[k-1], r << l);
mpz_add (S[k-1], S[k-1], S[k]);
mpz_mul (Q[k-1], Q[k-1], Q[k]);
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)
done = 1;
}
}
}
else /* special case p=1: the ith term being X^i/(2i+1) with X=1/2^r,
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 ++)
{
mpz_set_ui (Q[k + 1], 2 * i + 3);
mpz_mul_2exp (S[k], Q[k+1], r);
mpz_sub_ui (S[k], S[k], 1 + 2 * i);
mpz_mul_ui (Q[k], Q[k + 1], 1 + 2 * i);
log2_nb_terms[k] = 1; /* S[k]/Q[k] corresponds to 2 terms */
for (j = (i + 2) >> 1, l = 1; (j & 1) == 0; l++, j >>= 1, k --)
{
MPFR_ASSERTD (k > 0);
mpz_mul (S[k], S[k], Q[k-1]);
mpz_mul (S[k-1], S[k-1], Q[k]);
mpz_mul_2exp (S[k-1], S[k-1], r << l);
mpz_add (S[k-1], S[k-1], S[k]);
mpz_mul (Q[k-1], Q[k-1], Q[k]);
log2_nb_terms[k-1] = l + 1;
}
}
}
/* 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] */
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-1], S[k-1], Q[k]);
l += 1 << log2_nb_terms[k];
mpz_mul_2exp (S[k-1], S[k-1], r * l);
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;
expo = diff;
if (diff >= 0)
mpz_tdiv_q_2exp (S[0], S[0], diff);
else
mpz_mul_2exp (S[0], S[0], -diff);
MPFR_MPZ_SIZEINBASE2 (diff, Q[0]);
diff -= precy;
expo -= diff;
if (diff >= 0)
mpz_tdiv_q_2exp (Q[0], Q[0], diff);
else
mpz_mul_2exp (Q[0], Q[0], -diff);
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));
}
int
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;
mpfr_exp_t exptol;
mpfr_prec_t prec, realprec, est_lost, lost;
unsigned long twopoweri, log2p, red;
int comparaison, inexact;
int i, n0, oldn0;
MPFR_GROUP_DECL (group);
MPFR_SAVE_EXPO_DECL (expo);
MPFR_ZIV_DECL (loop);
MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd_mode),
("atan[%#R]=%R inexact=%d", atan, atan, inexact));
/* Singular cases */
if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
{
if (MPFR_IS_NAN (x))
{
MPFR_SET_NAN (atan);
MPFR_RET_NAN;
}
else if (MPFR_IS_INF (x))
{
MPFR_SAVE_EXPO_MARK (expo);
if (MPFR_IS_POS (x)) /* arctan(+inf) = Pi/2 */
inexact = mpfr_const_pi (atan, rnd_mode);
else /* arctan(-inf) = -Pi/2 */
{
inexact = -mpfr_const_pi (atan,
MPFR_INVERT_RND (rnd_mode));
MPFR_CHANGE_SIGN (atan);
}
mpfr_div_2ui (atan, atan, 1, rnd_mode); /* exact (no exceptions) */
MPFR_SAVE_EXPO_FREE (expo);
return mpfr_check_range (atan, inexact, rnd_mode);
}
else /* x is necessarily 0 */
{
MPFR_ASSERTD (MPFR_IS_ZERO (x));
MPFR_SET_ZERO (atan);
MPFR_SET_SAME_SIGN (atan, x);
MPFR_RET (0);
}
}
/* atan(x) = x - x^3/3 + x^5/5...
so the error is < 2^(3*EXP(x)-1)
so `EXP(x)-(3*EXP(x)-1)` = -2*EXP(x)+1 */
MPFR_FAST_COMPUTE_IF_SMALL_INPUT (atan, x, -2 * MPFR_GET_EXP (x), 1, 0,
rnd_mode, {});
/* Set x_p=|x| */
MPFR_TMP_INIT_ABS (xp, x);
MPFR_SAVE_EXPO_MARK (expo);
/* Other simple case arctan(-+1)=-+pi/4 */
comparaison = mpfr_cmp_ui (xp, 1);
if (MPFR_UNLIKELY (comparaison == 0))
{
int neg = MPFR_IS_NEG (x);
inexact = mpfr_const_pi (atan, MPFR_IS_POS (x) ? rnd_mode
: MPFR_INVERT_RND (rnd_mode));
if (neg)
{
inexact = -inexact;
MPFR_CHANGE_SIGN (atan);
}
mpfr_div_2ui (atan, atan, 2, rnd_mode); /* exact (no exceptions) */
MPFR_SAVE_EXPO_FREE (expo);
return mpfr_check_range (atan, inexact, rnd_mode);
}
realprec = MPFR_PREC (atan) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (atan)) + 4;
prec = realprec + GMP_NUMB_BITS;
/* Initialisation */
mpz_init (ukz);
MPFR_GROUP_INIT_4 (group, prec, sk, tmp, tmp2, arctgt);
oldn0 = 0;
tabz = (mpz_t *) 0;
MPFR_ZIV_INIT (loop, prec);
for (;;)
{
/* First, if |x| < 1, we need to have more prec to be able to round (sup)
n0 = ceil(log(prec_requested + 2 + 1+ln(2.4)/ln(2))/log(2)) */
mpfr_prec_t sup;
sup = MPFR_GET_EXP (xp) < 0 ? 2 - MPFR_GET_EXP (xp) : 1; /* sup >= 1 */
n0 = MPFR_INT_CEIL_LOG2 ((realprec + sup) + 3);
/* since realprec >= 4, n0 >= ceil(log2(8)) >= 3, thus 3*n0 > 2 */
prec = (realprec + sup) + 1 + MPFR_INT_CEIL_LOG2 (3*n0-2);
/* the number of lost bits due to argument reduction is
9 - 2 * EXP(sk), which we estimate by 9 + 2*ceil(log2(p))
since we manage that sk < 1/p */
if (MPFR_PREC (atan) > 100)
{
log2p = MPFR_INT_CEIL_LOG2(prec) / 2 - 3;
est_lost = 9 + 2 * log2p;
prec += est_lost;
}
else
log2p = est_lost = 0; /* don't reduce the argument */
/* 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);
}
/* The mpfr_ui_div below mustn't underflow. This is guaranteed by
MPFR_SAVE_EXPO_MARK, but let's check that for maintainability. */
MPFR_ASSERTD (__gmpfr_emax <= 1 - __gmpfr_emin);
if (comparaison > 0) /* use atan(xp) = Pi/2 - atan(1/xp) */
mpfr_ui_div (sk, 1, xp, MPFR_RNDN);
else
mpfr_set (sk, xp, MPFR_RNDN);
/* now 0 < sk <= 1 */
/* Argument reduction: atan(x) = 2 atan((sqrt(1+x^2)-1)/x).
We want |sk| < k/sqrt(p) where p is the target precision. */
lost = 0;
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_add_ui (tmp, tmp, 1, MPFR_RNDN);
mpfr_sqrt (tmp, tmp, MPFR_RNDN);
mpfr_sub_ui (tmp, tmp, 1, MPFR_RNDN);
if (red == 0 && comparaison > 0)
/* use xp = 1/sk */
mpfr_mul (sk, tmp, xp, MPFR_RNDN);
else
mpfr_div (sk, tmp, sk, MPFR_RNDN);
}
/* 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.
*/
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;
}
/* Assignation */
MPFR_SET_ZERO (arctgt);
twopoweri = 1 << 0;
MPFR_ASSERTD (n0 >= 4);
for (i = 0 ; i < n0; i++)
{
if (MPFR_UNLIKELY (MPFR_IS_ZERO (sk)))
break;
/* Calculation of trunc(tmp) --> mpz */
mpfr_mul_2ui (tmp, sk, twopoweri, MPFR_RNDN);
mpfr_trunc (tmp, tmp);
if (!MPFR_IS_ZERO (tmp))
{
/* tmp = ukz*2^exptol */
exptol = mpfr_get_z_2exp (ukz, tmp);
/* since the s_k are decreasing (see algorithms.tex),
and s_0 = min(|x|, 1/|x|) < 1, we have sk < 1,
thus exptol < 0 */
MPFR_ASSERTD (exptol < 0);
mpz_tdiv_q_2exp (ukz, ukz, (unsigned long int) (-exptol));
/* since tmp is a non-zero integer, and tmp = ukzold*2^exptol,
we now have ukz = tmp, thus ukz is non-zero */
/* Calculation of arctan(Ak) */
mpfr_set_z (tmp, ukz, MPFR_RNDN);
mpfr_div_2ui (tmp, tmp, twopoweri, MPFR_RNDN);
mpfr_atan_aux (tmp2, ukz, twopoweri, n0 - i, tabz);
mpfr_mul (tmp2, tmp2, tmp, MPFR_RNDN);
/* Addition */
mpfr_add (arctgt, arctgt, tmp2, MPFR_RNDN);
/* Next iteration */
mpfr_sub (tmp2, sk, tmp, MPFR_RNDN);
mpfr_mul (sk, sk, tmp, MPFR_RNDN);
mpfr_add_ui (sk, sk, 1, MPFR_RNDN);
mpfr_div (sk, tmp2, sk, MPFR_RNDN);
}
twopoweri <<= 1;
}
/* Add last step (Arctan(sk) ~= sk */
mpfr_add (arctgt, arctgt, sk, MPFR_RNDN);
/* argument reduction */
mpfr_mul_2exp (arctgt, arctgt, red, MPFR_RNDN);
if (comparaison > 0)
{ /* atan(x) = Pi/2-atan(1/x) for x > 0 */
mpfr_const_pi (tmp, MPFR_RNDN);
mpfr_div_2ui (tmp, tmp, 1, MPFR_RNDN);
mpfr_sub (arctgt, tmp, arctgt, MPFR_RNDN);
}
MPFR_SET_POS (arctgt);
can_round:
if (MPFR_LIKELY (MPFR_CAN_ROUND (arctgt, realprec + est_lost - lost,
MPFR_PREC (atan), rnd_mode)))
break;
MPFR_ZIV_NEXT (loop, realprec);
}
MPFR_ZIV_FREE (loop);
inexact = mpfr_set4 (atan, arctgt, rnd_mode, MPFR_SIGN (x));
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);
return mpfr_check_range (arctgt, inexact, rnd_mode);
}
Reference in New Issue View Git Blame Copy Permalink