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