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