178556Sobrien/* mpfr_cos -- cosine of a floating-point number 278556Sobrien 378556SobrienCopyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 478556SobrienContributed by the AriC and Caramel projects, INRIA. 578556Sobrien 678556SobrienThis file is part of the GNU MPFR Library. 7167974Sdelphij 8167974SdelphijThe GNU MPFR Library is free software; you can redistribute it and/or modify 9167974Sdelphijit under the terms of the GNU Lesser General Public License as published by 1078556Sobrienthe Free Software Foundation; either version 3 of the License, or (at your 11215041Sobrienoption) any later version. 12215041Sobrien 1378556SobrienThe GNU MPFR Library is distributed in the hope that it will be useful, but 14167974SdelphijWITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15167974Sdelphijor FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 1678556SobrienLicense for more details. 17167974Sdelphij 18167974SdelphijYou should have received a copy of the GNU Lesser General Public License 19167974Sdelphijalong with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 2078556Sobrienhttp://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2178556Sobrien51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 2278556Sobrien 2378556Sobrien#define MPFR_NEED_LONGLONG_H 2478556Sobrien#include "mpfr-impl.h" 2578556Sobrien 2678556Sobrienstatic int 2778556Sobrienmpfr_cos_fast (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 2878556Sobrien{ 2978556Sobrien int inex; 3078556Sobrien 3178556Sobrien inex = mpfr_sincos_fast (NULL, y, x, rnd_mode); 3278556Sobrien inex = inex >> 2; /* 0: exact, 1: rounded up, 2: rounded down */ 3378556Sobrien return (inex == 2) ? -1 : inex; 3478556Sobrien} 3578556Sobrien 3678556Sobrien/* f <- 1 - r/2! + r^2/4! + ... + (-1)^l r^l/(2l)! + ... 3778556Sobrien Assumes |r| < 1/2, and f, r have the same precision. 3878556Sobrien Returns e such that the error on f is bounded by 2^e ulps. 3978556Sobrien*/ 4078556Sobrienstatic int 4178556Sobrienmpfr_cos2_aux (mpfr_ptr f, mpfr_srcptr r) 4278556Sobrien{ 4378556Sobrien mpz_t x, t, s; 4478556Sobrien mpfr_exp_t ex, l, m; 4578556Sobrien mpfr_prec_t p, q; 4678556Sobrien unsigned long i, maxi, imax; 4778556Sobrien 4878556Sobrien MPFR_ASSERTD(mpfr_get_exp (r) <= -1); 4978556Sobrien 5078556Sobrien /* compute minimal i such that i*(i+1) does not fit in an unsigned long, 5178556Sobrien assuming that there are no padding bits. */ 5278556Sobrien maxi = 1UL << (CHAR_BIT * sizeof(unsigned long) / 2); 5378556Sobrien if (maxi * (maxi / 2) == 0) /* test checked at compile time */ 5478556Sobrien { 5578556Sobrien /* can occur only when there are padding bits. */ 5678556Sobrien /* maxi * (maxi-1) is representable iff maxi * (maxi / 2) != 0 */ 5778556Sobrien do 5878556Sobrien maxi /= 2; 5978556Sobrien while (maxi * (maxi / 2) == 0); 6078556Sobrien } 6178556Sobrien 6278556Sobrien mpz_init (x); 6378556Sobrien mpz_init (s); 6478556Sobrien mpz_init (t); 6578556Sobrien ex = mpfr_get_z_2exp (x, r); /* r = x*2^ex */ 6678556Sobrien 6778556Sobrien /* remove trailing zeroes */ 6878556Sobrien l = mpz_scan1 (x, 0); 6978556Sobrien ex += l; 7078556Sobrien mpz_fdiv_q_2exp (x, x, l); 7178556Sobrien 7278556Sobrien /* since |r| < 1, r = x*2^ex, and x is an integer, necessarily ex < 0 */ 7378556Sobrien 7478556Sobrien p = mpfr_get_prec (f); /* same than r */ 7578556Sobrien /* bound for number of iterations */ 7678556Sobrien imax = p / (-mpfr_get_exp (r)); 7778556Sobrien imax += (imax == 0); 7878556Sobrien q = 2 * MPFR_INT_CEIL_LOG2(imax) + 4; /* bound for (3l)^2 */ 7978556Sobrien 8078556Sobrien mpz_set_ui (s, 1); /* initialize sum with 1 */ 8178556Sobrien mpz_mul_2exp (s, s, p + q); /* scale all values by 2^(p+q) */ 8278556Sobrien mpz_set (t, s); /* invariant: t is previous term */ 8378556Sobrien for (i = 1; (m = mpz_sizeinbase (t, 2)) >= q; i += 2) 8478556Sobrien { 8578556Sobrien /* adjust precision of x to that of t */ 8678556Sobrien l = mpz_sizeinbase (x, 2); 8778556Sobrien if (l > m) 8878556Sobrien { 8978556Sobrien l -= m; 9078556Sobrien mpz_fdiv_q_2exp (x, x, l); 9178556Sobrien ex += l; 9278556Sobrien } 9378556Sobrien /* multiply t by r */ 9478556Sobrien mpz_mul (t, t, x); 9578556Sobrien mpz_fdiv_q_2exp (t, t, -ex); 9678556Sobrien /* divide t by i*(i+1) */ 9778556Sobrien if (i < maxi) 9878556Sobrien mpz_fdiv_q_ui (t, t, i * (i + 1)); 9978556Sobrien else 10078556Sobrien { 10178556Sobrien mpz_fdiv_q_ui (t, t, i); 10278556Sobrien mpz_fdiv_q_ui (t, t, i + 1); 10378556Sobrien } 10478556Sobrien /* if m is the (current) number of bits of t, we can consider that 10578556Sobrien all operations on t so far had precision >= m, so we can prove 10678556Sobrien by induction that the relative error on t is of the form 10778556Sobrien (1+u)^(3l)-1, where |u| <= 2^(-m), and l=(i+1)/2 is the # of loops. 10878556Sobrien Since |(1+x^2)^(1/x) - 1| <= 4x/3 for |x| <= 1/2, 10978556Sobrien for |u| <= 1/(3l)^2, the absolute error is bounded by 11078556Sobrien 4/3*(3l)*2^(-m)*t <= 4*l since |t| < 2^m. 11178556Sobrien Therefore the error on s is bounded by 2*l*(l+1). */ 11278556Sobrien /* add or subtract to s */ 11378556Sobrien if (i % 4 == 1) 11478556Sobrien mpz_sub (s, s, t); 11578556Sobrien else 11678556Sobrien mpz_add (s, s, t); 11778556Sobrien } 11878556Sobrien 11978556Sobrien mpfr_set_z (f, s, MPFR_RNDN); 12078556Sobrien mpfr_div_2ui (f, f, p + q, MPFR_RNDN); 12178556Sobrien 12278556Sobrien mpz_clear (x); 12378556Sobrien mpz_clear (s); 12478556Sobrien mpz_clear (t); 125146293Sobrien 126146293Sobrien l = (i - 1) / 2; /* number of iterations */ 127146293Sobrien return 2 * MPFR_INT_CEIL_LOG2 (l + 1) + 1; /* bound is 2l(l+1) */ 128146293Sobrien} 129146293Sobrien 130146293Sobrienint 131146293Sobrienmpfr_cos (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 132146293Sobrien{ 133146293Sobrien mpfr_prec_t K0, K, precy, m, k, l; 134146293Sobrien int inexact, reduce = 0; 135146293Sobrien mpfr_t r, s, xr, c; 136146293Sobrien mpfr_exp_t exps, cancel = 0, expx; 137146293Sobrien MPFR_ZIV_DECL (loop); 138146293Sobrien MPFR_SAVE_EXPO_DECL (expo); 139146293Sobrien MPFR_GROUP_DECL (group); 140146293Sobrien 141146293Sobrien MPFR_LOG_FUNC ( 142146293Sobrien ("x[%Pu]=%*.Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), 14378556Sobrien ("y[%Pu]=%*.Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, 14478556Sobrien inexact)); 14578556Sobrien 14678556Sobrien if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 14778556Sobrien { 14878556Sobrien if (MPFR_IS_NAN (x) || MPFR_IS_INF (x)) 14978556Sobrien { 15078556Sobrien MPFR_SET_NAN (y); 15178556Sobrien MPFR_RET_NAN; 15278556Sobrien } 15378556Sobrien else 15478556Sobrien { 15578556Sobrien MPFR_ASSERTD (MPFR_IS_ZERO (x)); 15678556Sobrien return mpfr_set_ui (y, 1, rnd_mode); 15778556Sobrien } 15878556Sobrien } 15978556Sobrien 16078556Sobrien MPFR_SAVE_EXPO_MARK (expo); 16178556Sobrien 16278556Sobrien /* cos(x) = 1-x^2/2 + ..., so error < 2^(2*EXP(x)-1) */ 16378556Sobrien expx = MPFR_GET_EXP (x); 16478556Sobrien MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, __gmpfr_one, -2 * expx, 16578556Sobrien 1, 0, rnd_mode, expo, {}); 16678556Sobrien 16778556Sobrien /* Compute initial precision */ 16878556Sobrien precy = MPFR_PREC (y); 16978556Sobrien 17078556Sobrien if (precy >= MPFR_SINCOS_THRESHOLD) 17178556Sobrien { 17278556Sobrien MPFR_SAVE_EXPO_FREE (expo); 17378556Sobrien return mpfr_cos_fast (y, x, rnd_mode); 17478556Sobrien } 17578556Sobrien 17678556Sobrien K0 = __gmpfr_isqrt (precy / 3); 17778556Sobrien m = precy + 2 * MPFR_INT_CEIL_LOG2 (precy) + 2 * K0; 17878556Sobrien 17978556Sobrien if (expx >= 3) 18078556Sobrien { 18178556Sobrien reduce = 1; 18278556Sobrien /* As expx + m - 1 will silently be converted into mpfr_prec_t 18378556Sobrien in the mpfr_init2 call, the assert below may be useful to 18478556Sobrien avoid undefined behavior. */ 18578556Sobrien MPFR_ASSERTN (expx + m - 1 <= MPFR_PREC_MAX); 18678556Sobrien mpfr_init2 (c, expx + m - 1); 18778556Sobrien mpfr_init2 (xr, m); 18878556Sobrien } 18978556Sobrien 19078556Sobrien MPFR_GROUP_INIT_2 (group, m, r, s); 19178556Sobrien MPFR_ZIV_INIT (loop, m); 19278556Sobrien for (;;) 19378556Sobrien { 19478556Sobrien /* If |x| >= 4, first reduce x cmod (2*Pi) into xr, using mpfr_remainder: 19578556Sobrien let e = EXP(x) >= 3, and m the target precision: 19678556Sobrien (1) c <- 2*Pi [precision e+m-1, nearest] 19778556Sobrien (2) xr <- remainder (x, c) [precision m, nearest] 19878556Sobrien We have |c - 2*Pi| <= 1/2ulp(c) = 2^(3-e-m) 19978556Sobrien |xr - x - k c| <= 1/2ulp(xr) <= 2^(1-m) 20078556Sobrien |k| <= |x|/(2*Pi) <= 2^(e-2) 20178556Sobrien Thus |xr - x - 2kPi| <= |k| |c - 2Pi| + 2^(1-m) <= 2^(2-m). 20278556Sobrien It follows |cos(xr) - cos(x)| <= 2^(2-m). */ 20378556Sobrien if (reduce) 20478556Sobrien { 20578556Sobrien 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