1/* mpfr_cos -- cosine of a floating-point number
2
3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4Contributed by the AriC and Caramel projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23#define MPFR_NEED_LONGLONG_H
24#include "mpfr-impl.h"
25
26static int
27mpfr_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*/
40static int
41mpfr_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
130int
131mpfr_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