contrib/mpfr/src/cos.c

   1 /* mpfr_cos -- cosine of a floating-point number
   2
   3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Caramel projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   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 3 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.LESSER.  If not, see
  20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 static int
  27 mpfr_cos_fast (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
  28 {
  29   int inex;
  30
  31   inex = mpfr_sincos_fast (NULL, y, x, rnd_mode);
  32   inex = inex >> 2; /* 0: exact, 1: rounded up, 2: rounded down */
  33   return (inex == 2) ? -1 : inex;
  34 }
  35
  36 /* f <- 1 - r/2! + r^2/4! + ... + (-1)^l r^l/(2l)! + ...
  37    Assumes |r| < 1/2, and f, r have the same precision.
  38    Returns e such that the error on f is bounded by 2^e ulps.
  39 */
  40 static int
  41 mpfr_cos2_aux (mpfr_ptr f, mpfr_srcptr r)
  42 {
  43   mpz_t x, t, s;
  44   mpfr_exp_t ex, l, m;
  45   mpfr_prec_t p, q;
  46   unsigned long i, maxi, imax;
  47
  48   MPFR_ASSERTD(mpfr_get_exp (r) <= -1);
  49
  50   /* compute minimal i such that i*(i+1) does not fit in an unsigned long,
  51      assuming that there are no padding bits. */
  52   maxi = 1UL << (CHAR_BIT * sizeof(unsigned long) / 2);
  53   if (maxi * (maxi / 2) == 0) /* test checked at compile time */
  54     {
  55       /* can occur only when there are padding bits. */
  56       /* maxi * (maxi-1) is representable iff maxi * (maxi / 2) != 0 */
  57       do
  58         maxi /= 2;
  59       while (maxi * (maxi / 2) == 0);
  60     }
  61
  62   mpz_init (x);
  63   mpz_init (s);
  64   mpz_init (t);
  65   ex = mpfr_get_z_2exp (x, r); /* r = x*2^ex */
  66
  67   /* remove trailing zeroes */
  68   l = mpz_scan1 (x, 0);
  69   ex += l;
  70   mpz_fdiv_q_2exp (x, x, l);
  71
  72   /* since |r| < 1, r = x*2^ex, and x is an integer, necessarily ex < 0 */
  73
  74   p = mpfr_get_prec (f); /* same than r */
  75   /* bound for number of iterations */
  76   imax = p / (-mpfr_get_exp (r));
  77   imax += (imax == 0);
  78   q = 2 * MPFR_INT_CEIL_LOG2(imax) + 4; /* bound for (3l)^2 */
  79
  80   mpz_set_ui (s, 1); /* initialize sum with 1 */
  81   mpz_mul_2exp (s, s, p + q); /* scale all values by 2^(p+q) */
  82   mpz_set (t, s); /* invariant: t is previous term */
  83   for (i = 1; (m = mpz_sizeinbase (t, 2)) >= q; i += 2)
  84     {
  85       /* adjust precision of x to that of t */
  86       l = mpz_sizeinbase (x, 2);
  87       if (l > m)
  88         {
  89           l -= m;
  90           mpz_fdiv_q_2exp (x, x, l);
  91           ex += l;
  92         }
  93       /* multiply t by r */
  94       mpz_mul (t, t, x);
  95       mpz_fdiv_q_2exp (t, t, -ex);
  96       /* divide t by i*(i+1) */
  97       if (i < maxi)
  98         mpz_fdiv_q_ui (t, t, i * (i + 1));
  99       else
 100         {
 101           mpz_fdiv_q_ui (t, t, i);
 102           mpz_fdiv_q_ui (t, t, i + 1);
 103         }
 104       /* if m is the (current) number of bits of t, we can consider that
 105          all operations on t so far had precision >= m, so we can prove
 106          by induction that the relative error on t is of the form
 107          (1+u)^(3l)-1, where |u| <= 2^(-m), and l=(i+1)/2 is the # of loops.
 108          Since |(1+x^2)^(1/x) - 1| <= 4x/3 for |x| <= 1/2,
 109          for |u| <= 1/(3l)^2, the absolute error is bounded by
 110          4/3*(3l)*2^(-m)*t <= 4*l since |t| < 2^m.
 111          Therefore the error on s is bounded by 2*l*(l+1). */
 112       /* add or subtract to s */
 113       if (i % 4 == 1)
 114         mpz_sub (s, s, t);
 115       else
 116         mpz_add (s, s, t);
 117     }
 118
 119   mpfr_set_z (f, s, MPFR_RNDN);
 120   mpfr_div_2ui (f, f, p + q, MPFR_RNDN);
 121
 122   mpz_clear (x);
 123   mpz_clear (s);
 124   mpz_clear (t);
 125
 126   l = (i - 1) / 2; /* number of iterations */
 127   return 2 * MPFR_INT_CEIL_LOG2 (l + 1) + 1; /* bound is 2l(l+1) */
 128 }
 129
 130 int
 131 mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
 132 {
 133   mpfr_prec_t K0, K, precy, m, k, l;
 134   int inexact, reduce = 0;
 135   mpfr_t r, s, xr, c;
 136   mpfr_exp_t exps, cancel = 0, expx;
 137   MPFR_ZIV_DECL (loop);
 138   MPFR_SAVE_EXPO_DECL (expo);
 139   MPFR_GROUP_DECL (group);
 140
 141   MPFR_LOG_FUNC (
 142     ("x[%Pu]=%*.Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
 143     ("y[%Pu]=%*.Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y,
 144      inexact));
 145
 146   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 147     {
 148       if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
 149         {
 150           MPFR_SET_NAN (y);
 151           MPFR_RET_NAN;
 152         }
 153       else
 154         {
 155           MPFR_ASSERTD (MPFR_IS_ZERO (x));
 156           return mpfr_set_ui (y, 1, rnd_mode);
 157         }
 158     }
 159
 160   MPFR_SAVE_EXPO_MARK (expo);
 161
 162   /* cos(x) = 1-x^2/2 + ..., so error < 2^(2*EXP(x)-1) */
 163   expx = MPFR_GET_EXP (x);
 164   MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, __gmpfr_one, -2 * expx,
 165                                     1, 0, rnd_mode, expo, {});
 166
 167   /* Compute initial precision */
 168   precy = MPFR_PREC (y);
 169
 170   if (precy >= MPFR_SINCOS_THRESHOLD)
 171     {
 172       MPFR_SAVE_EXPO_FREE (expo);
 173       return mpfr_cos_fast (y, x, rnd_mode);
 174     }
 175
 176   K0 = __gmpfr_isqrt (precy / 3);
 177   m = precy + 2 * MPFR_INT_CEIL_LOG2 (precy) + 2 * K0;
 178
 179   if (expx >= 3)
 180     {
 181       reduce = 1;
 182       /* As expx + m - 1 will silently be converted into mpfr_prec_t
 183          in the mpfr_init2 call, the assert below may be useful to
 184          avoid undefined behavior. */
 185       MPFR_ASSERTN (expx + m - 1 <= MPFR_PREC_MAX);
 186       mpfr_init2 (c, expx + m - 1);
 187       mpfr_init2 (xr, m);
 188     }
 189
 190   MPFR_GROUP_INIT_2 (group, m, r, s);
 191   MPFR_ZIV_INIT (loop, m);
 192   for (;;)
 193     {
 194       /* If |x| >= 4, first reduce x cmod (2*Pi) into xr, using mpfr_remainder:
 195          let e = EXP(x) >= 3, and m the target precision:
 196          (1) c <- 2*Pi              [precision e+m-1, nearest]
 197          (2) xr <- remainder (x, c) [precision m, nearest]
 198          We have |c - 2*Pi| <= 1/2ulp(c) = 2^(3-e-m)
 199                  |xr - x - k c| <= 1/2ulp(xr) <= 2^(1-m)
 200                  |k| <= |x|/(2*Pi) <= 2^(e-2)
 201          Thus |xr - x - 2kPi| <= |k| |c - 2Pi| + 2^(1-m) <= 2^(2-m).
 202          It follows |cos(xr) - cos(x)| <= 2^(2-m). */
 203       if (reduce)
 204         {
 205           mpfr_const_pi (c, MPFR_RNDN);
 206           mpfr_mul_2ui (c, c, 1, MPFR_RNDN); /* 2Pi */
 207           mpfr_remainder (xr, x, c, MPFR_RNDN);
 208           if (MPFR_IS_ZERO(xr))
 209             goto ziv_next;
 210           /* now |xr| <= 4, thus r <= 16 below */
 211           mpfr_mul (r, xr, xr, MPFR_RNDU); /* err <= 1 ulp */
 212         }
 213       else
 214         mpfr_mul (r, x, x, MPFR_RNDU); /* err <= 1 ulp */
 215
 216       /* now |x| < 4 (or xr if reduce = 1), thus |r| <= 16 */
 217
 218       /* we need |r| < 1/2 for mpfr_cos2_aux, i.e., EXP(r) - 2K <= -1 */
 219       K = K0 + 1 + MAX(0, MPFR_GET_EXP(r)) / 2;
 220       /* since K0 >= 0, if EXP(r) < 0, then K >= 1, thus EXP(r) - 2K <= -3;
 221          otherwise if EXP(r) >= 0, then K >= 1/2 + EXP(r)/2, thus
 222          EXP(r) - 2K <= -1 */
 223
 224       MPFR_SET_EXP (r, MPFR_GET_EXP (r) - 2 * K); /* Can't overflow! */
 225
 226       /* s <- 1 - r/2! + ... + (-1)^l r^l/(2l)! */
 227       l = mpfr_cos2_aux (s, r);
 228       /* l is the error bound in ulps on s */
 229       MPFR_SET_ONE (r);
 230       for (k = 0; k < K; k++)
 231         {
 232           mpfr_sqr (s, s, MPFR_RNDU);            /* err <= 2*olderr */
 233           MPFR_SET_EXP (s, MPFR_GET_EXP (s) + 1); /* Can't overflow */
 234           mpfr_sub (s, s, r, MPFR_RNDN);         /* err <= 4*olderr */
 235           if (MPFR_IS_ZERO(s))
 236             goto ziv_next;
 237           MPFR_ASSERTD (MPFR_GET_EXP (s) <= 1);
 238         }
 239
 240       /* The absolute error on s is bounded by (2l+1/3)*2^(2K-m)
 241          2l+1/3 <= 2l+1.
 242          If |x| >= 4, we need to add 2^(2-m) for the argument reduction
 243          by 2Pi: if K = 0, this amounts to add 4 to 2l+1/3, i.e., to add
 244          2 to l; if K >= 1, this amounts to add 1 to 2*l+1/3. */
 245       l = 2 * l + 1;
 246       if (reduce)
 247         l += (K == 0) ? 4 : 1;
 248       k = MPFR_INT_CEIL_LOG2 (l) + 2*K;
 249       /* now the error is bounded by 2^(k-m) = 2^(EXP(s)-err) */
 250
 251       exps = MPFR_GET_EXP (s);
 252       if (MPFR_LIKELY (MPFR_CAN_ROUND (s, exps + m - k, precy, rnd_mode)))
 253         break;
 254
 255       if (MPFR_UNLIKELY (exps == 1))
 256         /* s = 1 or -1, and except x=0 which was already checked above,
 257            cos(x) cannot be 1 or -1, so we can round if the error is less
 258            than 2^(-precy) for directed rounding, or 2^(-precy-1) for rounding
 259            to nearest. */
 260         {
 261           if (m > k && (m - k >= precy + (rnd_mode == MPFR_RNDN)))
 262             {
 263               /* If round to nearest or away, result is s = 1 or -1,
 264                  otherwise it is round(nexttoward (s, 0)). However in order to
 265                  have the inexact flag correctly set below, we set |s| to
 266                  1 - 2^(-m) in all cases. */
 267               mpfr_nexttozero (s);
 268               break;
 269             }
 270         }
 271
 272       if (exps < cancel)
 273         {
 274           m += cancel - exps;
 275           cancel = exps;
 276         }
 277
 278     ziv_next:
 279       MPFR_ZIV_NEXT (loop, m);
 280       MPFR_GROUP_REPREC_2 (group, m, r, s);
 281       if (reduce)
 282         {
 283           mpfr_set_prec (xr, m);
 284           mpfr_set_prec (c, expx + m - 1);
 285         }
 286     }
 287   MPFR_ZIV_FREE (loop);
 288   inexact = mpfr_set (y, s, rnd_mode);
 289   MPFR_GROUP_CLEAR (group);
 290   if (reduce)
 291     {
 292       mpfr_clear (xr);
 293       mpfr_clear (c);
 294     }
 295
 296   MPFR_SAVE_EXPO_FREE (expo);
 297   return mpfr_check_range (y, inexact, rnd_mode);
 298 }