1 /* mpfr_sin -- sine of a floating-point number
3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC 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
);