1 /* mpfr_sin_cos -- sine and cosine of a floating-point number
3 Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* (y, z) <- (sin(x), cos(x)), return value is 0 iff both results are exact
29 mpfr_sin_cos (mpfr_ptr y
, mpfr_ptr z
, mpfr_srcptr x
, mp_rnd_t rnd_mode
)
38 MPFR_SAVE_EXPO_DECL (expo
);
40 MPFR_ASSERTN (y
!= z
);
42 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
44 if (MPFR_IS_NAN(x
) || MPFR_IS_INF(x
))
52 MPFR_ASSERTD (MPFR_IS_ZERO (x
));
54 MPFR_SET_SAME_SIGN (y
, x
);
55 /* y = 0, thus exact, but z is inexact in case of underflow
57 return mpfr_set_ui (z
, 1, rnd_mode
);
61 MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x
, x
, rnd_mode
),
62 ("sin[%#R]=%R cos[%#R]=%R", y
, y
, z
, z
));
64 MPFR_SAVE_EXPO_MARK (expo
);
66 prec
= MAX (MPFR_PREC (y
), MPFR_PREC (z
));
67 m
= prec
+ MPFR_INT_CEIL_LOG2 (prec
) + 13;
68 expx
= MPFR_GET_EXP (x
);
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. */
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. */
83 /* y and x differ, thus we can safely try to compute y first */
85 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y
, x
, -2 * expx
, 2, 0, rnd_mode
,
91 /* we can go here only if we can round sin(x) */
92 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z
, __gmpfr_one
, -2 * expx
,
98 /* if we go here, one of the two MPFR_FAST_COMPUTE_IF_SMALL_INPUT
101 else /* y and x are the same variable: try to compute z first, which
102 necessarily differs */
104 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z
, __gmpfr_one
, -2 * expx
,
107 goto 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,
124 MPFR_ZIV_INIT (loop
, m
);
127 /* the following is copied from sin.c */
128 if (expx
>= 2) /* reduce the argument */
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
);
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
)
146 else /* the input argument is already reduced */
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
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
)))
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
);
176 MPFR_CHANGE_SIGN (c
);
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
)))
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
)
195 MPFR_ZIV_NEXT (loop
, m
);
196 mpfr_set_prec (c
, m
);
198 MPFR_ZIV_FREE (loop
);
200 inexy
= mpfr_set (y
, c
, rnd_mode
);
201 inexz
= mpfr_set (z
, xr
, rnd_mode
);
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 */