1 /* mpfr_j0, mpfr_j1, mpfr_jn -- Bessel functions of 1st kind, integer order.
2 https://pubs.opengroup.org/onlinepubs/9699919799/functions/j0.html
3
4 Copyright 2007-2023 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramba 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 https://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 /* Relations: j(-n,z) = (-1)^n j(n,z)
28 j(n,-z) = (-1)^n j(n,z)
29 */
30
31 static int mpfr_jn_asympt (mpfr_ptr, long, mpfr_srcptr, mpfr_rnd_t);
32
33 int
34 mpfr_j0 (mpfr_ptr res, mpfr_srcptr z, mpfr_rnd_t r)
35 {
36 return mpfr_jn (res, 0, z, r);
37 }
38
39 int
40 mpfr_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 */
49 static unsigned long
50 mpfr_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
82 int
83 mpfr_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[%Pd]=%.*Rg rnd=%d", n, mpfr_get_prec (z), mpfr_log_prec, z, r),
99 ("res[%Pd]=%.*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), VSIGN (_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_IS_POS (z))
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_sqr (y, 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 /* the expected number of lost bits is k0, if err is larger than k0
315 most probably there is a cancellation in the series, thus we add
316 err - k0 bits to prec */
317 if (err > k0)
318 MPFR_INC_PREC (prec, err - k0);
319 MPFR_ZIV_NEXT (loop, prec);
320 }
321 MPFR_ZIV_FREE (loop);
322
323 inex = ((n >= 0) || ((n & 1) == 0)) ? mpfr_set (res, s, r)
324 : mpfr_neg (res, s, r);
325
326 MPFR_GROUP_CLEAR (g);
327
328 end:
329 MPFR_SAVE_EXPO_FREE (expo);
330 return mpfr_check_range (res, inex, r);
331 }
332
333 #define MPFR_JN
334 #include "jyn_asympt.c"