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 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 << (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 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