1 /* mpfr_j0, mpfr_j1, mpfr_jn -- Bessel functions of 1st kind, integer order.
2 http://www.opengroup.org/onlinepubs/009695399/functions/j0.html
4 Copyright 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
5 Contributed by the Arenaire and Caramel projects, INRIA.
7 This file is part of the GNU MPFR Library.
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.
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.
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. */
24 #define MPFR_NEED_LONGLONG_H
25 #include "mpfr-impl.h"
27 /* Relations: j(-n,z) = (-1)^n j(n,z)
28 j(n,-z) = (-1)^n j(n,z)
31 static int mpfr_jn_asympt (mpfr_ptr, long, mpfr_srcptr, mpfr_rnd_t);
34 mpfr_j0 (mpfr_ptr res, mpfr_srcptr z, mpfr_rnd_t r)
36 return mpfr_jn (res, 0, z, r);
40 mpfr_j1 (mpfr_ptr res, mpfr_srcptr z, mpfr_rnd_t r)
42 return mpfr_jn (res, 1, z, r);
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).
50 mpfr_jn_k0 (unsigned long n, mpfr_srcptr z)
59 mpfr_abs (t, z, MPFR_RNDN); /* t = 2*k1 */
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 */
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);
83 mpfr_jn (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r)
88 mpfr_prec_t prec, pbound, err;
90 mpfr_exp_t exps, expT, diffexp;
92 unsigned long k, zz, k0;
94 MPFR_SAVE_EXPO_DECL (expo);
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));
102 absn = SAFE_ABS (unsigned long, n);
104 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (z)))
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
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 */
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);
125 return mpfr_set_ui (res, 0, r);
129 MPFR_SAVE_EXPO_MARK (expo);
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. */
134 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (res, __gmpfr_one, -2 * MPFR_GET_EXP (z),
135 2, 0, r, inex = _inexact; goto end);
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. */
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. */
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)))
168 mpfr_nexttoinf (res);
172 inex = inex2 != 0 ? inex2 : _inexact;
173 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
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)
185 inex = mpfr_jn_asympt (res, n, z, r);
190 MPFR_GROUP_INIT_3 (g, 32, y, s, t);
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. */
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);
204 mpfr_mul (y, y, z, MPFR_RNDD);
205 mpfr_neg (y, y, MPFR_RNDU);
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))
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))
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;
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;
240 MPFR_ZIV_INIT (loop, prec);
243 MPFR_BLOCK_DECL (flags);
245 MPFR_GROUP_REPREC_3 (g, prec, y, s, t);
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);
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);
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);
276 mpfr_div_ui (t, t, k, MPFR_RNDN);
277 mpfr_div_ui (t, t, k + absn, MPFR_RNDN);
280 exps = MPFR_IS_ZERO (s) ? MPFR_EMIN_MIN : MPFR_GET_EXP (t);
283 mpfr_add (s, s, t, MPFR_RNDN);
284 exps = MPFR_IS_ZERO (s) ? MPFR_EMIN_MIN : MPFR_GET_EXP (s);
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)
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);
301 if (MPFR_LIKELY (MPFR_CAN_ROUND (s, prec - err, MPFR_PREC(res), r)))
303 if (MPFR_LIKELY (! (MPFR_UNDERFLOW (flags) ||
304 MPFR_OVERFLOW (flags))))
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);
314 MPFR_ZIV_NEXT (loop, prec);
316 MPFR_ZIV_FREE (loop);
318 inex = ((n >= 0) || ((n & 1) == 0)) ? mpfr_set (res, s, r)
319 : mpfr_neg (res, s, r);
321 MPFR_GROUP_CLEAR (g);
324 MPFR_SAVE_EXPO_FREE (expo);
325 return mpfr_check_range (res, inex, r);
329 #include "jyn_asympt.c"