contrib/mpfr/cos.c

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