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 4 Copyright 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 5 Contributed by the Arenaire and Caramel projects, INRIA. 6 7 This file is part of the GNU MPFR Library. 8 9 The GNU MPFR Library is free software; you can redistribute it and/or modify 10 it under the terms of the GNU Lesser General Public License as published by 11 the Free Software Foundation; either version 3 of the License, or (at your 12 option) any later version. 13 14 The GNU MPFR Library is distributed in the hope that it will be useful, but 15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17 License for more details. 18 19 You should have received a copy of the GNU Lesser General Public License 20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 23 24 #define MPFR_NEED_LONGLONG_H 25 #include "mpfr-impl.h" 26 27 static int mpfr_yn_asympt (mpfr_ptr, long, mpfr_srcptr, mpfr_rnd_t); 28 29 int 30 mpfr_y0 (mpfr_ptr res, mpfr_srcptr z, mpfr_rnd_t r) 31 { 32 return mpfr_yn (res, 0, z, r); 33 } 34 35 int 36 mpfr_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 */ 45 static mpfr_exp_t 46 mpfr_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 */ 84 static mpfr_exp_t 85 mpfr_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 149 int 150 mpfr_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_SIGN(z) > 0) 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_SIGN(z) < 0) 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