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 (
86 y
, x
, -2 * expx
, 2, 0, rnd_mode
,
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
,
96 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, __gmpfr_flags
);
100 /* if we go here, one of the two MPFR_FAST_COMPUTE_IF_SMALL_INPUT
103 else /* y and x are the same variable: try to compute z first, which
104 necessarily differs */
106 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (
107 z
, __gmpfr_one
, -2 * expx
, 1, 0, rnd_mode
,
109 goto 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
,
117 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, __gmpfr_flags
);
127 MPFR_ZIV_INIT (loop
, m
);
130 /* the following is copied from sin.c */
131 if (expx
>= 2) /* reduce the argument */
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
);
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
)
149 else /* the input argument is already reduced */
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
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
)))
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)
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
);
184 MPFR_CHANGE_SIGN (c
);
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
)))
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
)
201 MPFR_ZIV_NEXT (loop
, m
);
202 mpfr_set_prec (c
, m
);
204 MPFR_ZIV_FREE (loop
);
206 inexy
= mpfr_set (y
, c
, rnd_mode
);
207 inexz
= mpfr_set (z
, xr
, rnd_mode
);
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 */