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