1/* mpfr_j0, mpfr_j1, mpfr_jn -- Bessel functions of 1st kind, integer order. 2 http://www.opengroup.org/onlinepubs/009695399/functions/j0.html 3 4Copyright 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 5Contributed by the AriC and Caramel 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 27/* Relations: j(-n,z) = (-1)^n j(n,z) 28 j(n,-z) = (-1)^n j(n,z) 29*/ 30 31static int mpfr_jn_asympt (mpfr_ptr, long, mpfr_srcptr, mpfr_rnd_t); 32 33int 34mpfr_j0 (mpfr_ptr res, mpfr_srcptr z, mpfr_rnd_t r) 35{ 36 return mpfr_jn (res, 0, z, r); 37} 38 39int 40mpfr_j1 (mpfr_ptr res, mpfr_srcptr z, mpfr_rnd_t r) 41{ 42 return mpfr_jn (res, 1, z, r); 43} 44 45/* Estimate k1 such that z^2/4 = k1 * (k1 + n) 46 i.e., k1 = (sqrt(n^2+z^2)-n)/2 = n/2 * (sqrt(1+(z/n)^2) - 1) if n != 0. 47 Return k0 = min(2*k1/log(2), ULONG_MAX). 48*/ 49static unsigned long 50mpfr_jn_k0 (unsigned long n, mpfr_srcptr z) 51{ 52 mpfr_t t, u; 53 unsigned long k0; 54 55 mpfr_init2 (t, 32); 56 mpfr_init2 (u, 32); 57 if (n == 0) 58 { 59 mpfr_abs (t, z, MPFR_RNDN); /* t = 2*k1 */ 60 } 61 else 62 { 63 mpfr_div_ui (t, z, n, MPFR_RNDN); 64 mpfr_sqr (t, t, MPFR_RNDN); 65 mpfr_add_ui (t, t, 1, MPFR_RNDN); 66 mpfr_sqrt (t, t, MPFR_RNDN); 67 mpfr_sub_ui (t, t, 1, MPFR_RNDN); 68 mpfr_mul_ui (t, t, n, MPFR_RNDN); /* t = 2*k1 */ 69 } 70 /* the following is a 32-bit approximation to nearest to 1/log(2) */ 71 mpfr_set_str_binary (u, "1.0111000101010100011101100101001"); 72 mpfr_mul (t, t, u, MPFR_RNDN); 73 if (mpfr_fits_ulong_p (t, MPFR_RNDN)) 74 k0 = mpfr_get_ui (t, MPFR_RNDN); 75 else 76 k0 = ULONG_MAX; 77 mpfr_clear (t); 78 mpfr_clear (u); 79 return k0; 80} 81 82int 83mpfr_jn (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r) 84{ 85 int inex; 86 int exception = 0; 87 unsigned long absn; 88 mpfr_prec_t prec, pbound, err; 89 mpfr_uprec_t uprec; 90 mpfr_exp_t exps, expT, diffexp; 91 mpfr_t y, s, t, absz; 92 unsigned long k, zz, k0; 93 MPFR_GROUP_DECL(g); 94 MPFR_SAVE_EXPO_DECL (expo); 95 MPFR_ZIV_DECL (loop); 96 97 MPFR_LOG_FUNC 98 (("n=%d x[%Pu]=%.*Rg rnd=%d", n, mpfr_get_prec (z), mpfr_log_prec, z, r), 99 ("res[%Pu]=%.*Rg inexact=%d", 100 mpfr_get_prec (res), mpfr_log_prec, res, inex)); 101 102 absn = SAFE_ABS (unsigned long, n); 103 104 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (z))) 105 { 106 if (MPFR_IS_NAN (z)) 107 { 108 MPFR_SET_NAN (res); 109 MPFR_RET_NAN; 110 } 111 /* j(n,z) tends to zero when z goes to +Inf or -Inf, oscillating around 112 0. We choose to return +0 in that case. */ 113 else if (MPFR_IS_INF (z)) /* FIXME: according to j(-n,z) = (-1)^n j(n,z) 114 we might want to give a sign depending on 115 z and n */ 116 return mpfr_set_ui (res, 0, r); 117 else /* z=0: j(0,0)=1, j(n odd,+/-0) = +/-0 if n > 0, -/+0 if n < 0, 118 j(n even,+/-0) = +0 */ 119 { 120 if (n == 0) 121 return mpfr_set_ui (res, 1, r); 122 else if (absn & 1) /* n odd */ 123 return (n > 0) ? mpfr_set (res, z, r) : mpfr_neg (res, z, r); 124 else /* n even */ 125 return mpfr_set_ui (res, 0, r); 126 } 127 } 128 129 MPFR_SAVE_EXPO_MARK (expo); 130 131 /* check for tiny input for j0: j0(z) = 1 - z^2/4 + ..., more precisely 132 |j0(z) - 1| <= z^2/4 for -1 <= z <= 1. */ 133 if (n == 0) 134 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (res, __gmpfr_one, -2 * MPFR_GET_EXP (z), 135 2, 0, r, inex = _inexact; goto end); 136 137 /* idem for j1: j1(z) = z/2 - z^3/16 + ..., more precisely 138 |j1(z) - z/2| <= |z^3|/16 for -1 <= z <= 1, with the sign of j1(z) - z/2 139 being the opposite of that of z. */ 140 /* TODO: add a test to trigger an error when 141 inex = _inexact; goto end 142 is forgotten in MPFR_FAST_COMPUTE_IF_SMALL_INPUT below. */ 143 if (n == 1) 144 { 145 /* We first compute 2j1(z) = z - z^3/8 + ..., then divide by 2 using 146 the "extra" argument of MPFR_FAST_COMPUTE_IF_SMALL_INPUT. But we 147 must also handle the underflow case (an overflow is not possible 148 for small inputs). If an underflow occurred in mpfr_round_near_x, 149 the rounding was to zero or equivalent, and the result is 0, so 150 that the division by 2 will give the wanted result. Otherwise... 151 The rounded result in unbounded exponent range is res/2. If the 152 division by 2 doesn't underflow, it is exact, and we can return 153 this result. And an underflow in the division is a real underflow. 154 In case of directed rounding mode, the result is correct. But in 155 case of rounding to nearest, there is a double rounding problem, 156 and the result is 0 iff the result before the division is the 157 minimum positive number and _inexact has the same sign as z; 158 but in rounding to nearest, res/2 will yield 0 iff |res| is the 159 minimum positive number, so that we just need to test the result 160 of the division and the sign of _inexact. */ 161 mpfr_clear_flags (); 162 MPFR_FAST_COMPUTE_IF_SMALL_INPUT 163 (res, z, -2 * MPFR_GET_EXP (z), 3, 0, r, { 164 int inex2 = mpfr_div_2ui (res, res, 1, r); 165 if (MPFR_UNLIKELY (r == MPFR_RNDN && MPFR_IS_ZERO (res)) && 166 (MPFR_ASSERTN (inex2 != 0), SIGN (_inexact) != MPFR_SIGN (z))) 167 { 168 mpfr_nexttoinf (res); 169 inex = - inex2; 170 } 171 else 172 inex = inex2 != 0 ? inex2 : _inexact; 173 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); 174 goto end; 175 }); 176 } 177 178 /* we can use the asymptotic expansion as soon as |z| > p log(2)/2, 179 but to get some margin we use it for |z| > p/2 */ 180 pbound = MPFR_PREC (res) / 2 + 3; 181 MPFR_ASSERTN (pbound <= ULONG_MAX); 182 MPFR_ALIAS (absz, z, 1, MPFR_EXP (z)); 183 if (mpfr_cmp_ui (absz, pbound) > 0) 184 { 185 inex = mpfr_jn_asympt (res, n, z, r); 186 if (inex != 0) 187 goto end; 188 } 189 190 MPFR_GROUP_INIT_3 (g, 32, y, s, t); 191 192 /* check underflow case: |j(n,z)| <= 1/sqrt(2 Pi n) (ze/2n)^n 193 (see algorithms.tex) */ 194 /* FIXME: the code below doesn't detect all the underflow cases. Either 195 this should be done, or the generic code should detect underflows. */ 196 if (absn > 0) 197 { 198 /* the following is an upper 32-bit approximation to exp(1)/2 */ 199 mpfr_set_str_binary (y, "1.0101101111110000101010001011001"); 200 if (MPFR_SIGN(z) > 0) 201 mpfr_mul (y, y, z, MPFR_RNDU); 202 else 203 { 204 mpfr_mul (y, y, z, MPFR_RNDD); 205 mpfr_neg (y, y, MPFR_RNDU); 206 } 207 mpfr_div_ui (y, y, absn, MPFR_RNDU); 208 /* now y is an upper approximation to |ze/2n|: y < 2^EXP(y), 209 thus |j(n,z)| < 1/2*y^n < 2^(n*EXP(y)-1). 210 If n*EXP(y) < emin then we have an underflow. 211 Note that if emin = MPFR_EMIN_MIN and j = 1, this inequality 212 will never be satisfied. 213 Warning: absn is an unsigned long. */ 214 if ((MPFR_GET_EXP (y) < 0 && absn > - expo.saved_emin) 215 || (absn <= - MPFR_EMIN_MIN && 216 MPFR_GET_EXP (y) < expo.saved_emin / (mpfr_exp_t) absn)) 217 { 218 MPFR_GROUP_CLEAR (g); 219 MPFR_SAVE_EXPO_FREE (expo); 220 return mpfr_underflow (res, (r == MPFR_RNDN) ? MPFR_RNDZ : r, 221 (n % 2) ? ((n > 0) ? MPFR_SIGN(z) : -MPFR_SIGN(z)) 222 : MPFR_SIGN_POS); 223 } 224 } 225 226 /* the logarithm of the ratio between the largest term in the series 227 and the first one is roughly bounded by k0, which we add to the 228 working precision to take into account this cancellation */ 229 /* The following operations avoid integer overflow and ensure that 230 prec <= MPFR_PREC_MAX (prec = MPFR_PREC_MAX won't prevent an abort, 231 but the failure should be handled cleanly). */ 232 k0 = mpfr_jn_k0 (absn, z); 233 MPFR_LOG_MSG (("k0 = %lu\n", k0)); 234 uprec = MPFR_PREC_MAX - 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC_MAX) - 3; 235 if (k0 < uprec) 236 uprec = k0; 237 uprec += MPFR_PREC (res) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (res)) + 3; 238 prec = uprec < MPFR_PREC_MAX ? (mpfr_prec_t) uprec : MPFR_PREC_MAX; 239 240 MPFR_ZIV_INIT (loop, prec); 241 for (;;) 242 { 243 MPFR_BLOCK_DECL (flags); 244 245 MPFR_GROUP_REPREC_3 (g, prec, y, s, t); 246 MPFR_BLOCK (flags, { 247 mpfr_pow_ui (t, z, absn, MPFR_RNDN); /* z^|n| */ 248 mpfr_mul (y, z, z, MPFR_RNDN); /* z^2 */ 249 mpfr_clear_erangeflag (); 250 zz = mpfr_get_ui (y, MPFR_RNDU); 251 /* FIXME: The error analysis is incorrect in case of range error. */ 252 MPFR_ASSERTN (! mpfr_erangeflag_p ()); /* since mpfr_clear_erangeflag */ 253 mpfr_div_2ui (y, y, 2, MPFR_RNDN); /* z^2/4 */ 254 mpfr_fac_ui (s, absn, MPFR_RNDN); /* |n|! */ 255 mpfr_div (t, t, s, MPFR_RNDN); 256 if (absn > 0) 257 mpfr_div_2ui (t, t, absn, MPFR_RNDN); 258 mpfr_set (s, t, MPFR_RNDN); 259 /* note: we assume here that the maximal error bound is proportional to 260 2^exps, which is true also in the case where s=0 */ 261 exps = MPFR_IS_ZERO (s) ? MPFR_EMIN_MIN : MPFR_GET_EXP (s); 262 expT = exps; 263 for (k = 1; ; k++) 264 { 265 MPFR_LOG_MSG (("loop on k, k = %lu\n", k)); 266 mpfr_mul (t, t, y, MPFR_RNDN); 267 mpfr_neg (t, t, MPFR_RNDN); 268 /* Mathematically: absn <= LONG_MAX + 1 <= (ULONG_MAX + 1) / 2, 269 and in practice, k is not very large, so that one should have 270 k + absn <= ULONG_MAX. */ 271 MPFR_ASSERTN (absn <= ULONG_MAX - k); 272 if (k + absn <= ULONG_MAX / k) 273 mpfr_div_ui (t, t, k * (k + absn), MPFR_RNDN); 274 else 275 { 276 mpfr_div_ui (t, t, k, MPFR_RNDN); 277 mpfr_div_ui (t, t, k + absn, MPFR_RNDN); 278 } 279 /* see above note */ 280 exps = MPFR_IS_ZERO (s) ? MPFR_EMIN_MIN : MPFR_GET_EXP (t); 281 if (exps > expT) 282 expT = exps; 283 mpfr_add (s, s, t, MPFR_RNDN); 284 exps = MPFR_IS_ZERO (s) ? MPFR_EMIN_MIN : MPFR_GET_EXP (s); 285 if (exps > expT) 286 expT = exps; 287 /* Above it has been checked that k + absn <= ULONG_MAX. */ 288 if (MPFR_GET_EXP (t) + (mpfr_exp_t) prec <= exps && 289 zz / (2 * k) < k + absn) 290 break; 291 } 292 }); 293 /* the error is bounded by (4k^2+21/2k+7) ulp(s)*2^(expT-exps) 294 <= (k+2)^2 ulp(s)*2^(2+expT-exps) */ 295 diffexp = expT - exps; 296 err = 2 * MPFR_INT_CEIL_LOG2(k + 2) + 2; 297 /* FIXME: Can an overflow occur in the following sum? */ 298 MPFR_ASSERTN (diffexp >= 0 && err >= 0 && 299 diffexp <= MPFR_PREC_MAX - err); 300 err += diffexp; 301 if (MPFR_LIKELY (MPFR_CAN_ROUND (s, prec - err, MPFR_PREC(res), r))) 302 { 303 if (MPFR_LIKELY (! (MPFR_UNDERFLOW (flags) || 304 MPFR_OVERFLOW (flags)))) 305 break; 306 /* The error analysis is incorrect in case of exception. 307 If an underflow or overflow occurred, try once more in 308 a larger precision, and if this happens a second time, 309 then abort to avoid a probable infinite loop. This is 310 a problem that must be fixed! */ 311 MPFR_ASSERTN (! exception); 312 exception = 1; 313 } 314 MPFR_ZIV_NEXT (loop, prec); 315 } 316 MPFR_ZIV_FREE (loop); 317 318 inex = ((n >= 0) || ((n & 1) == 0)) ? mpfr_set (res, s, r) 319 : mpfr_neg (res, s, r); 320 321 MPFR_GROUP_CLEAR (g); 322 323 end: 324 MPFR_SAVE_EXPO_FREE (expo); 325 return mpfr_check_range (res, inex, r); 326} 327 328#define MPFR_JN 329#include "jyn_asympt.c" 330