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 (
  86             y, x, -2 * expx, 2, 0, rnd_mode,
  87             { inexy = _inexact;
  88               goto small_input; });
  89           if (0)
  90             {
  91             small_input:
  92               /* we can go here only if we can round sin(x) */
  93               MPFR_FAST_COMPUTE_IF_SMALL_INPUT (
  94                 z, __gmpfr_one, -2 * expx, 1, 0, rnd_mode,
  95                 { inexz = _inexact;
  96                   MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
  97                   goto end; });
  98             }
  99
 100           /* if we go here, one of the two MPFR_FAST_COMPUTE_IF_SMALL_INPUT
 101              calls failed */
 102         }
 103       else /* y and x are the same variable: try to compute z first, which
 104               necessarily differs */
 105         {
 106           MPFR_FAST_COMPUTE_IF_SMALL_INPUT (
 107             z, __gmpfr_one, -2 * expx, 1, 0, rnd_mode,
 108             { inexz = _inexact;
 109               goto small_input2; });
 110           if (0)
 111             {
 112             small_input2:
 113               /* we can go here only if we can round cos(x) */
 114               MPFR_FAST_COMPUTE_IF_SMALL_INPUT (
 115                 y, x, -2 * expx, 2, 0, rnd_mode,
 116                 { inexy = _inexact;
 117                   MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
 118                   goto end; });
 119             }
 120         }
 121       m += 2 * (-expx);
 122     }
 123
 124   mpfr_init (c);
 125   mpfr_init (xr);
 126
 127   MPFR_ZIV_INIT (loop, m);
 128   for (;;)
 129     {
 130       /* the following is copied from sin.c */
 131       if (expx >= 2) /* reduce the argument */
 132         {
 133           reduce = 1;
 134           mpfr_set_prec (c, expx + m - 1);
 135           mpfr_set_prec (xr, m);
 136           mpfr_const_pi (c, GMP_RNDN);
 137           mpfr_mul_2ui (c, c, 1, GMP_RNDN);
 138           mpfr_remainder (xr, x, c, GMP_RNDN);
 139           mpfr_div_2ui (c, c, 1, GMP_RNDN);
 140           if (MPFR_SIGN (xr) > 0)
 141             mpfr_sub (c, c, xr, GMP_RNDZ);
 142           else
 143             mpfr_add (c, c, xr, GMP_RNDZ);
 144           if (MPFR_IS_ZERO(xr) || MPFR_EXP(xr) < (mp_exp_t) 3 - (mp_exp_t) m
 145               || MPFR_EXP(c) < (mp_exp_t) 3 - (mp_exp_t) m)
 146             goto next_step;
 147           xx = xr;
 148         }
 149       else /* the input argument is already reduced */
 150         {
 151           reduce = 0;
 152           xx = x;
 153         }
 154
 155       neg = MPFR_IS_NEG (xx); /* gives sign of sin(x) */
 156       mpfr_set_prec (c, m);
 157       mpfr_cos (c, xx, GMP_RNDZ);
 158       /* If no argument reduction was performed, the error is at most ulp(c),
 159          otherwise it is at most ulp(c) + 2^(2-m). Since |c| < 1, we have
 160          ulp(c) <= 2^(-m), thus the error is bounded by 2^(3-m) in that later
 161          case. */
 162       if (reduce == 0)
 163         err = m;
 164       else
 165         err = MPFR_GET_EXP (c) + (mp_exp_t) (m - 3);
 166       if (!mpfr_can_round (c, err, GMP_RNDN, GMP_RNDZ,
 167                            MPFR_PREC (z) + (rnd_mode == GMP_RNDN)))
 168         goto next_step;
 169
 170       /* we can't set z now, because in case z = x, and the mpfr_can_round()
 171          call below fails, we will have clobbered the input */
 172       mpfr_set_prec (xr, MPFR_PREC(c));
 173       mpfr_swap (xr, c); /* save the approximation of the cosine in xr */
 174       mpfr_sqr (c, xr, GMP_RNDU); /* the absolute error is bounded by
 175                                      2^(5-m) if reduce=1, and by 2^(2-m)
 176                                      otherwise */
 177       mpfr_ui_sub (c, 1, c, GMP_RNDN); /* error bounded by 2^(6-m) if reduce
 178                                           is 1, and 2^(3-m) otherwise */
 179       mpfr_sqrt (c, c, GMP_RNDN); /* the absolute error is bounded by
 180                                      2^(6-m-Exp(c)) if reduce=1, and
 181                                      2^(3-m-Exp(c)) otherwise */
 182       err = 3 + 3 * reduce - MPFR_GET_EXP (c);
 183       if (neg)
 184         MPFR_CHANGE_SIGN (c);
 185
 186       /* the absolute error on c is at most 2^(err-m), which we must put
 187          in the form 2^(EXP(c)-err). */
 188       err = MPFR_GET_EXP (c) + (mp_exp_t) m - err;
 189       if (mpfr_can_round (c, err, GMP_RNDN, GMP_RNDZ,
 190                           MPFR_PREC (y) + (rnd_mode == GMP_RNDN)))
 191         break;
 192       /* check for huge cancellation */
 193       if (err < (mp_exp_t) MPFR_PREC (y))
 194         m += MPFR_PREC (y) - err;
 195       /* Check if near 1 */
 196       if (MPFR_GET_EXP (c) == 1
 197           && MPFR_MANT (c)[MPFR_LIMB_SIZE (c)-1] == MPFR_LIMB_HIGHBIT)
 198         m += m;
 199
 200     next_step:
 201       MPFR_ZIV_NEXT (loop, m);
 202       mpfr_set_prec (c, m);
 203     }
 204   MPFR_ZIV_FREE (loop);
 205
 206   inexy = mpfr_set (y, c, rnd_mode);
 207   inexz = mpfr_set (z, xr, rnd_mode);
 208
 209   mpfr_clear (c);
 210   mpfr_clear (xr);
 211
 212  end:
 213   MPFR_SAVE_EXPO_FREE (expo);
 214   mpfr_check_range (y, inexy, rnd_mode);
 215   mpfr_check_range (z, inexz, rnd_mode);
 216   MPFR_RET (1); /* Always inexact */
 217 }