1 /* mpfr_cos -- cosine of a floating-point number 2 3 Copyright 2001-2020 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 than 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[%Pu]=%*.Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), 145 ("y[%Pu]=%*.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 to 267 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 } 302