1 /* mpfr_sin -- sine of a floating-point number
3 Copyright 2001-2015 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
)
80 inexact
= mpfr_sin_fast (y
, x
, rnd_mode
);
84 m
= precy
+ MPFR_INT_CEIL_LOG2 (precy
) + 13;
85 expx
= MPFR_GET_EXP (x
);
90 MPFR_ZIV_INIT (loop
, m
);
93 /* first perform argument reduction modulo 2*Pi (if needed),
94 also helps to determine the sign of sin(x) */
95 if (expx
>= 2) /* If Pi < x < 4, we need to reduce too, to determine
96 the sign of sin(x). For 2 <= |x| < Pi, we could avoid
100 /* As expx + m - 1 will silently be converted into mpfr_prec_t
101 in the mpfr_set_prec call, the assert below may be useful to
102 avoid undefined behavior. */
103 MPFR_ASSERTN (expx
+ m
- 1 <= MPFR_PREC_MAX
);
104 mpfr_set_prec (c
, expx
+ m
- 1);
105 mpfr_set_prec (xr
, m
);
106 mpfr_const_pi (c
, MPFR_RNDN
);
107 mpfr_mul_2ui (c
, c
, 1, MPFR_RNDN
);
108 mpfr_remainder (xr
, x
, c
, MPFR_RNDN
);
109 /* The analysis is similar to that of cos.c:
110 |xr - x - 2kPi| <= 2^(2-m). Thus we can decide the sign
111 of sin(x) if xr is at distance at least 2^(2-m) of both
113 mpfr_div_2ui (c
, c
, 1, MPFR_RNDN
);
114 /* Since c approximates Pi with an error <= 2^(2-expx-m) <= 2^(-m),
115 it suffices to check that c - |xr| >= 2^(2-m). */
116 if (MPFR_SIGN (xr
) > 0)
117 mpfr_sub (c
, c
, xr
, MPFR_RNDZ
);
119 mpfr_add (c
, c
, xr
, MPFR_RNDZ
);
121 || MPFR_GET_EXP(xr
) < (mpfr_exp_t
) 3 - (mpfr_exp_t
) m
123 || MPFR_GET_EXP(c
) < (mpfr_exp_t
) 3 - (mpfr_exp_t
) m
)
126 /* |xr - x - 2kPi| <= 2^(2-m), thus |sin(xr) - sin(x)| <= 2^(2-m) */
129 else /* the input argument is already reduced */
135 sign
= MPFR_SIGN(xx
);
136 /* now that the argument is reduced, precision m is enough */
137 mpfr_set_prec (c
, m
);
138 mpfr_cos (c
, xx
, MPFR_RNDZ
); /* can't be exact */
139 mpfr_nexttoinf (c
); /* now c = cos(x) rounded away */
140 mpfr_mul (c
, c
, c
, MPFR_RNDU
); /* away */
141 mpfr_ui_sub (c
, 1, c
, MPFR_RNDZ
);
142 mpfr_sqrt (c
, c
, MPFR_RNDZ
);
143 if (MPFR_IS_NEG_SIGN(sign
))
146 /* Warning: c may be 0! */
147 if (MPFR_UNLIKELY (MPFR_IS_ZERO (c
)))
149 /* Huge cancellation: increase prec a lot! */
150 m
= MAX (m
, MPFR_PREC (x
));
155 /* the absolute error on c is at most 2^(3-m-EXP(c)),
156 plus 2^(2-m) if there was an argument reduction.
157 Since EXP(c) <= 1, 3-m-EXP(c) >= 2-m, thus the error
158 is at most 2^(3-m-EXP(c)) in case of argument reduction. */
159 err
= 2 * MPFR_GET_EXP (c
) + (mpfr_exp_t
) m
- 3 - (reduce
!= 0);
160 if (MPFR_CAN_ROUND (c
, err
, precy
, rnd_mode
))
163 /* check for huge cancellation (Near 0) */
164 if (err
< (mpfr_exp_t
) MPFR_PREC (y
))
165 m
+= MPFR_PREC (y
) - err
;
166 /* Check if near 1 */
167 if (MPFR_GET_EXP (c
) == 1)
172 /* Else generic increase */
173 MPFR_ZIV_NEXT (loop
, m
);
175 MPFR_ZIV_FREE (loop
);
177 inexact
= mpfr_set (y
, c
, rnd_mode
);
178 /* inexact cannot be 0, since this would mean that c was representable
179 within the target precision, but in that case mpfr_can_round will fail */
185 MPFR_SAVE_EXPO_FREE (expo
);
186 return mpfr_check_range (y
, inexact
, rnd_mode
);