yn.c revision 1.1.1.3
1/* mpfr_y0, mpfr_y1, mpfr_yn -- Bessel functions of 2nd kind, integer order. 2 http://www.opengroup.org/onlinepubs/009695399/functions/y0.html 3 4Copyright 2007-2018 Free Software Foundation, Inc. 5Contributed by the AriC and Caramba projects, INRIA. 6 7This file is part of the GNU MPFR Library. 8 9The GNU MPFR Library is free software; you can redistribute it and/or modify 10it under the terms of the GNU Lesser General Public License as published by 11the Free Software Foundation; either version 3 of the License, or (at your 12option) any later version. 13 14The GNU MPFR Library is distributed in the hope that it will be useful, but 15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17License for more details. 18 19You should have received a copy of the GNU Lesser General Public License 20along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 21http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2251 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 23 24#define MPFR_NEED_LONGLONG_H 25#include "mpfr-impl.h" 26 27static int mpfr_yn_asympt (mpfr_ptr, long, mpfr_srcptr, mpfr_rnd_t); 28 29int 30mpfr_y0 (mpfr_ptr res, mpfr_srcptr z, mpfr_rnd_t r) 31{ 32 return mpfr_yn (res, 0, z, r); 33} 34 35int 36mpfr_y1 (mpfr_ptr res, mpfr_srcptr z, mpfr_rnd_t r) 37{ 38 return mpfr_yn (res, 1, z, r); 39} 40 41/* compute in s an approximation of S1 = sum((n-k)!/k!*y^k,k=0..n) 42 return e >= 0 the exponent difference between the maximal value of |s| 43 during the for loop and the final value of |s|. 44*/ 45static mpfr_exp_t 46mpfr_yn_s1 (mpfr_ptr s, mpfr_srcptr y, unsigned long n) 47{ 48 unsigned long k; 49 mpz_t f; 50 mpfr_exp_t e, emax; 51 52 mpz_init_set_ui (f, 1); 53 /* we compute n!*S1 = sum(a[k]*y^k,k=0..n) where a[k] = n!*(n-k)!/k!, 54 a[0] = (n!)^2, a[1] = n!*(n-1)!, ..., a[n-1] = n, a[n] = 1 */ 55 mpfr_set_ui (s, 1, MPFR_RNDN); /* a[n] */ 56 emax = MPFR_EXP(s); 57 for (k = n; k-- > 0;) 58 { 59 /* a[k]/a[k+1] = (n-k)!/k!/(n-(k+1))!*(k+1)! = (k+1)*(n-k) */ 60 mpfr_mul (s, s, y, MPFR_RNDN); 61 mpz_mul_ui (f, f, n - k); 62 mpz_mul_ui (f, f, k + 1); 63 /* invariant: f = a[k] */ 64 mpfr_add_z (s, s, f, MPFR_RNDN); 65 e = MPFR_EXP(s); 66 if (e > emax) 67 emax = e; 68 } 69 /* now we have f = (n!)^2 */ 70 mpz_sqrt (f, f); 71 mpfr_div_z (s, s, f, MPFR_RNDN); 72 mpz_clear (f); 73 return emax - MPFR_EXP(s); 74} 75 76/* compute in s an approximation of 77 S3 = c*sum((h(k)+h(n+k))*y^k/k!/(n+k)!,k=0..infinity) 78 where h(k) = 1 + 1/2 + ... + 1/k 79 k=0: h(n) 80 k=1: 1+h(n+1) 81 k=2: 3/2+h(n+2) 82 Returns e such that the error is bounded by 2^e ulp(s). 83*/ 84static mpfr_exp_t 85mpfr_yn_s3 (mpfr_ptr s, mpfr_srcptr y, mpfr_srcptr c, unsigned long n) 86{ 87 unsigned long k, zz; 88 mpfr_t t, u; 89 mpz_t p, q; /* p/q will store h(k)+h(n+k) */ 90 mpfr_exp_t exps, expU; 91 92 zz = mpfr_get_ui (y, MPFR_RNDU); /* y = z^2/4 */ 93 MPFR_ASSERTN (zz < ULONG_MAX - 2); 94 zz += 2; /* z^2 <= 2^zz */ 95 mpz_init_set_ui (p, 0); 96 mpz_init_set_ui (q, 1); 97 /* initialize p/q to h(n) */ 98 for (k = 1; k <= n; k++) 99 { 100 /* p/q + 1/k = (k*p+q)/(q*k) */ 101 mpz_mul_ui (p, p, k); 102 mpz_add (p, p, q); 103 mpz_mul_ui (q, q, k); 104 } 105 mpfr_init2 (t, MPFR_PREC(s)); 106 mpfr_init2 (u, MPFR_PREC(s)); 107 mpfr_fac_ui (t, n, MPFR_RNDN); 108 mpfr_div (t, c, t, MPFR_RNDN); /* c/n! */ 109 mpfr_mul_z (u, t, p, MPFR_RNDN); 110 mpfr_div_z (s, u, q, MPFR_RNDN); 111 exps = MPFR_EXP (s); 112 expU = exps; 113 for (k = 1; ;k ++) 114 { 115 /* update t */ 116 mpfr_mul (t, t, y, MPFR_RNDN); 117 mpfr_div_ui (t, t, k, MPFR_RNDN); 118 mpfr_div_ui (t, t, n + k, MPFR_RNDN); 119 /* update p/q: 120 p/q + 1/k + 1/(n+k) = [p*k*(n+k) + q*(n+k) + q*k]/(q*k*(n+k)) */ 121 mpz_mul_ui (p, p, k); 122 mpz_mul_ui (p, p, n + k); 123 mpz_addmul_ui (p, q, n + 2 * k); 124 mpz_mul_ui (q, q, k); 125 mpz_mul_ui (q, q, n + k); 126 mpfr_mul_z (u, t, p, MPFR_RNDN); 127 mpfr_div_z (u, u, q, MPFR_RNDN); 128 exps = MPFR_EXP (u); 129 if (exps > expU) 130 expU = exps; 131 mpfr_add (s, s, u, MPFR_RNDN); 132 exps = MPFR_EXP (s); 133 if (exps > expU) 134 expU = exps; 135 if (MPFR_EXP (u) + (mpfr_exp_t) MPFR_PREC (u) < MPFR_EXP (s) && 136 zz / (2 * k) < k + n) 137 break; 138 } 139 mpfr_clear (t); 140 mpfr_clear (u); 141 mpz_clear (p); 142 mpz_clear (q); 143 exps = expU - MPFR_EXP (s); 144 /* the error is bounded by (6k^2+33/2k+11) 2^exps ulps 145 <= 8*(k+2)^2 2^exps ulps */ 146 return 3 + 2 * MPFR_INT_CEIL_LOG2(k + 2) + exps; 147} 148 149int 150mpfr_yn (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r) 151{ 152 int inex; 153 unsigned long absn; 154 MPFR_SAVE_EXPO_DECL (expo); 155 156 MPFR_LOG_FUNC 157 (("n=%ld x[%Pu]=%.*Rg rnd=%d", n, mpfr_get_prec (z), mpfr_log_prec, z, r), 158 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (res), mpfr_log_prec, res, inex)); 159 160 absn = SAFE_ABS (unsigned long, n); 161 162 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (z))) 163 { 164 if (MPFR_IS_NAN (z)) 165 { 166 MPFR_SET_NAN (res); /* y(n,NaN) = NaN */ 167 MPFR_RET_NAN; 168 } 169 /* y(n,z) tends to zero when z goes to +Inf, oscillating around 170 0. We choose to return +0 in that case. */ 171 else if (MPFR_IS_INF (z)) 172 { 173 if (MPFR_IS_POS (z)) 174 return mpfr_set_ui (res, 0, r); 175 else /* y(n,-Inf) = NaN */ 176 { 177 MPFR_SET_NAN (res); 178 MPFR_RET_NAN; 179 } 180 } 181 else /* y(n,z) tends to -Inf for n >= 0 or n even, to +Inf otherwise, 182 when z goes to zero */ 183 { 184 MPFR_SET_INF(res); 185 if (n >= 0 || ((unsigned long) n & 1) == 0) 186 MPFR_SET_NEG(res); 187 else 188 MPFR_SET_POS(res); 189 MPFR_SET_DIVBY0 (); 190 MPFR_RET(0); 191 } 192 } 193 194 /* for z < 0, y(n,z) is imaginary except when j(n,|z|) = 0, which we 195 assume does not happen for a rational z. */ 196 if (MPFR_IS_NEG (z)) 197 { 198 MPFR_SET_NAN (res); 199 MPFR_RET_NAN; 200 } 201 202 /* now z is not singular, and z > 0 */ 203 204 MPFR_SAVE_EXPO_MARK (expo); 205 206 /* Deal with tiny arguments. We have: 207 y0(z) = 2 log(z)/Pi + 2 (euler - log(2))/Pi + O(log(z)*z^2), more 208 precisely for 0 <= z <= 1/2, with g(z) = 2/Pi + 2(euler-log(2))/Pi/log(z), 209 g(z) - 0.41*z^2 < y0(z)/log(z) < g(z) 210 thus since log(z) is negative: 211 g(z)*log(z) < y0(z) < (g(z) - z^2/2)*log(z) 212 and since |g(z)| >= 0.63 for 0 <= z <= 1/2, the relative error on 213 y0(z)/log(z) is bounded by 0.41*z^2/0.63 <= 0.66*z^2. 214 Note: we use both the main term in log(z) and the constant term, because 215 otherwise the relative error would be only in 1/log(|log(z)|). 216 */ 217 if (n == 0 && MPFR_EXP(z) < - (mpfr_exp_t) (MPFR_PREC(res) / 2)) 218 { 219 mpfr_t l, h, t, logz; 220 mpfr_prec_t prec; 221 int ok, inex2; 222 223 prec = MPFR_PREC(res) + 10; 224 mpfr_init2 (l, prec); 225 mpfr_init2 (h, prec); 226 mpfr_init2 (t, prec); 227 mpfr_init2 (logz, prec); 228 /* first enclose log(z) + euler - log(2) = log(z/2) + euler */ 229 mpfr_log (logz, z, MPFR_RNDD); /* lower bound of log(z) */ 230 mpfr_set (h, logz, MPFR_RNDU); /* exact */ 231 mpfr_nextabove (h); /* upper bound of log(z) */ 232 mpfr_const_euler (t, MPFR_RNDD); /* lower bound of euler */ 233 mpfr_add (l, logz, t, MPFR_RNDD); /* lower bound of log(z) + euler */ 234 mpfr_nextabove (t); /* upper bound of euler */ 235 mpfr_add (h, h, t, MPFR_RNDU); /* upper bound of log(z) + euler */ 236 mpfr_const_log2 (t, MPFR_RNDU); /* upper bound of log(2) */ 237 mpfr_sub (l, l, t, MPFR_RNDD); /* lower bound of log(z/2) + euler */ 238 mpfr_nextbelow (t); /* lower bound of log(2) */ 239 mpfr_sub (h, h, t, MPFR_RNDU); /* upper bound of log(z/2) + euler */ 240 mpfr_const_pi (t, MPFR_RNDU); /* upper bound of Pi */ 241 mpfr_div (l, l, t, MPFR_RNDD); /* lower bound of (log(z/2)+euler)/Pi */ 242 mpfr_nextbelow (t); /* lower bound of Pi */ 243 mpfr_div (h, h, t, MPFR_RNDD); /* upper bound of (log(z/2)+euler)/Pi */ 244 mpfr_mul_2ui (l, l, 1, MPFR_RNDD); /* lower bound on g(z)*log(z) */ 245 mpfr_mul_2ui (h, h, 1, MPFR_RNDU); /* upper bound on g(z)*log(z) */ 246 /* we now have l <= g(z)*log(z) <= h, and we need to add -z^2/2*log(z) 247 to h */ 248 mpfr_mul (t, z, z, MPFR_RNDU); /* upper bound on z^2 */ 249 /* since logz is negative, a lower bound corresponds to an upper bound 250 for its absolute value */ 251 mpfr_neg (t, t, MPFR_RNDD); 252 mpfr_div_2ui (t, t, 1, MPFR_RNDD); 253 mpfr_mul (t, t, logz, MPFR_RNDU); /* upper bound on z^2/2*log(z) */ 254 mpfr_add (h, h, t, MPFR_RNDU); 255 inex = mpfr_prec_round (l, MPFR_PREC(res), r); 256 inex2 = mpfr_prec_round (h, MPFR_PREC(res), r); 257 /* we need h=l and inex=inex2 */ 258 ok = (inex == inex2) && mpfr_equal_p (l, h); 259 if (ok) 260 mpfr_set (res, h, r); /* exact */ 261 mpfr_clear (l); 262 mpfr_clear (h); 263 mpfr_clear (t); 264 mpfr_clear (logz); 265 if (ok) 266 goto end; 267 } 268 269 /* small argument check for y1(z) = -2/Pi/z + O(log(z)): 270 for 0 <= z <= 1, |y1(z) + 2/Pi/z| <= 0.25 */ 271 if (n == 1 && MPFR_EXP(z) + 1 < - (mpfr_exp_t) MPFR_PREC(res)) 272 { 273 mpfr_t y; 274 mpfr_prec_t prec; 275 mpfr_exp_t err1; 276 int ok; 277 MPFR_BLOCK_DECL (flags); 278 279 /* since 2/Pi > 0.5, and |y1(z)| >= |2/Pi/z|, if z <= 2^(-emax-1), 280 then |y1(z)| > 2^emax */ 281 prec = MPFR_PREC(res) + 10; 282 mpfr_init2 (y, prec); 283 mpfr_const_pi (y, MPFR_RNDU); /* Pi*(1+u)^2, where here and below u 284 represents a quantity <= 1/2^prec */ 285 mpfr_mul (y, y, z, MPFR_RNDU); /* Pi*z * (1+u)^4, upper bound */ 286 MPFR_BLOCK (flags, mpfr_ui_div (y, 2, y, MPFR_RNDZ)); 287 /* 2/Pi/z * (1+u)^6, lower bound, with possible overflow */ 288 if (MPFR_OVERFLOW (flags)) 289 { 290 mpfr_clear (y); 291 MPFR_SAVE_EXPO_FREE (expo); 292 return mpfr_overflow (res, r, -1); 293 } 294 mpfr_neg (y, y, MPFR_RNDN); 295 /* (1+u)^6 can be written 1+7u [for another value of u], thus the 296 error on 2/Pi/z is less than 7ulp(y). The truncation error is less 297 than 1/4, thus if ulp(y)>=1/4, the total error is less than 8ulp(y), 298 otherwise it is less than 1/4+7/8 <= 2. */ 299 if (MPFR_EXP(y) + 2 >= MPFR_PREC(y)) /* ulp(y) >= 1/4 */ 300 err1 = 3; 301 else /* ulp(y) <= 1/8 */ 302 err1 = (mpfr_exp_t) MPFR_PREC(y) - MPFR_EXP(y) + 1; 303 ok = MPFR_CAN_ROUND (y, prec - err1, MPFR_PREC(res), r); 304 if (ok) 305 inex = mpfr_set (res, y, r); 306 mpfr_clear (y); 307 if (ok) 308 goto end; 309 } 310 311 /* we can use the asymptotic expansion as soon as z > p log(2)/2, 312 but to get some margin we use it for z > p/2 */ 313 if (mpfr_cmp_ui (z, MPFR_PREC(res) / 2 + 3) > 0) 314 { 315 inex = mpfr_yn_asympt (res, n, z, r); 316 if (inex != 0) 317 goto end; 318 } 319 320 /* General case */ 321 { 322 mpfr_prec_t prec; 323 mpfr_exp_t err1, err2, err3; 324 mpfr_t y, s1, s2, s3; 325 MPFR_ZIV_DECL (loop); 326 327 mpfr_init (y); 328 mpfr_init (s1); 329 mpfr_init (s2); 330 mpfr_init (s3); 331 332 prec = MPFR_PREC(res) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (res)) + 13; 333 MPFR_ZIV_INIT (loop, prec); 334 for (;;) 335 { 336 mpfr_set_prec (y, prec); 337 mpfr_set_prec (s1, prec); 338 mpfr_set_prec (s2, prec); 339 mpfr_set_prec (s3, prec); 340 341 mpfr_mul (y, z, z, MPFR_RNDN); 342 mpfr_div_2ui (y, y, 2, MPFR_RNDN); /* z^2/4 */ 343 344 /* store (z/2)^n temporarily in s2 */ 345 mpfr_pow_ui (s2, z, absn, MPFR_RNDN); 346 mpfr_div_2si (s2, s2, absn, MPFR_RNDN); 347 348 /* compute S1 * (z/2)^(-n) */ 349 if (n == 0) 350 { 351 mpfr_set_ui (s1, 0, MPFR_RNDN); 352 err1 = 0; 353 } 354 else 355 err1 = mpfr_yn_s1 (s1, y, absn - 1); 356 mpfr_div (s1, s1, s2, MPFR_RNDN); /* (z/2)^(-n) * S1 */ 357 /* See algorithms.tex: the relative error on s1 is bounded by 358 (3n+3)*2^(e+1-prec). */ 359 err1 = MPFR_INT_CEIL_LOG2 (3 * absn + 3) + err1 + 1; 360 /* rel_err(s1) <= 2^(err1-prec), thus err(s1) <= 2^err1 ulps */ 361 362 /* compute (z/2)^n * S3 */ 363 mpfr_neg (y, y, MPFR_RNDN); /* -z^2/4 */ 364 err3 = mpfr_yn_s3 (s3, y, s2, absn); /* (z/2)^n * S3 */ 365 /* the error on s3 is bounded by 2^err3 ulps */ 366 367 /* add s1+s3 */ 368 err1 += MPFR_EXP(s1); 369 mpfr_add (s1, s1, s3, MPFR_RNDN); 370 /* the error is bounded by 1/2 + 2^err1*2^(- EXP(s1)) 371 + 2^err3*2^(EXP(s3) - EXP(s1)) */ 372 err3 += MPFR_EXP(s3); 373 err1 = (err3 > err1) ? err3 + 1 : err1 + 1; 374 err1 -= MPFR_EXP(s1); 375 err1 = (err1 >= 0) ? err1 + 1 : 1; 376 /* now the error on s1 is bounded by 2^err1*ulp(s1) */ 377 378 /* compute S2 */ 379 mpfr_div_2ui (s2, z, 1, MPFR_RNDN); /* z/2 */ 380 mpfr_log (s2, s2, MPFR_RNDN); /* log(z/2) */ 381 mpfr_const_euler (s3, MPFR_RNDN); 382 err2 = MPFR_EXP(s2) > MPFR_EXP(s3) ? MPFR_EXP(s2) : MPFR_EXP(s3); 383 mpfr_add (s2, s2, s3, MPFR_RNDN); /* log(z/2) + gamma */ 384 err2 -= MPFR_EXP(s2); 385 mpfr_mul_2ui (s2, s2, 1, MPFR_RNDN); /* 2*(log(z/2) + gamma) */ 386 mpfr_jn (s3, absn, z, MPFR_RNDN); /* Jn(z) */ 387 mpfr_mul (s2, s2, s3, MPFR_RNDN); /* 2*(log(z/2) + gamma)*Jn(z) */ 388 err2 += 4; /* the error on s2 is bounded by 2^err2 ulps, see 389 algorithms.tex */ 390 391 /* add all three sums */ 392 err1 += MPFR_EXP(s1); /* the error on s1 is bounded by 2^err1 */ 393 err2 += MPFR_EXP(s2); /* the error on s2 is bounded by 2^err2 */ 394 mpfr_sub (s2, s2, s1, MPFR_RNDN); /* s2 - (s1+s3) */ 395 err2 = (err1 > err2) ? err1 + 1 : err2 + 1; 396 err2 -= MPFR_EXP(s2); 397 err2 = (err2 >= 0) ? err2 + 1 : 1; 398 /* now the error on s2 is bounded by 2^err2*ulp(s2) */ 399 mpfr_const_pi (y, MPFR_RNDN); /* error bounded by 1 ulp */ 400 mpfr_div (s2, s2, y, MPFR_RNDN); /* error bounded by 401 2^(err2+1)*ulp(s2) */ 402 err2 ++; 403 404 if (MPFR_LIKELY (MPFR_CAN_ROUND (s2, prec - err2, MPFR_PREC(res), r))) 405 break; 406 MPFR_ZIV_NEXT (loop, prec); 407 } 408 MPFR_ZIV_FREE (loop); 409 410 /* Assume two's complement for the test n & 1 */ 411 inex = mpfr_set4 (res, s2, r, n >= 0 || (n & 1) == 0 ? 412 MPFR_SIGN (s2) : - MPFR_SIGN (s2)); 413 414 mpfr_clear (y); 415 mpfr_clear (s1); 416 mpfr_clear (s2); 417 mpfr_clear (s3); 418 } 419 420 end: 421 MPFR_SAVE_EXPO_FREE (expo); 422 return mpfr_check_range (res, inex, r); 423} 424 425#define MPFR_YN 426#include "jyn_asympt.c" 427