1 /* mpfr_y0, mpfr_y1, mpfr_yn -- Bessel functions of 2nd kind, integer order.
2 https://pubs.opengroup.org/onlinepubs/9699919799/functions/y0.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 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[%Pd]=%.*Rg rnd=%d", n, mpfr_get_prec (z), mpfr_log_prec, z, r),
158 ("y[%Pd]=%.*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_sqr (t, 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 prec = MPFR_PREC(res) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (res)) + 13;
328
329 mpfr_init2 (y, prec);
330 mpfr_init2 (s1, prec);
331 mpfr_init2 (s2, prec);
332 mpfr_init2 (s3, prec);
333
334 MPFR_ZIV_INIT (loop, prec);
335 for (;;)
336 {
337 mpfr_set_prec (y, prec);
338 mpfr_set_prec (s1, prec);
339 mpfr_set_prec (s2, prec);
340 mpfr_set_prec (s3, prec);
341
342 mpfr_sqr (y, z, MPFR_RNDN);
343 mpfr_div_2ui (y, y, 2, MPFR_RNDN); /* z^2/4 */
344
345 /* store (z/2)^n temporarily in s2 */
346 mpfr_pow_ui (s2, z, absn, MPFR_RNDN);
347 mpfr_div_2si (s2, s2, absn, MPFR_RNDN);
348
349 /* compute S1 * (z/2)^(-n) */
350 if (n == 0)
351 {
352 mpfr_set_ui (s1, 0, MPFR_RNDN);
353 err1 = 0;
354 }
355 else
356 err1 = mpfr_yn_s1 (s1, y, absn - 1);
357 mpfr_div (s1, s1, s2, MPFR_RNDN); /* (z/2)^(-n) * S1 */
358 /* See algorithms.tex: the relative error on s1 is bounded by
359 (3n+3)*2^(e+1-prec). */
360 err1 = MPFR_INT_CEIL_LOG2 (3 * absn + 3) + err1 + 1;
361 /* rel_err(s1) <= 2^(err1-prec), thus err(s1) <= 2^err1 ulps */
362
363 /* compute (z/2)^n * S3 */
364 mpfr_neg (y, y, MPFR_RNDN); /* -z^2/4 */
365 err3 = mpfr_yn_s3 (s3, y, s2, absn); /* (z/2)^n * S3 */
366 /* the error on s3 is bounded by 2^err3 ulps */
367
368 /* add s1+s3 */
369 err1 += MPFR_EXP(s1);
370 mpfr_add (s1, s1, s3, MPFR_RNDN);
371 /* the error is bounded by 1/2 + 2^err1*2^(- EXP(s1))
372 + 2^err3*2^(EXP(s3) - EXP(s1)) */
373 err3 += MPFR_EXP(s3);
374 err1 = (err3 > err1) ? err3 + 1 : err1 + 1;
375 err1 -= MPFR_EXP(s1);
376 err1 = (err1 >= 0) ? err1 + 1 : 1;
377 /* now the error on s1 is bounded by 2^err1*ulp(s1) */
378
379 /* compute S2 */
380 mpfr_div_2ui (s2, z, 1, MPFR_RNDN); /* z/2 */
381 mpfr_log (s2, s2, MPFR_RNDN); /* log(z/2) */
382 mpfr_const_euler (s3, MPFR_RNDN);
383 err2 = MPFR_EXP(s2) > MPFR_EXP(s3) ? MPFR_EXP(s2) : MPFR_EXP(s3);
384 mpfr_add (s2, s2, s3, MPFR_RNDN); /* log(z/2) + gamma */
385 err2 -= MPFR_EXP(s2);
386 mpfr_mul_2ui (s2, s2, 1, MPFR_RNDN); /* 2*(log(z/2) + gamma) */
387 mpfr_jn (s3, absn, z, MPFR_RNDN); /* Jn(z) */
388 mpfr_mul (s2, s2, s3, MPFR_RNDN); /* 2*(log(z/2) + gamma)*Jn(z) */
389 err2 += 4; /* the error on s2 is bounded by 2^err2 ulps, see
390 algorithms.tex */
391
392 /* add all three sums */
393 err1 += MPFR_EXP(s1); /* the error on s1 is bounded by 2^err1 */
394 err2 += MPFR_EXP(s2); /* the error on s2 is bounded by 2^err2 */
395 mpfr_sub (s2, s2, s1, MPFR_RNDN); /* s2 - (s1+s3) */
396 err2 = (err1 > err2) ? err1 + 1 : err2 + 1;
397 err2 -= MPFR_EXP(s2);
398 err2 = (err2 >= 0) ? err2 + 1 : 1;
399 /* now the error on s2 is bounded by 2^err2*ulp(s2) */
400 mpfr_const_pi (y, MPFR_RNDN); /* error bounded by 1 ulp */
401 mpfr_div (s2, s2, y, MPFR_RNDN); /* error bounded by
402 2^(err2+1)*ulp(s2) */
403 err2 ++;
404
405 if (MPFR_LIKELY (MPFR_CAN_ROUND (s2, prec - err2, MPFR_PREC(res), r)))
406 break;
407 MPFR_ZIV_NEXT (loop, prec);
408 }
409 MPFR_ZIV_FREE (loop);
410
411 /* Assume two's complement for the test n & 1 */
412 inex = mpfr_set4 (res, s2, r, n >= 0 || (n & 1) == 0 ?
413 MPFR_SIGN (s2) : - MPFR_SIGN (s2));
414
415 mpfr_clear (y);
416 mpfr_clear (s1);
417 mpfr_clear (s2);
418 mpfr_clear (s3);
419 }
420
421 end:
422 MPFR_SAVE_EXPO_FREE (expo);
423 return mpfr_check_range (res, inex, r);
424 }
425
426 #define MPFR_YN
427 #include "jyn_asympt.c"