1 /* mpfr_cos -- cosine 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_cos_fast (mpfr_ptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd_mode
)
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
;
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.
41 mpfr_cos2_aux (mpfr_ptr f
, mpfr_srcptr r
)
46 unsigned long i
, maxi
, imax
;
48 MPFR_ASSERTD(mpfr_get_exp (r
) <= -1);
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 */
55 /* can occur only when there are padding bits. */
56 /* maxi * (maxi-1) is representable iff maxi * (maxi / 2) != 0 */
59 while (maxi
* (maxi
/ 2) == 0);
65 ex
= mpfr_get_z_2exp (x
, r
); /* r = x*2^ex */
67 /* remove trailing zeroes */
70 mpz_fdiv_q_2exp (x
, x
, l
);
72 /* since |r| < 1, r = x*2^ex, and x is an integer, necessarily ex < 0 */
74 p
= mpfr_get_prec (f
); /* same than r */
75 /* bound for number of iterations */
76 imax
= p
/ (-mpfr_get_exp (r
));
78 q
= 2 * MPFR_INT_CEIL_LOG2(imax
) + 4; /* bound for (3l)^2 */
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)
85 /* adjust precision of x to that of t */
86 l
= mpz_sizeinbase (x
, 2);
90 mpz_fdiv_q_2exp (x
, x
, l
);
95 mpz_fdiv_q_2exp (t
, t
, -ex
);
96 /* divide t by i*(i+1) */
98 mpz_fdiv_q_ui (t
, t
, i
* (i
+ 1));
101 mpz_fdiv_q_ui (t
, t
, i
);
102 mpz_fdiv_q_ui (t
, t
, i
+ 1);
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 */
119 mpfr_set_z (f
, s
, MPFR_RNDN
);
120 mpfr_div_2ui (f
, f
, p
+ q
, MPFR_RNDN
);
126 l
= (i
- 1) / 2; /* number of iterations */
127 return 2 * MPFR_INT_CEIL_LOG2 (l
+ 1) + 1; /* bound is 2l(l+1) */
131 mpfr_cos (mpfr_ptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd_mode
)
133 mpfr_prec_t K0
, K
, precy
, m
, k
, l
;
134 int inexact
, reduce
= 0;
136 mpfr_exp_t exps
, cancel
= 0, expx
;
137 MPFR_ZIV_DECL (loop
);
138 MPFR_SAVE_EXPO_DECL (expo
);
139 MPFR_GROUP_DECL (group
);
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
,
146 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
148 if (MPFR_IS_NAN (x
) || MPFR_IS_INF (x
))
155 MPFR_ASSERTD (MPFR_IS_ZERO (x
));
156 return mpfr_set_ui (y
, 1, rnd_mode
);
160 MPFR_SAVE_EXPO_MARK (expo
);
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
, {});
167 /* Compute initial precision */
168 precy
= MPFR_PREC (y
);
170 if (precy
>= MPFR_SINCOS_THRESHOLD
)
172 inexact
= mpfr_cos_fast (y
, x
, rnd_mode
);
176 K0
= __gmpfr_isqrt (precy
/ 3);
177 m
= precy
+ 2 * MPFR_INT_CEIL_LOG2 (precy
) + 2 * K0
;
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);
190 MPFR_GROUP_INIT_2 (group
, m
, r
, s
);
191 MPFR_ZIV_INIT (loop
, m
);
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). */
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
))
210 /* now |xr| <= 4, thus r <= 16 below */
211 mpfr_mul (r
, xr
, xr
, MPFR_RNDU
); /* err <= 1 ulp */
214 mpfr_mul (r
, x
, x
, MPFR_RNDU
); /* err <= 1 ulp */
216 /* now |x| < 4 (or xr if reduce = 1), thus |r| <= 16 */
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
224 MPFR_SET_EXP (r
, MPFR_GET_EXP (r
) - 2 * K
); /* Can't overflow! */
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 */
230 for (k
= 0; k
< K
; k
++)
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 */
237 MPFR_ASSERTD (MPFR_GET_EXP (s
) <= 1);
240 /* The absolute error on s is bounded by (2l+1/3)*2^(2K-m)
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. */
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) */
251 exps
= MPFR_GET_EXP (s
);
252 if (MPFR_LIKELY (MPFR_CAN_ROUND (s
, exps
+ m
- k
, precy
, rnd_mode
)))
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
261 if (m
> k
&& (m
- k
>= precy
+ (rnd_mode
== MPFR_RNDN
)))
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. */
279 MPFR_ZIV_NEXT (loop
, m
);
280 MPFR_GROUP_REPREC_2 (group
, m
, r
, s
);
283 mpfr_set_prec (xr
, m
);
284 mpfr_set_prec (c
, expx
+ m
- 1);
287 MPFR_ZIV_FREE (loop
);
288 inexact
= mpfr_set (y
, s
, rnd_mode
);
289 MPFR_GROUP_CLEAR (group
);
297 MPFR_SAVE_EXPO_FREE (expo
);
298 return mpfr_check_range (y
, inexact
, rnd_mode
);