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