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       /* As expx + m - 1 will silently be converted into mpfr_prec_t
 164          in the mpfr_init2 call, the assert below may be useful to
 165          avoid undefined behavior. */
 166       MPFR_ASSERTN (expx + m - 1 <= MPFR_PREC_MAX);
 167       mpfr_init2 (c, expx + m - 1);
 168       mpfr_init2 (xr, m);
 169     }
 170
 171   MPFR_GROUP_INIT_2 (group, m, r, s);
 172   MPFR_ZIV_INIT (loop, m);
 173   for (;;)
 174     {
 175       /* If |x| >= 4, first reduce x cmod (2*Pi) into xr, using mpfr_remainder:
 176          let e = EXP(x) >= 3, and m the target precision:
 177          (1) c <- 2*Pi              [precision e+m-1, nearest]
 178          (2) xr <- remainder (x, c) [precision m, nearest]
 179          We have |c - 2*Pi| <= 1/2ulp(c) = 2^(3-e-m)
 180                  |xr - x - k c| <= 1/2ulp(xr) <= 2^(1-m)
 181                  |k| <= |x|/(2*Pi) <= 2^(e-2)
 182          Thus |xr - x - 2kPi| <= |k| |c - 2Pi| + 2^(1-m) <= 2^(2-m).
 183          It follows |cos(xr) - cos(x)| <= 2^(2-m). */
 184       if (reduce)
 185         {
 186           mpfr_const_pi (c, GMP_RNDN);
 187           mpfr_mul_2ui (c, c, 1, GMP_RNDN); /* 2Pi */
 188           mpfr_remainder (xr, x, c, GMP_RNDN);
 189           if (MPFR_IS_ZERO(xr))
 190             goto ziv_next;
 191           /* now |xr| <= 4, thus r <= 16 below */
 192           mpfr_mul (r, xr, xr, GMP_RNDU); /* err <= 1 ulp */
 193         }
 194       else
 195         mpfr_mul (r, x, x, GMP_RNDU); /* err <= 1 ulp */
 196
 197       /* now |x| < 4 (or xr if reduce = 1), thus |r| <= 16 */
 198
 199       /* we need |r| < 1/2 for mpfr_cos2_aux, i.e., EXP(r) - 2K <= -1 */
 200       K = K0 + 1 + MAX(0, MPFR_EXP(r)) / 2;
 201       /* since K0 >= 0, if EXP(r) < 0, then K >= 1, thus EXP(r) - 2K <= -3;
 202          otherwise if EXP(r) >= 0, then K >= 1/2 + EXP(r)/2, thus
 203          EXP(r) - 2K <= -1 */
 204
 205       MPFR_SET_EXP (r, MPFR_GET_EXP (r) - 2 * K); /* Can't overflow! */
 206
 207       /* s <- 1 - r/2! + ... + (-1)^l r^l/(2l)! */
 208       l = mpfr_cos2_aux (s, r);
 209       /* l is the error bound in ulps on s */
 210       MPFR_SET_ONE (r);
 211       for (k = 0; k < K; k++)
 212         {
 213           mpfr_sqr (s, s, GMP_RNDU);            /* err <= 2*olderr */
 214           MPFR_SET_EXP (s, MPFR_GET_EXP (s) + 1); /* Can't overflow */
 215           mpfr_sub (s, s, r, GMP_RNDN);         /* err <= 4*olderr */
 216           if (MPFR_IS_ZERO(s))
 217             goto ziv_next;
 218           MPFR_ASSERTD (MPFR_GET_EXP (s) <= 1);
 219         }
 220
 221       /* The absolute error on s is bounded by (2l+1/3)*2^(2K-m)
 222          2l+1/3 <= 2l+1.
 223          If |x| >= 4, we need to add 2^(2-m) for the argument reduction
 224          by 2Pi: if K = 0, this amounts to add 4 to 2l+1/3, i.e., to add
 225          2 to l; if K >= 1, this amounts to add 1 to 2*l+1/3. */
 226       l = 2 * l + 1;
 227       if (reduce)
 228         l += (K == 0) ? 4 : 1;
 229       k = MPFR_INT_CEIL_LOG2 (l) + 2*K;
 230       /* now the error is bounded by 2^(k-m) = 2^(EXP(s)-err) */
 231
 232       exps = MPFR_GET_EXP (s);
 233       if (MPFR_LIKELY (MPFR_CAN_ROUND (s, exps + m - k, precy, rnd_mode)))
 234         break;
 235
 236       if (MPFR_UNLIKELY (exps == 1))
 237         /* s = 1 or -1, and except x=0 which was already checked above,
 238            cos(x) cannot be 1 or -1, so we can round if the error is less
 239            than 2^(-precy) for directed rounding, or 2^(-precy-1) for rounding
 240            to nearest. */
 241         {
 242           if (m > k && (m - k >= precy + (rnd_mode == GMP_RNDN)))
 243             {
 244               /* If round to nearest or away, result is s = 1 or -1,
 245                  otherwise it is round(nexttoward (s, 0)). However in order to
 246                  have the inexact flag correctly set below, we set |s| to
 247                  1 - 2^(-m) in all cases. */
 248               mpfr_nexttozero (s);
 249               break;
 250             }
 251         }
 252
 253       if (exps < cancel)
 254         {
 255           m += cancel - exps;
 256           cancel = exps;
 257         }
 258
 259     ziv_next:
 260       MPFR_ZIV_NEXT (loop, m);
 261       MPFR_GROUP_REPREC_2 (group, m, r, s);
 262       if (reduce)
 263         {
 264           mpfr_set_prec (xr, m);
 265           mpfr_set_prec (c, expx + m - 1);
 266         }
 267     }
 268   MPFR_ZIV_FREE (loop);
 269   inexact = mpfr_set (y, s, rnd_mode);
 270   MPFR_GROUP_CLEAR (group);
 271   if (reduce)
 272     {
 273       mpfr_clear (xr);
 274       mpfr_clear (c);
 275     }
 276
 277   MPFR_SAVE_EXPO_FREE (expo);
 278   return mpfr_check_range (y, inexact, rnd_mode);
 279 }