1 /* mpfr_cos -- cosine of a floating-point number
2
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.
5
6 This file is part of the GNU MPFR Library.
7
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.
12
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.
17
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. */
22
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
25
26 static int
mpfr_cos_fast(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)27 mpfr_cos_fast (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
28 {
29 int inex;
30
31 inex = mpfr_sincos_fast (NULL, y, x, rnd_mode);
32 inex = inex >> 2; /* 0: exact, 1: rounded up, 2: rounded down */
33 return (inex == 2) ? -1 : inex;
34 }
35
36 /* f <- 1 - r/2! + r^2/4! + ... + (-1)^l r^l/(2l)! + ...
37 Assumes |r| < 1/2, and f, r have the same precision.
38 Returns e such that the error on f is bounded by 2^e ulps.
39 */
40 static int
mpfr_cos2_aux(mpfr_ptr f,mpfr_srcptr r)41 mpfr_cos2_aux (mpfr_ptr f, mpfr_srcptr r)
42 {
43 mpz_t x, t, s;
44 mpfr_exp_t ex, l, m;
45 mpfr_prec_t p, q;
46 unsigned long i, maxi, imax;
47
48 MPFR_ASSERTD(mpfr_get_exp (r) <= -1);
49
50 /* compute minimal i such that i*(i+1) does not fit in an unsigned long,
51 assuming that there are no padding bits. */
52 maxi = 1UL << (CHAR_BIT * sizeof(unsigned long) / 2);
53 if (maxi * (maxi / 2) == 0) /* test checked at compile time */
54 {
55 /* can occur only when there are padding bits. */
56 /* maxi * (maxi-1) is representable iff maxi * (maxi / 2) != 0 */
57 do
58 maxi /= 2;
59 while (maxi * (maxi / 2) == 0);
60 }
61
62 mpz_init (x);
63 mpz_init (s);
64 mpz_init (t);
65 ex = mpfr_get_z_2exp (x, r); /* r = x*2^ex */
66
67 /* remove trailing zeroes */
68 l = mpz_scan1 (x, 0);
69 ex += l;
70 mpz_fdiv_q_2exp (x, x, l);
71
72 /* since |r| < 1, r = x*2^ex, and x is an integer, necessarily ex < 0 */
73
74 p = mpfr_get_prec (f); /* same than r */
75 /* bound for number of iterations */
76 imax = p / (-mpfr_get_exp (r));
77 imax += (imax == 0);
78 q = 2 * MPFR_INT_CEIL_LOG2(imax) + 4; /* bound for (3l)^2 */
79
80 mpz_set_ui (s, 1); /* initialize sum with 1 */
81 mpz_mul_2exp (s, s, p + q); /* scale all values by 2^(p+q) */
82 mpz_set (t, s); /* invariant: t is previous term */
83 for (i = 1; (m = mpz_sizeinbase (t, 2)) >= q; i += 2)
84 {
85 /* adjust precision of x to that of t */
86 l = mpz_sizeinbase (x, 2);
87 if (l > m)
88 {
89 l -= m;
90 mpz_fdiv_q_2exp (x, x, l);
91 ex += l;
92 }
93 /* multiply t by r */
94 mpz_mul (t, t, x);
95 mpz_fdiv_q_2exp (t, t, -ex);
96 /* divide t by i*(i+1) */
97 if (i < maxi)
98 mpz_fdiv_q_ui (t, t, i * (i + 1));
99 else
100 {
101 mpz_fdiv_q_ui (t, t, i);
102 mpz_fdiv_q_ui (t, t, i + 1);
103 }
104 /* if m is the (current) number of bits of t, we can consider that
105 all operations on t so far had precision >= m, so we can prove
106 by induction that the relative error on t is of the form
107 (1+u)^(3l)-1, where |u| <= 2^(-m), and l=(i+1)/2 is the # of loops.
108 Since |(1+x^2)^(1/x) - 1| <= 4x/3 for |x| <= 1/2,
109 for |u| <= 1/(3l)^2, the absolute error is bounded by
110 4/3*(3l)*2^(-m)*t <= 4*l since |t| < 2^m.
111 Therefore the error on s is bounded by 2*l*(l+1). */
112 /* add or subtract to s */
113 if (i % 4 == 1)
114 mpz_sub (s, s, t);
115 else
116 mpz_add (s, s, t);
117 }
118
119 mpfr_set_z (f, s, MPFR_RNDN);
120 mpfr_div_2ui (f, f, p + q, MPFR_RNDN);
121
122 mpz_clear (x);
123 mpz_clear (s);
124 mpz_clear (t);
125
126 l = (i - 1) / 2; /* number of iterations */
127 return 2 * MPFR_INT_CEIL_LOG2 (l + 1) + 1; /* bound is 2l(l+1) */
128 }
129
130 int
mpfr_cos(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)131 mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
132 {
133 mpfr_prec_t K0, K, precy, m, k, l;
134 int inexact, reduce = 0;
135 mpfr_t r, s, xr, c;
136 mpfr_exp_t exps, cancel = 0, expx;
137 MPFR_ZIV_DECL (loop);
138 MPFR_SAVE_EXPO_DECL (expo);
139 MPFR_GROUP_DECL (group);
140
141 MPFR_LOG_FUNC (
142 ("x[%Pu]=%*.Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
143 ("y[%Pu]=%*.Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y,
144 inexact));
145
146 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
147 {
148 if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
149 {
150 MPFR_SET_NAN (y);
151 MPFR_RET_NAN;
152 }
153 else
154 {
155 MPFR_ASSERTD (MPFR_IS_ZERO (x));
156 return mpfr_set_ui (y, 1, rnd_mode);
157 }
158 }
159
160 MPFR_SAVE_EXPO_MARK (expo);
161
162 /* cos(x) = 1-x^2/2 + ..., so error < 2^(2*EXP(x)-1) */
163 expx = MPFR_GET_EXP (x);
164 MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, __gmpfr_one, -2 * expx,
165 1, 0, rnd_mode, expo, {});
166
167 /* Compute initial precision */
168 precy = MPFR_PREC (y);
169
170 if (precy >= MPFR_SINCOS_THRESHOLD)
171 {
172 MPFR_SAVE_EXPO_FREE (expo);
173 return mpfr_cos_fast (y, x, rnd_mode);
174 }
175
176 K0 = __gmpfr_isqrt (precy / 3);
177 m = precy + 2 * MPFR_INT_CEIL_LOG2 (precy) + 2 * K0;
178
179 if (expx >= 3)
180 {
181 reduce = 1;
182 /* As expx + m - 1 will silently be converted into mpfr_prec_t
183 in the mpfr_init2 call, the assert below may be useful to
184 avoid undefined behavior. */
185 MPFR_ASSERTN (expx + m - 1 <= MPFR_PREC_MAX);
186 mpfr_init2 (c, expx + m - 1);
187 mpfr_init2 (xr, m);
188 }
189
190 MPFR_GROUP_INIT_2 (group, m, r, s);
191 MPFR_ZIV_INIT (loop, m);
192 for (;;)
193 {
194 /* If |x| >= 4, first reduce x cmod (2*Pi) into xr, using mpfr_remainder:
195 let e = EXP(x) >= 3, and m the target precision:
196 (1) c <- 2*Pi [precision e+m-1, nearest]
197 (2) xr <- remainder (x, c) [precision m, nearest]
198 We have |c - 2*Pi| <= 1/2ulp(c) = 2^(3-e-m)
199 |xr - x - k c| <= 1/2ulp(xr) <= 2^(1-m)
200 |k| <= |x|/(2*Pi) <= 2^(e-2)
201 Thus |xr - x - 2kPi| <= |k| |c - 2Pi| + 2^(1-m) <= 2^(2-m).
202 It follows |cos(xr) - cos(x)| <= 2^(2-m). */
203 if (reduce)
204 {
205 mpfr_const_pi (c, MPFR_RNDN);
206 mpfr_mul_2ui (c, c, 1, MPFR_RNDN); /* 2Pi */
207 mpfr_remainder (xr, x, c, MPFR_RNDN);
208 if (MPFR_IS_ZERO(xr))
209 goto ziv_next;
210 /* now |xr| <= 4, thus r <= 16 below */
211 mpfr_mul (r, xr, xr, MPFR_RNDU); /* err <= 1 ulp */
212 }
213 else
214 mpfr_mul (r, x, x, MPFR_RNDU); /* err <= 1 ulp */
215
216 /* now |x| < 4 (or xr if reduce = 1), thus |r| <= 16 */
217
218 /* we need |r| < 1/2 for mpfr_cos2_aux, i.e., EXP(r) - 2K <= -1 */
219 K = K0 + 1 + MAX(0, MPFR_GET_EXP(r)) / 2;
220 /* since K0 >= 0, if EXP(r) < 0, then K >= 1, thus EXP(r) - 2K <= -3;
221 otherwise if EXP(r) >= 0, then K >= 1/2 + EXP(r)/2, thus
222 EXP(r) - 2K <= -1 */
223
224 MPFR_SET_EXP (r, MPFR_GET_EXP (r) - 2 * K); /* Can't overflow! */
225
226 /* s <- 1 - r/2! + ... + (-1)^l r^l/(2l)! */
227 l = mpfr_cos2_aux (s, r);
228 /* l is the error bound in ulps on s */
229 MPFR_SET_ONE (r);
230 for (k = 0; k < K; k++)
231 {
232 mpfr_sqr (s, s, MPFR_RNDU); /* err <= 2*olderr */
233 MPFR_SET_EXP (s, MPFR_GET_EXP (s) + 1); /* Can't overflow */
234 mpfr_sub (s, s, r, MPFR_RNDN); /* err <= 4*olderr */
235 if (MPFR_IS_ZERO(s))
236 goto ziv_next;
237 MPFR_ASSERTD (MPFR_GET_EXP (s) <= 1);
238 }
239
240 /* The absolute error on s is bounded by (2l+1/3)*2^(2K-m)
241 2l+1/3 <= 2l+1.
242 If |x| >= 4, we need to add 2^(2-m) for the argument reduction
243 by 2Pi: if K = 0, this amounts to add 4 to 2l+1/3, i.e., to add
244 2 to l; if K >= 1, this amounts to add 1 to 2*l+1/3. */
245 l = 2 * l + 1;
246 if (reduce)
247 l += (K == 0) ? 4 : 1;
248 k = MPFR_INT_CEIL_LOG2 (l) + 2*K;
249 /* now the error is bounded by 2^(k-m) = 2^(EXP(s)-err) */
250
251 exps = MPFR_GET_EXP (s);
252 if (MPFR_LIKELY (MPFR_CAN_ROUND (s, exps + m - k, precy, rnd_mode)))
253 break;
254
255 if (MPFR_UNLIKELY (exps == 1))
256 /* s = 1 or -1, and except x=0 which was already checked above,
257 cos(x) cannot be 1 or -1, so we can round if the error is less
258 than 2^(-precy) for directed rounding, or 2^(-precy-1) for rounding
259 to nearest. */
260 {
261 if (m > k && (m - k >= precy + (rnd_mode == MPFR_RNDN)))
262 {
263 /* If round to nearest or away, result is s = 1 or -1,
264 otherwise it is round(nexttoward (s, 0)). However in order to
265 have the inexact flag correctly set below, we set |s| to
266 1 - 2^(-m) in all cases. */
267 mpfr_nexttozero (s);
268 break;
269 }
270 }
271
272 if (exps < cancel)
273 {
274 m += cancel - exps;
275 cancel = exps;
276 }
277
278 ziv_next:
279 MPFR_ZIV_NEXT (loop, m);
280 MPFR_GROUP_REPREC_2 (group, m, r, s);
281 if (reduce)
282 {
283 mpfr_set_prec (xr, m);
284 mpfr_set_prec (c, expx + m - 1);
285 }
286 }
287 MPFR_ZIV_FREE (loop);
288 inexact = mpfr_set (y, s, rnd_mode);
289 MPFR_GROUP_CLEAR (group);
290 if (reduce)
291 {
292 mpfr_clear (xr);
293 mpfr_clear (c);
294 }
295
296 MPFR_SAVE_EXPO_FREE (expo);
297 return mpfr_check_range (y, inexact, rnd_mode);
298 }
299