1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1998-2012 Ross Ihaka and the R Core team.
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, a copy is available at
17 * http://www.r-project.org/Licenses/
18 */
19
20 /* DESCRIPTION --> see below */
21
22
23 /* From http://www.netlib.org/specfun/rybesl Fortran translated by f2c,...
24 * ------------------------------=#---- Martin Maechler, ETH Zurich
25 */
26 #include "nmath.h"
27 #include "bessel.h"
28
29 #ifndef MATHLIB_STANDALONE
30 #include <R_ext/Memory.h>
31 #endif
32
33 #define min0(x, y) (((x) <= (y)) ? (x) : (y))
34
35 static void Y_bessel(double *x, double *alpha, long *nb,
36 double *by, long *ncalc);
37
bessel_y(double x,double alpha)38 double bessel_y(double x, double alpha)
39 {
40 long nb, ncalc;
41 double na, *by;
42 #ifndef MATHLIB_STANDALONE
43 const void *vmax;
44 #endif
45
46 #ifdef IEEE_754
47 /* NaNs propagated correctly */
48 if (ISNAN(x) || ISNAN(alpha)) return x + alpha;
49 #endif
50 if (x < 0) {
51 ML_ERROR(ME_RANGE, "bessel_y");
52 return ML_NAN;
53 }
54 na = floor(alpha);
55 if (alpha < 0) {
56 /* Using Abramowitz & Stegun 9.1.2
57 * this may not be quite optimal (CPU and accuracy wise) */
58 return(bessel_y(x, -alpha) * cos(M_PI * alpha) -
59 ((alpha == na) ? 0 :
60 bessel_j(x, -alpha) * sin(M_PI * alpha)));
61 }
62 nb = 1+ (long)na;/* nb-1 <= alpha < nb */
63 alpha -= (double)(nb-1);
64 #ifdef MATHLIB_STANDALONE
65 by = (double *) calloc(nb, sizeof(double));
66 if (!by) MATHLIB_ERROR("%s", _("bessel_y allocation error"));
67 #else
68 vmax = vmaxget();
69 by = (double *) R_alloc((size_t) nb, sizeof(double));
70 #endif
71 Y_bessel(&x, &alpha, &nb, by, &ncalc);
72 if(ncalc != nb) {/* error input */
73 if(ncalc == -1)
74 return ML_POSINF;
75 else if(ncalc < -1)
76 MATHLIB_WARNING4(_("bessel_y(%g): ncalc (=%ld) != nb (=%ld); alpha=%g. Arg. out of range?\n"),
77 x, ncalc, nb, alpha);
78 else /* ncalc >= 0 */
79 MATHLIB_WARNING2(_("bessel_y(%g,nu=%g): precision lost in result\n"),
80 x, alpha+(double)nb-1);
81 }
82 x = by[nb-1];
83 #ifdef MATHLIB_STANDALONE
84 free(by);
85 #else
86 vmaxset(vmax);
87 #endif
88 return x;
89 }
90
91 /* modified version of bessel_y that accepts a work array instead of
92 allocating one. */
bessel_y_ex(double x,double alpha,double * by)93 double bessel_y_ex(double x, double alpha, double *by)
94 {
95 long nb, ncalc;
96 double na;
97
98 #ifdef IEEE_754
99 /* NaNs propagated correctly */
100 if (ISNAN(x) || ISNAN(alpha)) return x + alpha;
101 #endif
102 if (x < 0) {
103 ML_ERROR(ME_RANGE, "bessel_y");
104 return ML_NAN;
105 }
106 na = floor(alpha);
107 if (alpha < 0) {
108 /* Using Abramowitz & Stegun 9.1.2
109 * this may not be quite optimal (CPU and accuracy wise) */
110 return(bessel_y_ex(x, -alpha, by) * cos(M_PI * alpha) -
111 ((alpha == na) ? 0 :
112 bessel_j_ex(x, -alpha, by) * sin(M_PI * alpha)));
113 }
114 nb = 1+ (long)na;/* nb-1 <= alpha < nb */
115 alpha -= (double)(nb-1);
116 Y_bessel(&x, &alpha, &nb, by, &ncalc);
117 if(ncalc != nb) {/* error input */
118 if(ncalc == -1)
119 return ML_POSINF;
120 else if(ncalc < -1)
121 MATHLIB_WARNING4(_("bessel_y(%g): ncalc (=%ld) != nb (=%ld); alpha=%g. Arg. out of range?\n"),
122 x, ncalc, nb, alpha);
123 else /* ncalc >= 0 */
124 MATHLIB_WARNING2(_("bessel_y(%g,nu=%g): precision lost in result\n"),
125 x, alpha+(double)nb-1);
126 }
127 x = by[nb-1];
128 return x;
129 }
130
Y_bessel(double * x,double * alpha,long * nb,double * by,long * ncalc)131 static void Y_bessel(double *x, double *alpha, long *nb,
132 double *by, long *ncalc)
133 {
134 /* ----------------------------------------------------------------------
135
136 This routine calculates Bessel functions Y_(N+ALPHA) (X)
137 for non-negative argument X, and non-negative order N+ALPHA.
138
139
140 Explanation of variables in the calling sequence
141
142 X - Non-negative argument for which
143 Y's are to be calculated.
144 ALPHA - Fractional part of order for which
145 Y's are to be calculated. 0 <= ALPHA < 1.0.
146 NB - Number of functions to be calculated, NB > 0.
147 The first function calculated is of order ALPHA, and the
148 last is of order (NB - 1 + ALPHA).
149 BY - Output vector of length NB. If the
150 routine terminates normally (NCALC=NB), the vector BY
151 contains the functions Y(ALPHA,X), ... , Y(NB-1+ALPHA,X),
152 If (0 < NCALC < NB), BY(I) contains correct function
153 values for I <= NCALC, and contains the ratios
154 Y(ALPHA+I-1,X)/Y(ALPHA+I-2,X) for the rest of the array.
155 NCALC - Output variable indicating possible errors.
156 Before using the vector BY, the user should check that
157 NCALC=NB, i.e., all orders have been calculated to
158 the desired accuracy. See error returns below.
159
160
161 *******************************************************************
162
163 Error returns
164
165 In case of an error, NCALC != NB, and not all Y's are
166 calculated to the desired accuracy.
167
168 NCALC < -1: An argument is out of range. For example,
169 NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >=
170 XMAX. In this case, BY[0] = 0.0, the remainder of the
171 BY-vector is not calculated, and NCALC is set to
172 MIN0(NB,0)-2 so that NCALC != NB.
173 NCALC = -1: Y(ALPHA,X) >= XINF. The requested function
174 values are set to 0.0.
175 1 < NCALC < NB: Not all requested function values could
176 be calculated accurately. BY(I) contains correct function
177 values for I <= NCALC, and and the remaining NB-NCALC
178 array elements contain 0.0.
179
180
181 Intrinsic functions required are:
182
183 DBLE, EXP, INT, MAX, MIN, REAL, SQRT
184
185
186 Acknowledgement
187
188 This program draws heavily on Temme's Algol program for Y(a,x)
189 and Y(a+1,x) and on Campbell's programs for Y_nu(x). Temme's
190 scheme is used for x < THRESH, and Campbell's scheme is used
191 in the asymptotic region. Segments of code from both sources
192 have been translated into Fortran 77, merged, and heavily modified.
193 Modifications include parameterization of machine dependencies,
194 use of a new approximation for ln(gamma(x)), and built-in
195 protection against over/underflow.
196
197 References: "Bessel functions J_nu(x) and Y_nu(x) of float
198 order and float argument," Campbell, J. B.,
199 Comp. Phy. Comm. 18, 1979, pp. 133-142.
200
201 "On the numerical evaluation of the ordinary
202 Bessel function of the second kind," Temme,
203 N. M., J. Comput. Phys. 21, 1976, pp. 343-350.
204
205 Latest modification: March 19, 1990
206
207 Modified by: W. J. Cody
208 Applied Mathematics Division
209 Argonne National Laboratory
210 Argonne, IL 60439
211 ----------------------------------------------------------------------*/
212
213
214 /* ----------------------------------------------------------------------
215 Mathematical constants
216 FIVPI = 5*PI
217 PIM5 = 5*PI - 15
218 ----------------------------------------------------------------------*/
219 const static double fivpi = 15.707963267948966192;
220 const static double pim5 = .70796326794896619231;
221
222 /*----------------------------------------------------------------------
223 Coefficients for Chebyshev polynomial expansion of
224 1/gamma(1-x), abs(x) <= .5
225 ----------------------------------------------------------------------*/
226 const static double ch[21] = { -6.7735241822398840964e-24,
227 -6.1455180116049879894e-23,2.9017595056104745456e-21,
228 1.3639417919073099464e-19,2.3826220476859635824e-18,
229 -9.0642907957550702534e-18,-1.4943667065169001769e-15,
230 -3.3919078305362211264e-14,-1.7023776642512729175e-13,
231 9.1609750938768647911e-12,2.4230957900482704055e-10,
232 1.7451364971382984243e-9,-3.3126119768180852711e-8,
233 -8.6592079961391259661e-7,-4.9717367041957398581e-6,
234 7.6309597585908126618e-5,.0012719271366545622927,
235 .0017063050710955562222,-.07685284084478667369,
236 -.28387654227602353814,.92187029365045265648 };
237
238 /* Local variables */
239 long i, k, na;
240
241 double alfa, div, ddiv, even, gamma, term, cosmu, sinmu,
242 b, c, d, e, f, g, h, p, q, r, s, d1, d2, q0, pa,pa1, qa,qa1,
243 en, en1, nu, ex, ya,ya1, twobyx, den, odd, aye, dmu, x2, xna;
244
245 en1 = ya = ya1 = 0; /* -Wall */
246
247 ex = *x;
248 nu = *alpha;
249 if (*nb > 0 && 0. <= nu && nu < 1.) {
250 if(ex < DBL_MIN || ex > xlrg_BESS_Y) {
251 /* Warning is not really appropriate, give
252 * proper limit:
253 * ML_ERROR(ME_RANGE, "Y_bessel"); */
254 *ncalc = *nb;
255 if(ex > xlrg_BESS_Y) by[0]= 0.; /*was ML_POSINF */
256 else if(ex < DBL_MIN) by[0]=ML_NEGINF;
257 for(i=0; i < *nb; i++)
258 by[i] = by[0];
259 return;
260 }
261 xna = ftrunc(nu + .5);
262 na = (long) xna;
263 if (na == 1) {/* <==> .5 <= *alpha < 1 <==> -5. <= nu < 0 */
264 nu -= xna;
265 }
266 if (nu == -.5) {
267 p = M_SQRT_2dPI / sqrt(ex);
268 ya = p * sin(ex);
269 ya1 = -p * cos(ex);
270 } else if (ex < 3.) {
271 /* -------------------------------------------------------------
272 Use Temme's scheme for small X
273 ------------------------------------------------------------- */
274 b = ex * .5;
275 d = -log(b);
276 f = nu * d;
277 e = pow(b, -nu);
278 if (fabs(nu) < M_eps_sinc)
279 c = M_1_PI;
280 else
281 c = nu / sin(nu * M_PI);
282
283 /* ------------------------------------------------------------
284 Computation of sinh(f)/f
285 ------------------------------------------------------------ */
286 if (fabs(f) < 1.) {
287 x2 = f * f;
288 en = 19.;
289 s = 1.;
290 for (i = 1; i <= 9; ++i) {
291 s = s * x2 / en / (en - 1.) + 1.;
292 en -= 2.;
293 }
294 } else {
295 s = (e - 1. / e) * .5 / f;
296 }
297 /* --------------------------------------------------------
298 Computation of 1/gamma(1-a) using Chebyshev polynomials */
299 x2 = nu * nu * 8.;
300 aye = ch[0];
301 even = 0.;
302 alfa = ch[1];
303 odd = 0.;
304 for (i = 3; i <= 19; i += 2) {
305 even = -(aye + aye + even);
306 aye = -even * x2 - aye + ch[i - 1];
307 odd = -(alfa + alfa + odd);
308 alfa = -odd * x2 - alfa + ch[i];
309 }
310 even = (even * .5 + aye) * x2 - aye + ch[20];
311 odd = (odd + alfa) * 2.;
312 gamma = odd * nu + even;
313 /* End of computation of 1/gamma(1-a)
314 ----------------------------------------------------------- */
315 g = e * gamma;
316 e = (e + 1. / e) * .5;
317 f = 2. * c * (odd * e + even * s * d);
318 e = nu * nu;
319 p = g * c;
320 q = M_1_PI / g;
321 c = nu * M_PI_2;
322 if (fabs(c) < M_eps_sinc)
323 r = 1.;
324 else
325 r = sin(c) / c;
326
327 r = M_PI * c * r * r;
328 c = 1.;
329 d = -b * b;
330 h = 0.;
331 ya = f + r * q;
332 ya1 = p;
333 en = 1.;
334
335 while (fabs(g / (1. + fabs(ya))) +
336 fabs(h / (1. + fabs(ya1))) > DBL_EPSILON) {
337 f = (f * en + p + q) / (en * en - e);
338 c *= (d / en);
339 p /= en - nu;
340 q /= en + nu;
341 g = c * (f + r * q);
342 h = c * p - en * g;
343 ya += g;
344 ya1+= h;
345 en += 1.;
346 }
347 ya = -ya;
348 ya1 = -ya1 / b;
349 } else if (ex < thresh_BESS_Y) {
350 /* --------------------------------------------------------------
351 Use Temme's scheme for moderate X : 3 <= x < 16
352 -------------------------------------------------------------- */
353 c = (.5 - nu) * (.5 + nu);
354 b = ex + ex;
355 e = ex * M_1_PI * cos(nu * M_PI) / DBL_EPSILON;
356 e *= e;
357 p = 1.;
358 q = -ex;
359 r = 1. + ex * ex;
360 s = r;
361 en = 2.;
362 while (r * en * en < e) {
363 en1 = en + 1.;
364 d = (en - 1. + c / en) / s;
365 p = (en + en - p * d) / en1;
366 q = (-b + q * d) / en1;
367 s = p * p + q * q;
368 r *= s;
369 en = en1;
370 }
371 f = p / s;
372 p = f;
373 g = -q / s;
374 q = g;
375 L220:
376 en -= 1.;
377 if (en > 0.) {
378 r = en1 * (2. - p) - 2.;
379 s = b + en1 * q;
380 d = (en - 1. + c / en) / (r * r + s * s);
381 p = d * r;
382 q = d * s;
383 e = f + 1.;
384 f = p * e - g * q;
385 g = q * e + p * g;
386 en1 = en;
387 goto L220;
388 }
389 f = 1. + f;
390 d = f * f + g * g;
391 pa = f / d;
392 qa = -g / d;
393 d = nu + .5 - p;
394 q += ex;
395 pa1 = (pa * q - qa * d) / ex;
396 qa1 = (qa * q + pa * d) / ex;
397 b = ex - M_PI_2 * (nu + .5);
398 c = cos(b);
399 s = sin(b);
400 d = M_SQRT_2dPI / sqrt(ex);
401 ya = d * (pa * s + qa * c);
402 ya1 = d * (qa1 * s - pa1 * c);
403 } else { /* x > thresh_BESS_Y */
404 /* ----------------------------------------------------------
405 Use Campbell's asymptotic scheme.
406 ---------------------------------------------------------- */
407 na = 0;
408 d1 = ftrunc(ex / fivpi);
409 i = (long) d1;
410 dmu = ex - 15. * d1 - d1 * pim5 - (*alpha + .5) * M_PI_2;
411 if (i - (i / 2 << 1) == 0) {
412 cosmu = cos(dmu);
413 sinmu = sin(dmu);
414 } else {
415 cosmu = -cos(dmu);
416 sinmu = -sin(dmu);
417 }
418 ddiv = 8. * ex;
419 dmu = *alpha;
420 den = sqrt(ex);
421 for (k = 1; k <= 2; ++k) {
422 p = cosmu;
423 cosmu = sinmu;
424 sinmu = -p;
425 d1 = (2. * dmu - 1.) * (2. * dmu + 1.);
426 d2 = 0.;
427 div = ddiv;
428 p = 0.;
429 q = 0.;
430 q0 = d1 / div;
431 term = q0;
432 for (i = 2; i <= 20; ++i) {
433 d2 += 8.;
434 d1 -= d2;
435 div += ddiv;
436 term = -term * d1 / div;
437 p += term;
438 d2 += 8.;
439 d1 -= d2;
440 div += ddiv;
441 term *= (d1 / div);
442 q += term;
443 if (fabs(term) <= DBL_EPSILON) {
444 break;
445 }
446 }
447 p += 1.;
448 q += q0;
449 if (k == 1)
450 ya = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den;
451 else
452 ya1 = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den;
453 dmu += 1.;
454 }
455 }
456 if (na == 1) {
457 h = 2. * (nu + 1.) / ex;
458 if (h > 1.) {
459 if (fabs(ya1) > DBL_MAX / h) {
460 h = 0.;
461 ya = 0.;
462 }
463 }
464 h = h * ya1 - ya;
465 ya = ya1;
466 ya1 = h;
467 }
468
469 /* ---------------------------------------------------------------
470 Now have first one or two Y's
471 --------------------------------------------------------------- */
472 by[0] = ya;
473 *ncalc = 1;
474 if(*nb > 1) {
475 by[1] = ya1;
476 if (ya1 != 0.) {
477 aye = 1. + *alpha;
478 twobyx = 2. / ex;
479 *ncalc = 2;
480 for (i = 2; i < *nb; ++i) {
481 if (twobyx < 1.) {
482 if (fabs(by[i - 1]) * twobyx >= DBL_MAX / aye)
483 goto L450;
484 } else {
485 if (fabs(by[i - 1]) >= DBL_MAX / aye / twobyx)
486 goto L450;
487 }
488 by[i] = twobyx * aye * by[i - 1] - by[i - 2];
489 aye += 1.;
490 ++(*ncalc);
491 }
492 }
493 }
494 L450:
495 for (i = *ncalc; i < *nb; ++i)
496 by[i] = ML_NEGINF;/* was 0 */
497
498 } else {
499 by[0] = 0.;
500 *ncalc = min0(*nb,0) - 1;
501 }
502 }
503
504