1 /* mpfr_sin -- sine of a floating-point number
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.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
27 mpfr_sin_fast (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
31 inex = mpfr_sincos_fast (y, NULL, x, rnd_mode);
32 inex = inex & 3; /* 0: exact, 1: rounded up, 2: rounded down */
33 return (inex == 2) ? -1 : inex;
37 mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
43 int inexact, sign, reduce;
45 MPFR_SAVE_EXPO_DECL (expo);
48 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
49 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y,
52 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
54 if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
62 MPFR_ASSERTD (MPFR_IS_ZERO (x));
64 MPFR_SET_SAME_SIGN (y, x);
69 /* sin(x) = x - x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
70 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2 * MPFR_GET_EXP (x), 2, 0,
73 MPFR_SAVE_EXPO_MARK (expo);
75 /* Compute initial precision */
76 precy = MPFR_PREC (y);
78 if (precy >= MPFR_SINCOS_THRESHOLD)
79 return mpfr_sin_fast (y, x, rnd_mode);
81 m = precy + MPFR_INT_CEIL_LOG2 (precy) + 13;
82 expx = MPFR_GET_EXP (x);
87 MPFR_ZIV_INIT (loop, m);
90 /* first perform argument reduction modulo 2*Pi (if needed),
91 also helps to determine the sign of sin(x) */
92 if (expx >= 2) /* If Pi < x < 4, we need to reduce too, to determine
93 the sign of sin(x). For 2 <= |x| < Pi, we could avoid
97 /* As expx + m - 1 will silently be converted into mpfr_prec_t
98 in the mpfr_set_prec call, the assert below may be useful to
99 avoid undefined behavior. */
100 MPFR_ASSERTN (expx + m - 1 <= MPFR_PREC_MAX);
101 mpfr_set_prec (c, expx + m - 1);
102 mpfr_set_prec (xr, m);
103 mpfr_const_pi (c, MPFR_RNDN);
104 mpfr_mul_2ui (c, c, 1, MPFR_RNDN);
105 mpfr_remainder (xr, x, c, MPFR_RNDN);
106 /* The analysis is similar to that of cos.c:
107 |xr - x - 2kPi| <= 2^(2-m). Thus we can decide the sign
108 of sin(x) if xr is at distance at least 2^(2-m) of both
110 mpfr_div_2ui (c, c, 1, MPFR_RNDN);
111 /* Since c approximates Pi with an error <= 2^(2-expx-m) <= 2^(-m),
112 it suffices to check that c - |xr| >= 2^(2-m). */
113 if (MPFR_SIGN (xr) > 0)
114 mpfr_sub (c, c, xr, MPFR_RNDZ);
116 mpfr_add (c, c, xr, MPFR_RNDZ);
118 || MPFR_GET_EXP(xr) < (mpfr_exp_t) 3 - (mpfr_exp_t) m
120 || MPFR_GET_EXP(c) < (mpfr_exp_t) 3 - (mpfr_exp_t) m)
123 /* |xr - x - 2kPi| <= 2^(2-m), thus |sin(xr) - sin(x)| <= 2^(2-m) */
126 else /* the input argument is already reduced */
132 sign = MPFR_SIGN(xx);
133 /* now that the argument is reduced, precision m is enough */
134 mpfr_set_prec (c, m);
135 mpfr_cos (c, xx, MPFR_RNDZ); /* can't be exact */
136 mpfr_nexttoinf (c); /* now c = cos(x) rounded away */
137 mpfr_mul (c, c, c, MPFR_RNDU); /* away */
138 mpfr_ui_sub (c, 1, c, MPFR_RNDZ);
139 mpfr_sqrt (c, c, MPFR_RNDZ);
140 if (MPFR_IS_NEG_SIGN(sign))
143 /* Warning: c may be 0! */
144 if (MPFR_UNLIKELY (MPFR_IS_ZERO (c)))
146 /* Huge cancellation: increase prec a lot! */
147 m = MAX (m, MPFR_PREC (x));
152 /* the absolute error on c is at most 2^(3-m-EXP(c)),
153 plus 2^(2-m) if there was an argument reduction.
154 Since EXP(c) <= 1, 3-m-EXP(c) >= 2-m, thus the error
155 is at most 2^(3-m-EXP(c)) in case of argument reduction. */
156 err = 2 * MPFR_GET_EXP (c) + (mpfr_exp_t) m - 3 - (reduce != 0);
157 if (MPFR_CAN_ROUND (c, err, precy, rnd_mode))
160 /* check for huge cancellation (Near 0) */
161 if (err < (mpfr_exp_t) MPFR_PREC (y))
162 m += MPFR_PREC (y) - err;
163 /* Check if near 1 */
164 if (MPFR_GET_EXP (c) == 1)
169 /* Else generic increase */
170 MPFR_ZIV_NEXT (loop, m);
172 MPFR_ZIV_FREE (loop);
174 inexact = mpfr_set (y, c, rnd_mode);
175 /* inexact cannot be 0, since this would mean that c was representable
176 within the target precision, but in that case mpfr_can_round will fail */
181 MPFR_SAVE_EXPO_FREE (expo);
182 return mpfr_check_range (y, inexact, rnd_mode);