1 /* mpfr_jn_asympt, mpfr_yn_asympt -- shared code for mpfr_jn and mpfr_yn 2 3 Copyright 2007-2018 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 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 #ifdef MPFR_JN 24 # define FUNCTION mpfr_jn_asympt 25 #else 26 # ifdef MPFR_YN 27 # define FUNCTION mpfr_yn_asympt 28 # else 29 # error "neither MPFR_JN nor MPFR_YN is defined" 30 # endif 31 #endif 32 33 /* Implements asymptotic expansion for jn or yn (formulae 9.2.5 and 9.2.6 34 from Abramowitz & Stegun). 35 Assumes |z| > p log(2)/2, where p is the target precision 36 (z can be negative only for jn). 37 Return 0 if the expansion does not converge enough (the value 0 as inexact 38 flag should not happen for normal input). 39 Note: for MPFR_RNDF, it returns 0 if the expansion failed, and a non-zero 40 value otherwise (with no other meaning). 41 */ 42 static int 43 FUNCTION (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r) 44 { 45 mpfr_t s, c, P, Q, t, iz, err_t, err_s, err_u; 46 mpfr_prec_t w; 47 long k; 48 int inex, stop, diverge = 0; 49 mpfr_exp_t err2, err; 50 MPFR_ZIV_DECL (loop); 51 52 mpfr_init (c); 53 54 w = MPFR_PREC(res) + MPFR_INT_CEIL_LOG2(MPFR_PREC(res)) + 4; 55 56 MPFR_ZIV_INIT (loop, w); 57 for (;;) 58 { 59 mpfr_set_prec (c, w); 60 mpfr_init2 (s, w); 61 mpfr_init2 (P, w); 62 mpfr_init2 (Q, w); 63 mpfr_init2 (t, w); 64 mpfr_init2 (iz, w); 65 mpfr_init2 (err_t, 31); 66 mpfr_init2 (err_s, 31); 67 mpfr_init2 (err_u, 31); 68 69 /* Approximate sin(z) and cos(z). In the following, err <= k means that 70 the approximate value y and the true value x are related by 71 y = x * (1 + u)^k with |u| <= 2^(-w), following Higham's method. */ 72 mpfr_sin_cos (s, c, z, MPFR_RNDN); 73 if (MPFR_IS_NEG(z)) 74 mpfr_neg (s, s, MPFR_RNDN); /* compute jn/yn(|z|), fix sign later */ 75 /* The absolute error on s/c is bounded by 1/2 ulp(1/2) <= 2^(-w-1). */ 76 mpfr_add (t, s, c, MPFR_RNDN); 77 mpfr_sub (c, s, c, MPFR_RNDN); 78 mpfr_swap (s, t); 79 /* now s approximates sin(z)+cos(z), and c approximates sin(z)-cos(z), 80 with total absolute error bounded by 2^(1-w). */ 81 82 /* precompute 1/(8|z|) */ 83 mpfr_si_div (iz, MPFR_IS_POS(z) ? 1 : -1, z, MPFR_RNDN); /* err <= 1 */ 84 mpfr_div_2ui (iz, iz, 3, MPFR_RNDN); 85 86 /* compute P and Q */ 87 mpfr_set_ui (P, 1, MPFR_RNDN); 88 mpfr_set_ui (Q, 0, MPFR_RNDN); 89 mpfr_set_ui (t, 1, MPFR_RNDN); /* current term */ 90 mpfr_set_ui (err_t, 0, MPFR_RNDN); /* error on t */ 91 mpfr_set_ui (err_s, 0, MPFR_RNDN); /* error on P and Q (sum of errors) */ 92 for (k = 1, stop = 0; stop < 4; k++) 93 { 94 /* compute next term: t(k)/t(k-1) = (2n+2k-1)(2n-2k+1)/(8kz) */ 95 mpfr_mul_si (t, t, 2 * (n + k) - 1, MPFR_RNDN); /* err <= err_k + 1 */ 96 mpfr_mul_si (t, t, 2 * (n - k) + 1, MPFR_RNDN); /* err <= err_k + 2 */ 97 mpfr_div_ui (t, t, k, MPFR_RNDN); /* err <= err_k + 3 */ 98 mpfr_mul (t, t, iz, MPFR_RNDN); /* err <= err_k + 5 */ 99 /* the relative error on t is bounded by (1+u)^(5k)-1, which is 100 bounded by 6ku for 6ku <= 0.02: first |5 log(1+u)| <= |5.5u| 101 for |u| <= 0.15, then |exp(5.5u)-1| <= 6u for |u| <= 0.02. */ 102 mpfr_mul_ui (err_t, t, 6 * k, MPFR_IS_POS(t) ? MPFR_RNDU : MPFR_RNDD); 103 mpfr_abs (err_t, err_t, MPFR_RNDN); /* exact */ 104 /* the absolute error on t is bounded by err_t * 2^(-w) */ 105 mpfr_abs (err_u, t, MPFR_RNDU); 106 mpfr_mul_2ui (err_u, err_u, w, MPFR_RNDU); /* t * 2^w */ 107 mpfr_add (err_u, err_u, err_t, MPFR_RNDU); /* max|t| * 2^w */ 108 if (stop >= 2) 109 { 110 /* take into account the neglected terms: t * 2^w */ 111 mpfr_div_2ui (err_s, err_s, w, MPFR_RNDU); 112 if (MPFR_IS_POS(t)) 113 mpfr_add (err_s, err_s, t, MPFR_RNDU); 114 else 115 mpfr_sub (err_s, err_s, t, MPFR_RNDU); 116 mpfr_mul_2ui (err_s, err_s, w, MPFR_RNDU); 117 stop ++; 118 } 119 /* if k is odd, add to Q, otherwise to P */ 120 else if (k & 1) 121 { 122 /* if k = 1 mod 4, add, otherwise subtract */ 123 if ((k & 2) == 0) 124 mpfr_add (Q, Q, t, MPFR_RNDN); 125 else 126 mpfr_sub (Q, Q, t, MPFR_RNDN); 127 /* check if the next term is smaller than ulp(Q): if EXP(err_u) 128 <= EXP(Q), since the current term is bounded by 129 err_u * 2^(-w), it is bounded by ulp(Q) */ 130 if (MPFR_EXP(err_u) <= MPFR_EXP(Q)) 131 stop ++; 132 else 133 stop = 0; 134 } 135 else 136 { 137 /* if k = 0 mod 4, add, otherwise subtract */ 138 if ((k & 2) == 0) 139 mpfr_add (P, P, t, MPFR_RNDN); 140 else 141 mpfr_sub (P, P, t, MPFR_RNDN); 142 /* check if the next term is smaller than ulp(P) */ 143 if (MPFR_EXP(err_u) <= MPFR_EXP(P)) 144 stop ++; 145 else 146 stop = 0; 147 } 148 mpfr_add (err_s, err_s, err_t, MPFR_RNDU); 149 /* the sum of the rounding errors on P and Q is bounded by 150 err_s * 2^(-w) */ 151 152 /* stop when start to diverge */ 153 if (stop < 2 && 154 ((MPFR_IS_POS(z) && mpfr_cmp_ui (z, (k + 1) / 2) < 0) || 155 (MPFR_IS_NEG(z) && mpfr_cmp_si (z, - ((k + 1) / 2)) > 0))) 156 { 157 /* if we have to stop the series because it diverges, then 158 increasing the precision will most probably fail, since 159 we will stop to the same point, and thus compute a very 160 similar approximation */ 161 diverge = 1; 162 stop = 2; /* force stop */ 163 } 164 } 165 /* the sum of the total errors on P and Q is bounded by err_s * 2^(-w) */ 166 167 /* Now combine: the sum of the rounding errors on P and Q is bounded by 168 err_s * 2^(-w), and the absolute error on s/c is bounded by 2^(1-w) */ 169 if ((n & 1) == 0) /* n even: P * (sin + cos) + Q (cos - sin) for jn 170 Q * (sin + cos) + P (sin - cos) for yn */ 171 { 172 #ifdef MPFR_JN 173 mpfr_mul (c, c, Q, MPFR_RNDN); /* Q * (sin - cos) */ 174 mpfr_mul (s, s, P, MPFR_RNDN); /* P * (sin + cos) */ 175 #else 176 mpfr_mul (c, c, P, MPFR_RNDN); /* P * (sin - cos) */ 177 mpfr_mul (s, s, Q, MPFR_RNDN); /* Q * (sin + cos) */ 178 #endif 179 err = MPFR_EXP(c); 180 if (MPFR_EXP(s) > err) 181 err = MPFR_EXP(s); 182 #ifdef MPFR_JN 183 mpfr_sub (s, s, c, MPFR_RNDN); 184 #else 185 mpfr_add (s, s, c, MPFR_RNDN); 186 #endif 187 } 188 else /* n odd: P * (sin - cos) + Q (cos + sin) for jn, 189 Q * (sin - cos) - P (cos + sin) for yn */ 190 { 191 #ifdef MPFR_JN 192 mpfr_mul (c, c, P, MPFR_RNDN); /* P * (sin - cos) */ 193 mpfr_mul (s, s, Q, MPFR_RNDN); /* Q * (sin + cos) */ 194 #else 195 mpfr_mul (c, c, Q, MPFR_RNDN); /* Q * (sin - cos) */ 196 mpfr_mul (s, s, P, MPFR_RNDN); /* P * (sin + cos) */ 197 #endif 198 err = MPFR_EXP(c); 199 if (MPFR_EXP(s) > err) 200 err = MPFR_EXP(s); 201 #ifdef MPFR_JN 202 mpfr_add (s, s, c, MPFR_RNDN); 203 #else 204 mpfr_sub (s, c, s, MPFR_RNDN); 205 #endif 206 } 207 if ((n & 2) != 0) 208 mpfr_neg (s, s, MPFR_RNDN); 209 if (MPFR_EXP(s) > err) 210 err = MPFR_EXP(s); 211 /* the absolute error on s is bounded by P*err(s/c) + Q*err(s/c) 212 + err(P)*(s/c) + err(Q)*(s/c) + 3 * 2^(err - w - 1) 213 <= (|P|+|Q|) * 2^(1-w) + err_s * 2^(1-w) + 2^err * 2^(1-w), 214 since |c|, |old_s| <= 2. */ 215 err2 = (MPFR_EXP(P) >= MPFR_EXP(Q)) ? MPFR_EXP(P) + 2 : MPFR_EXP(Q) + 2; 216 /* (|P| + |Q|) * 2^(1 - w) <= 2^(err2 - w) */ 217 err = MPFR_EXP(err_s) >= err ? MPFR_EXP(err_s) + 2 : err + 2; 218 /* err_s * 2^(1-w) + 2^old_err * 2^(1-w) <= 2^err * 2^(-w) */ 219 err2 = (err >= err2) ? err + 1 : err2 + 1; 220 /* now the absolute error on s is bounded by 2^(err2 - w) */ 221 222 /* multiply by sqrt(1/(Pi*z)) */ 223 mpfr_const_pi (c, MPFR_RNDN); /* Pi, err <= 1 */ 224 mpfr_mul (c, c, z, MPFR_RNDN); /* err <= 2 */ 225 mpfr_si_div (c, MPFR_IS_POS(z) ? 1 : -1, c, MPFR_RNDN); /* err <= 3 */ 226 mpfr_sqrt (c, c, MPFR_RNDN); /* err<=5/2, thus the absolute error is 227 bounded by 3*u*|c| for |u| <= 0.25 */ 228 mpfr_mul (err_t, c, s, MPFR_SIGN(c)==MPFR_SIGN(s) ? MPFR_RNDU : MPFR_RNDD); 229 mpfr_abs (err_t, err_t, MPFR_RNDU); 230 mpfr_mul_ui (err_t, err_t, 3, MPFR_RNDU); 231 /* 3*2^(-w)*|old_c|*|s| [see below] is bounded by err_t * 2^(-w) */ 232 err2 += MPFR_EXP(c); 233 /* |old_c| * 2^(err2 - w) [see below] is bounded by 2^(err2-w) */ 234 mpfr_mul (c, c, s, MPFR_RNDN); /* the absolute error on c is bounded by 235 1/2 ulp(c) + 3*2^(-w)*|old_c|*|s| 236 + |old_c| * 2^(err2 - w) */ 237 /* compute err_t * 2^(-w) + 1/2 ulp(c) = (err_t + 2^EXP(c)) * 2^(-w) */ 238 err = (MPFR_EXP(err_t) > MPFR_EXP(c)) ? MPFR_EXP(err_t) + 1 : MPFR_EXP(c) + 1; 239 /* err_t * 2^(-w) + 1/2 ulp(c) <= 2^(err - w) */ 240 /* now err_t * 2^(-w) bounds 1/2 ulp(c) + 3*2^(-w)*|old_c|*|s| */ 241 err = (err >= err2) ? err + 1 : err2 + 1; 242 /* the absolute error on c is bounded by 2^(err - w) */ 243 244 mpfr_clear (s); 245 mpfr_clear (P); 246 mpfr_clear (Q); 247 mpfr_clear (t); 248 mpfr_clear (iz); 249 mpfr_clear (err_t); 250 mpfr_clear (err_s); 251 mpfr_clear (err_u); 252 253 err -= MPFR_EXP(c); 254 if (MPFR_LIKELY (MPFR_CAN_ROUND (c, w - err, MPFR_PREC(res), r))) 255 break; 256 if (diverge != 0) 257 { 258 MPFR_ZIV_FREE (loop); 259 mpfr_clear (c); 260 return 0; /* means that the asymptotic expansion failed */ 261 } 262 MPFR_ZIV_NEXT (loop, w); 263 } 264 MPFR_ZIV_FREE (loop); 265 266 inex = (MPFR_IS_POS(z) || ((n & 1) == 0)) ? mpfr_set (res, c, r) 267 : mpfr_neg (res, c, r); 268 mpfr_clear (c); 269 270 /* for RNDF, mpfr_set or mpfr_neg may return 0, but if we return 0, it 271 would mean the asymptotic expansion failed, thus we return 1 instead */ 272 return (r != MPFR_RNDF) ? inex : 1; 273 } 274