contrib/mpfr/sin_cos.c

   1 /* mpfr_sin_cos -- sine and cosine of a floating-point number
   2
   3 Copyright 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 /* (y, z) <- (sin(x), cos(x)), return value is 0 iff both results are exact
  27    ie, iff x = 0 */
  28 int
  29 mpfr_sin_cos (mpfr_ptr y, mpfr_ptr z, mpfr_srcptr x, mp_rnd_t rnd_mode)
  30 {
  31   mp_prec_t prec, m;
  32   int neg, reduce;
  33   mpfr_t c, xr;
  34   mpfr_srcptr xx;
  35   mp_exp_t err, expx;
  36   int inexy, inexz;
  37   MPFR_ZIV_DECL (loop);
  38   MPFR_SAVE_EXPO_DECL (expo);
  39
  40   MPFR_ASSERTN (y != z);
  41
  42   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
  43     {
  44       if (MPFR_IS_NAN(x) || MPFR_IS_INF(x))
  45         {
  46           MPFR_SET_NAN (y);
  47           MPFR_SET_NAN (z);
  48           MPFR_RET_NAN;
  49         }
  50       else /* x is zero */
  51         {
  52           MPFR_ASSERTD (MPFR_IS_ZERO (x));
  53           MPFR_SET_ZERO (y);
  54           MPFR_SET_SAME_SIGN (y, x);
  55           /* y = 0, thus exact, but z is inexact in case of underflow
  56              or overflow */
  57           return mpfr_set_ui (z, 1, rnd_mode);
  58         }
  59     }
  60
  61   MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd_mode),
  62                   ("sin[%#R]=%R cos[%#R]=%R", y, y, z, z));
  63
  64   MPFR_SAVE_EXPO_MARK (expo);
  65
  66   prec = MAX (MPFR_PREC (y), MPFR_PREC (z));
  67   m = prec + MPFR_INT_CEIL_LOG2 (prec) + 13;
  68   expx = MPFR_GET_EXP (x);
  69
  70   /* When x is close to 0, say 2^(-k), then there is a cancellation of about
  71      2k bits in 1-cos(x)^2. FIXME: in that case, it would be more efficient
  72      to compute sin(x) directly. VL: This is partly done by using
  73      MPFR_FAST_COMPUTE_IF_SMALL_INPUT from the mpfr_sin and mpfr_cos
  74      functions. Moreover, any overflow on m is avoided. */
  75   if (expx < 0)
  76     {
  77       /* Warning: in case y = x, and the first call to
  78          MPFR_FAST_COMPUTE_IF_SMALL_INPUT succeeds but the second fails,
  79          we will have clobbered the original value of x.
  80          The workaround is to first compute z = cos(x) in that case, since
  81          y and z are different. */
  82       if (y != x)
  83         /* y and x differ, thus we can safely try to compute y first */
  84         {
  85           MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2 * expx, 2, 0, rnd_mode,
  86                                             { inexy = _inexact;
  87                                               goto small_input; });
  88           if (0)
  89             {
  90             small_input:
  91               /* we can go here only if we can round sin(x) */
  92               MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, -2 * expx,
  93                                                 1, 0, rnd_mode,
  94                                                 { inexz = _inexact;
  95                                                   goto end; });
  96             }
  97
  98           /* if we go here, one of the two MPFR_FAST_COMPUTE_IF_SMALL_INPUT
  99              calls failed */
 100         }
 101       else /* y and x are the same variable: try to compute z first, which
 102               necessarily differs */
 103         {
 104           MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, -2 * expx,
 105                                             1, 0, rnd_mode,
 106                                             { inexz = _inexact;
 107                                               goto small_input2; });
 108           if (0)
 109             {
 110             small_input2:
 111               /* we can go here only if we can round cos(x) */
 112               MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2 * expx, 2, 0,
 113                                                 rnd_mode,
 114                                                 { inexy = _inexact;
 115                                                   goto end; });
 116             }
 117         }
 118       m += 2 * (-expx);
 119     }
 120
 121   mpfr_init (c);
 122   mpfr_init (xr);
 123
 124   MPFR_ZIV_INIT (loop, m);
 125   for (;;)
 126     {
 127       /* the following is copied from sin.c */
 128       if (expx >= 2) /* reduce the argument */
 129         {
 130           reduce = 1;
 131           mpfr_set_prec (c, expx + m - 1);
 132           mpfr_set_prec (xr, m);
 133           mpfr_const_pi (c, GMP_RNDN);
 134           mpfr_mul_2ui (c, c, 1, GMP_RNDN);
 135           mpfr_remainder (xr, x, c, GMP_RNDN);
 136           mpfr_div_2ui (c, c, 1, GMP_RNDN);
 137           if (MPFR_SIGN (xr) > 0)
 138             mpfr_sub (c, c, xr, GMP_RNDZ);
 139           else
 140             mpfr_add (c, c, xr, GMP_RNDZ);
 141           if (MPFR_IS_ZERO(xr) || MPFR_EXP(xr) < (mp_exp_t) 3 - (mp_exp_t) m
 142               || MPFR_EXP(c) < (mp_exp_t) 3 - (mp_exp_t) m)
 143             goto next_step;
 144           xx = xr;
 145         }
 146       else /* the input argument is already reduced */
 147         {
 148           reduce = 0;
 149           xx = x;
 150         }
 151
 152       neg = MPFR_IS_NEG (xx); /* gives sign of sin(x) */
 153       mpfr_set_prec (c, m);
 154       mpfr_cos (c, xx, GMP_RNDZ);
 155       /* If no argument reduction was performed, the error is at most ulp(c),
 156          otherwise it is at most ulp(c) + 2^(2-m). Since |c| < 1, we have
 157          ulp(c) <= 2^(-m), thus the error is bounded by 2^(3-m) in that later
 158          case. */
 159       if (reduce == 0)
 160         err = m;
 161       else
 162         err = MPFR_GET_EXP (c) + (mp_exp_t) (m - 3);
 163       if (!mpfr_can_round (c, err, GMP_RNDN, rnd_mode,
 164                            MPFR_PREC (z) + (rnd_mode == GMP_RNDN)))
 165         goto next_step;
 166
 167       /* we can't set z now, because in case z = x, and the mpfr_can_round()
 168          call below fails, we will have clobbered the input */
 169       mpfr_set_prec (xr, MPFR_PREC(c));
 170       mpfr_swap (xr, c); /* save the approximation of the cosine in xr */
 171       mpfr_sqr (c, xr, GMP_RNDU);
 172       mpfr_ui_sub (c, 1, c, GMP_RNDN);
 173       err = 2 + (- MPFR_GET_EXP (c)) / 2;
 174       mpfr_sqrt (c, c, GMP_RNDN);
 175       if (neg)
 176         MPFR_CHANGE_SIGN (c);
 177
 178       /* the absolute error on c is at most 2^(err-m), which we must put
 179          in the form 2^(EXP(c)-err). If there was an argument reduction,
 180          we need to add 2^(2-m); since err >= 2, the error is bounded by
 181          2^(err+1-m) in that case. */
 182       err = MPFR_GET_EXP (c) + (mp_exp_t) m - (err + reduce);
 183       if (mpfr_can_round (c, err, GMP_RNDN, rnd_mode,
 184                           MPFR_PREC (y) + (rnd_mode == GMP_RNDN)))
 185         break;
 186       /* check for huge cancellation */
 187       if (err < (mp_exp_t) MPFR_PREC (y))
 188         m += MPFR_PREC (y) - err;
 189       /* Check if near 1 */
 190       if (MPFR_GET_EXP (c) == 1
 191           && MPFR_MANT (c)[MPFR_LIMB_SIZE (c)-1] == MPFR_LIMB_HIGHBIT)
 192         m += m;
 193
 194     next_step:
 195       MPFR_ZIV_NEXT (loop, m);
 196       mpfr_set_prec (c, m);
 197     }
 198   MPFR_ZIV_FREE (loop);
 199
 200   inexy = mpfr_set (y, c, rnd_mode);
 201   inexz = mpfr_set (z, xr, rnd_mode);
 202
 203   mpfr_clear (c);
 204   mpfr_clear (xr);
 205
 206  end:
 207   /* FIXME: update the underflow flag if need be. */
 208   MPFR_SAVE_EXPO_FREE (expo);
 209   mpfr_check_range (y, inexy, rnd_mode);
 210   mpfr_check_range (z, inexz, rnd_mode);
 211   MPFR_RET (1); /* Always inexact */
 212 }