contrib/mpfr/atan.c

   1 /* mpfr_atan -- arc-tangent of a floating-point number
   2
   3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library, and was contributed by Mathieu Dutour.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 2.1 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 /* If x = p/2^r, put in y an approximation of atan(x)/x using 2^m terms
  27    for the series expansion, with an error of at most 1 ulp.
  28    Assumes |x| < 1.
  29
  30    If X=x^2, we want 1 - X/3 + X^2/5 - ... + (-1)^k*X^k/(2k+1) + ...
  31
  32    Assume p is non-zero.
  33 */
  34 static void
  35 mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, long r, int m, mpz_t *tab)
  36 {
  37   mpz_t *S, *Q, *ptoj;
  38   unsigned long n, i, k, j, l;
  39   mp_exp_t diff, expo;
  40   int im, done;
  41   mp_prec_t mult, *accu, *log2_nb_terms;
  42   mp_prec_t precy = MPFR_PREC(y);
  43
  44   MPFR_ASSERTD(mpz_cmp_ui (p, 0) != 0);
  45
  46   accu = (mp_prec_t*) (*__gmp_allocate_func) ((2 * m + 2) * sizeof (mp_prec_t));
  47   log2_nb_terms = accu + m + 1;
  48
  49   /* Set Tables */
  50   S    = tab;           /* S */
  51   ptoj = S + 1*(m+1);   /* p^2^j Precomputed table */
  52   Q    = S + 2*(m+1);   /* Product of Odd integer  table  */
  53
  54   /* From p to p^2, and r to 2r */
  55   mpz_mul (p, p, p);
  56   MPFR_ASSERTD (2 * r > r);
  57   r = 2 * r;
  58
  59   /* Normalize p */
  60   n = mpz_scan1 (p, 0);
  61   mpz_tdiv_q_2exp (p, p, n); /* exact */
  62   MPFR_ASSERTD (r > n);
  63   r -= n;
  64   /* since |p/2^r| < 1, and p is a non-zero integer, necessarily r > 0 */
  65
  66   MPFR_ASSERTD (mpz_sgn (p) > 0);
  67   MPFR_ASSERTD (m > 0);
  68
  69   /* check if p=1 (special case) */
  70   l = 0;
  71   /*
  72     We compute by binary splitting, with X = x^2 = p/2^r:
  73     P(a,b) = p if a+1=b, P(a,c)*P(c,b) otherwise
  74     Q(a,b) = (2a+1)*2^r if a+1=b [except Q(0,1)=1], Q(a,c)*Q(c,b) otherwise
  75     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
  76     Then atan(x)/x ~ S(0,i)/Q(0,i) for i so that (p/2^r)^i/i is small enough.
  77     The factor 2^(r*(b-a)) in Q(a,b) is implicit, thus we have to take it
  78     into account when we compute with Q.
  79   */
  80   accu[0] = 0; /* accu[k] = Mult[0] + ... + Mult[k], where Mult[j] is the
  81                   number of bits of the corresponding term S[j]/Q[j] */
  82   if (mpz_cmp_ui (p, 1) != 0)
  83     {
  84       /* p <> 1: precompute ptoj table */
  85       mpz_set (ptoj[0], p);
  86       for (im = 1 ; im <= m ; im ++)
  87         mpz_mul (ptoj[im], ptoj[im - 1], ptoj[im - 1]);
  88       /* main loop */
  89       n = 1UL << m;
  90       /* the ith term being X^i/(2i+1) with X=p/2^r, we can stop when
  91          p^i/2^(r*i) < 2^(-precy), i.e. r*i > precy + log2(p^i) */
  92       for (i = k = done = 0; (i < n) && (done == 0); i += 2, k ++)
  93         {
  94           /* initialize both S[k],Q[k] and S[k+1],Q[k+1] */
  95           mpz_set_ui (Q[k+1], 2 * i + 3); /* Q(i+1,i+2) */
  96           mpz_mul_ui (S[k+1], p, 2 * i + 1); /* S(i+1,i+2) */
  97           mpz_mul_2exp (S[k], Q[k+1], r);
  98           mpz_sub (S[k], S[k], S[k+1]); /* S(i,i+2) */
  99           mpz_mul_ui (Q[k], Q[k+1], 2 * i + 1); /* Q(i,i+2) */
 100           log2_nb_terms[k] = 1; /* S[k]/Q[k] corresponds to 2 terms */
 101           for (j = (i + 2) >> 1, l = 1; (j & 1) == 0; l ++, j >>= 1, k --)
 102             {
 103               /* invariant: S[k-1]/Q[k-1] and S[k]/Q[k] correspond
 104                  to 2^l terms each. We combine them into S[k-1]/Q[k-1] */
 105               MPFR_ASSERTD (k > 0);
 106               mpz_mul (S[k], S[k], Q[k-1]);
 107               mpz_mul (S[k], S[k], ptoj[l]);
 108               mpz_mul (S[k-1], S[k-1], Q[k]);
 109               mpz_mul_2exp (S[k-1], S[k-1], r << l);
 110               mpz_add (S[k-1], S[k-1], S[k]);
 111               mpz_mul (Q[k-1], Q[k-1], Q[k]);
 112               log2_nb_terms[k-1] = l + 1;
 113               /* now S[k-1]/Q[k-1] corresponds to 2^(l+1) terms */
 114               MPFR_MPZ_SIZEINBASE2(mult, ptoj[l+1]);
 115               /* FIXME: precompute bits(ptoj[l+1]) outside the loop? */
 116               mult = (r << (l + 1)) - mult - 1;
 117               accu[k-1] = (k == 1) ? mult : accu[k-2] + mult;
 118               if (accu[k-1] > precy)
 119                 done = 1;
 120             }
 121         }
 122     }
 123   else /* special case p=1: the ith term being X^i/(2i+1) with X=1/2^r,
 124           we can stop when r*i > precy i.e. i > precy/r */
 125     {
 126       n = 1UL << m;
 127       for (i = k = 0; (i < n) && (i <= precy / r); i += 2, k ++)
 128         {
 129           mpz_set_ui (Q[k + 1], 2 * i + 3);
 130           mpz_mul_2exp (S[k], Q[k+1], r);
 131           mpz_sub_ui (S[k], S[k], 1 + 2 * i);
 132           mpz_mul_ui (Q[k], Q[k + 1], 1 + 2 * i);
 133           log2_nb_terms[k] = 1; /* S[k]/Q[k] corresponds to 2 terms */
 134           for (j = (i + 2) >> 1, l = 1; (j & 1) == 0; l++, j >>= 1, k --)
 135             {
 136               MPFR_ASSERTD (k > 0);
 137               mpz_mul (S[k], S[k], Q[k-1]);
 138               mpz_mul (S[k-1], S[k-1], Q[k]);
 139               mpz_mul_2exp (S[k-1], S[k-1], r << l);
 140               mpz_add (S[k-1], S[k-1], S[k]);
 141               mpz_mul (Q[k-1], Q[k-1], Q[k]);
 142               log2_nb_terms[k-1] = l + 1;
 143             }
 144         }
 145     }
 146
 147   /* we need to combine S[0]/Q[0]...S[k-1]/Q[k-1] */
 148   l = 0; /* number of terms accumulated in S[k]/Q[k] */
 149   while (k > 1)
 150     {
 151       k --;
 152       /* combine S[k-1]/Q[k-1] and S[k]/Q[k] */
 153       j = log2_nb_terms[k-1];
 154       mpz_mul (S[k], S[k], Q[k-1]);
 155       if (mpz_cmp_ui (p, 1) != 0)
 156         mpz_mul (S[k], S[k], ptoj[j]);
 157       mpz_mul (S[k-1], S[k-1], Q[k]);
 158       l += 1 << log2_nb_terms[k];
 159       mpz_mul_2exp (S[k-1], S[k-1], r * l);
 160       mpz_add (S[k-1], S[k-1], S[k]);
 161       mpz_mul (Q[k-1], Q[k-1], Q[k]);
 162     }
 163   (*__gmp_free_func) (accu, (2 * m + 2) * sizeof (mp_prec_t));
 164
 165   MPFR_MPZ_SIZEINBASE2 (diff, S[0]);
 166   diff -= 2 * precy;
 167   expo = diff;
 168   if (diff >= 0)
 169     mpz_tdiv_q_2exp (S[0], S[0], diff);
 170   else
 171     mpz_mul_2exp (S[0], S[0], -diff);
 172
 173   MPFR_MPZ_SIZEINBASE2 (diff, Q[0]);
 174   diff -= precy;
 175   expo -= diff;
 176   if (diff >= 0)
 177     mpz_tdiv_q_2exp (Q[0], Q[0], diff);
 178   else
 179     mpz_mul_2exp (Q[0], Q[0], -diff);
 180
 181   mpz_tdiv_q (S[0], S[0], Q[0]);
 182   mpfr_set_z (y, S[0], GMP_RNDD);
 183   MPFR_SET_EXP (y, MPFR_EXP(y) + expo - r * (i - 1));
 184 }
 185
 186 int
 187 mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mp_rnd_t rnd_mode)
 188 {
 189   mpfr_t xp, arctgt, sk, tmp, tmp2;
 190   mpz_t  ukz;
 191   mpz_t *tabz;
 192   mp_exp_t exptol;
 193   mp_prec_t prec, realprec;
 194   unsigned long twopoweri;
 195   int comparaison, inexact;
 196   int i, n0, oldn0;
 197   MPFR_GROUP_DECL (group);
 198   MPFR_SAVE_EXPO_DECL (expo);
 199   MPFR_ZIV_DECL (loop);
 200
 201   MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd_mode),
 202                  ("atan[%#R]=%R inexact=%d", atan, atan, inexact));
 203
 204   /* Singular cases */
 205   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 206     {
 207       if (MPFR_IS_NAN (x))
 208         {
 209           MPFR_SET_NAN (atan);
 210           MPFR_RET_NAN;
 211         }
 212       else if (MPFR_IS_INF (x))
 213         {
 214           MPFR_SAVE_EXPO_MARK (expo);
 215           if (MPFR_IS_POS (x))  /* arctan(+inf) = Pi/2 */
 216             inexact = mpfr_const_pi (atan, rnd_mode);
 217           else /* arctan(-inf) = -Pi/2 */
 218             {
 219               inexact = -mpfr_const_pi (atan,
 220                                         MPFR_INVERT_RND (rnd_mode));
 221               MPFR_CHANGE_SIGN (atan);
 222             }
 223           mpfr_div_2ui (atan, atan, 1, rnd_mode);  /* exact (no exceptions) */
 224           MPFR_SAVE_EXPO_FREE (expo);
 225           return mpfr_check_range (atan, inexact, rnd_mode);
 226         }
 227       else /* x is necessarily 0 */
 228         {
 229           MPFR_ASSERTD (MPFR_IS_ZERO (x));
 230           MPFR_SET_ZERO (atan);
 231           MPFR_SET_SAME_SIGN (atan, x);
 232           MPFR_RET (0);
 233         }
 234     }
 235
 236   /* atan(x) = x - x^3/3 + x^5/5...
 237      so the error is < 2^(3*EXP(x)-1)
 238      so `EXP(x)-(3*EXP(x)-1)` = -2*EXP(x)+1 */
 239   MPFR_FAST_COMPUTE_IF_SMALL_INPUT (atan, x, -2 * MPFR_GET_EXP (x), 1, 0,
 240                                     rnd_mode, {});
 241
 242   /* Set x_p=|x| */
 243   MPFR_TMP_INIT_ABS (xp, x);
 244
 245   MPFR_SAVE_EXPO_MARK (expo);
 246
 247   /* Other simple case arctan(-+1)=-+pi/4 */
 248   comparaison = mpfr_cmp_ui (xp, 1);
 249   if (MPFR_UNLIKELY (comparaison == 0))
 250     {
 251       int neg = MPFR_IS_NEG (x);
 252       inexact = mpfr_const_pi (atan, MPFR_IS_POS (x) ? rnd_mode
 253                                : MPFR_INVERT_RND (rnd_mode));
 254       if (neg)
 255         {
 256           inexact = -inexact;
 257           MPFR_CHANGE_SIGN (atan);
 258         }
 259       mpfr_div_2ui (atan, atan, 2, rnd_mode);  /* exact (no exceptions) */
 260       MPFR_SAVE_EXPO_FREE (expo);
 261       return mpfr_check_range (atan, inexact, rnd_mode);
 262     }
 263
 264   realprec = MPFR_PREC (atan) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (atan)) + 4;
 265   prec = realprec + BITS_PER_MP_LIMB;
 266
 267   /* Initialisation    */
 268   mpz_init (ukz);
 269   MPFR_GROUP_INIT_4 (group, prec, sk, tmp, tmp2, arctgt);
 270   oldn0 = 0;
 271   tabz = (mpz_t *) 0;
 272
 273   MPFR_ZIV_INIT (loop, prec);
 274   for (;;)
 275     {
 276       /* First, if |x| < 1, we need to have more prec to be able to round (sup)
 277          n0 = ceil(log(prec_requested + 2 + 1+ln(2.4)/ln(2))/log(2)) */
 278       mp_prec_t sup;
 279 #if 0
 280       sup = 1;
 281       if (MPFR_GET_EXP (xp) < 0
 282           && (mpfr_uexp_t) (2-MPFR_GET_EXP (xp)) > realprec)
 283         sup = (mpfr_uexp_t) (2-MPFR_GET_EXP (xp)) - realprec;
 284 #else
 285       sup = MPFR_GET_EXP (xp) < 0 ? 2-MPFR_GET_EXP (xp) : 1;
 286 #endif
 287       n0 = MPFR_INT_CEIL_LOG2 ((realprec + sup) + 3);
 288       MPFR_ASSERTD (3*n0 > 2);
 289       prec = (realprec + sup) + 1 + MPFR_INT_CEIL_LOG2 (3*n0-2);
 290
 291       /* Initialisation */
 292       MPFR_GROUP_REPREC_4 (group, prec, sk, tmp, tmp2, arctgt);
 293       if (MPFR_LIKELY (oldn0 == 0))
 294         {
 295           oldn0 = 3*(n0+1);
 296           tabz = (mpz_t *) (*__gmp_allocate_func) (oldn0*sizeof (mpz_t));
 297           for (i = 0; i < oldn0; i++)
 298             mpz_init (tabz[i]);
 299         }
 300       else if (MPFR_UNLIKELY (oldn0 < 3 * (n0 + 1)))
 301         {
 302           tabz = (mpz_t *) (*__gmp_reallocate_func)
 303             (tabz, oldn0*sizeof (mpz_t), 3*(n0+1)*sizeof (mpz_t));
 304           for (i = oldn0; i < 3*(n0+1); i++)
 305             mpz_init (tabz[i]);
 306           oldn0 = 3*(n0+1);
 307         }
 308
 309       if (comparaison > 0)
 310         mpfr_ui_div (sk, 1, xp, GMP_RNDN);
 311       else
 312         mpfr_set (sk, xp, GMP_RNDN);
 313
 314       /* sk is 1/|x| if |x| > 1, and |x| otherwise, i.e. min(|x|, 1/|x|) */
 315
 316       /* If sk=1, then if |x| < 1, we have 1 - 2^(-prec-1) <= |x| < 1,
 317          or if |x| > 1, we have 1 - 2^(-prec-1) <= 1/|x| < 1, thus in all
 318          cases ||x| - 1| <= 2^(-prec), from which it follows
 319          |atan|x| - Pi/4| <= 2^(-prec), given the Taylor expansion
 320          atan(1+x) = Pi/4 + x/2 - x^2/4 + ...
 321          Since Pi/4 = 0.785..., the error is at most one ulp.
 322       */
 323       if (MPFR_UNLIKELY(mpfr_cmp_ui (sk, 1) == 0))
 324         {
 325           mpfr_const_pi (arctgt, GMP_RNDN); /* 1/2 ulp extra error */
 326           mpfr_div_2ui (arctgt, arctgt, 2, GMP_RNDN); /* exact */
 327           realprec = prec - 2;
 328           goto can_round;
 329         }
 330
 331       /* Assignation  */
 332       MPFR_SET_ZERO (arctgt);
 333       twopoweri = 1 << 0;
 334       MPFR_ASSERTD (n0 >= 4);
 335       /* FIXME: further reduce the argument so that it is less than
 336          1/n where n is the output precision. In such a way, the
 337          first calls to mpfr_atan_aux will not be too expensive,
 338          since the number of needed terms will be n/log(n), so the
 339          factorial contribution will be O(n). */
 340       for (i = 0 ; i < n0; i++)
 341         {
 342           if (MPFR_UNLIKELY (MPFR_IS_ZERO (sk)))
 343             break;
 344           /* Calculation of trunc(tmp) --> mpz */
 345           mpfr_mul_2ui (tmp, sk, twopoweri, GMP_RNDN);
 346           mpfr_trunc (tmp, tmp);
 347           if (!MPFR_IS_ZERO (tmp))
 348             {
 349               /* tmp = ukz*2^exptol */
 350               exptol = mpfr_get_z_exp (ukz, tmp);
 351               /* since the s_k are decreasing (see algorithms.tex),
 352                  and s_0 = min(|x|, 1/|x|) < 1, we have sk < 1,
 353                  thus exptol < 0 */
 354               MPFR_ASSERTD (exptol < 0);
 355               mpz_tdiv_q_2exp (ukz, ukz, (unsigned long int) (-exptol));
 356               /* since tmp is a non-zero integer, and tmp = ukzold*2^exptol,
 357                  we now have ukz = tmp, thus ukz is non-zero */
 358               /* Calculation of arctan(Ak) */
 359               mpfr_set_z (tmp, ukz, GMP_RNDN);
 360               mpfr_div_2ui (tmp, tmp, twopoweri, GMP_RNDN);
 361               mpfr_atan_aux (tmp2, ukz, twopoweri, n0 - i, tabz);
 362               mpfr_mul (tmp2, tmp2, tmp, GMP_RNDN);
 363               /* Addition */
 364               mpfr_add (arctgt, arctgt, tmp2, GMP_RNDN);
 365               /* Next iteration */
 366               mpfr_sub (tmp2, sk, tmp, GMP_RNDN);
 367               mpfr_mul (sk, sk, tmp, GMP_RNDN);
 368               mpfr_add_ui (sk, sk, 1, GMP_RNDN);
 369               mpfr_div (sk, tmp2, sk, GMP_RNDN);
 370             }
 371           twopoweri <<= 1;
 372         }
 373       /* Add last step (Arctan(sk) ~= sk */
 374       mpfr_add (arctgt, arctgt, sk, GMP_RNDN);
 375       if (comparaison > 0)
 376         {
 377           mpfr_const_pi (tmp, GMP_RNDN);
 378           mpfr_div_2ui (tmp, tmp, 1, GMP_RNDN);
 379           mpfr_sub (arctgt, tmp, arctgt, GMP_RNDN);
 380         }
 381       MPFR_SET_POS (arctgt);
 382
 383     can_round:
 384       if (MPFR_LIKELY (MPFR_CAN_ROUND (arctgt, realprec, MPFR_PREC (atan),
 385                                        rnd_mode)))
 386         break;
 387       MPFR_ZIV_NEXT (loop, realprec);
 388     }
 389   MPFR_ZIV_FREE (loop);
 390
 391   inexact = mpfr_set4 (atan, arctgt, rnd_mode, MPFR_SIGN (x));
 392
 393   for (i = 0 ; i < oldn0 ; i++)
 394     mpz_clear (tabz[i]);
 395   mpz_clear (ukz);
 396   (*__gmp_free_func) (tabz, oldn0*sizeof (mpz_t));
 397   MPFR_GROUP_CLEAR (group);
 398
 399   MPFR_SAVE_EXPO_FREE (expo);
 400   return mpfr_check_range (arctgt, inexact, rnd_mode);
 401 }