125c28e83SPiotr Jasiukajtis /*
225c28e83SPiotr Jasiukajtis * CDDL HEADER START
325c28e83SPiotr Jasiukajtis *
425c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e83SPiotr Jasiukajtis *
825c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e83SPiotr Jasiukajtis * and limitations under the License.
1225c28e83SPiotr Jasiukajtis *
1325c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e83SPiotr Jasiukajtis *
1925c28e83SPiotr Jasiukajtis * CDDL HEADER END
2025c28e83SPiotr Jasiukajtis */
2125c28e83SPiotr Jasiukajtis
2225c28e83SPiotr Jasiukajtis /*
2325c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
2425c28e83SPiotr Jasiukajtis */
2525c28e83SPiotr Jasiukajtis /*
2625c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
2725c28e83SPiotr Jasiukajtis * Use is subject to license terms.
2825c28e83SPiotr Jasiukajtis */
2925c28e83SPiotr Jasiukajtis
30ddc0e0b5SRichard Lowe #pragma weak __jn = jn
31ddc0e0b5SRichard Lowe #pragma weak __yn = yn
3225c28e83SPiotr Jasiukajtis
3325c28e83SPiotr Jasiukajtis /*
3425c28e83SPiotr Jasiukajtis * floating point Bessel's function of the 1st and 2nd kind
3525c28e83SPiotr Jasiukajtis * of order n: jn(n,x),yn(n,x);
3625c28e83SPiotr Jasiukajtis *
3725c28e83SPiotr Jasiukajtis * Special cases:
3825c28e83SPiotr Jasiukajtis * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
3925c28e83SPiotr Jasiukajtis * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
4025c28e83SPiotr Jasiukajtis * Note 2. About jn(n,x), yn(n,x)
4125c28e83SPiotr Jasiukajtis * For n=0, j0(x) is called,
4225c28e83SPiotr Jasiukajtis * for n=1, j1(x) is called,
4325c28e83SPiotr Jasiukajtis * for n<x, forward recursion us used starting
4425c28e83SPiotr Jasiukajtis * from values of j0(x) and j1(x).
4525c28e83SPiotr Jasiukajtis * for n>x, a continued fraction approximation to
4625c28e83SPiotr Jasiukajtis * j(n,x)/j(n-1,x) is evaluated and then backward
4725c28e83SPiotr Jasiukajtis * recursion is used starting from a supposed value
4825c28e83SPiotr Jasiukajtis * for j(n,x). The resulting value of j(0,x) is
4925c28e83SPiotr Jasiukajtis * compared with the actual value to correct the
5025c28e83SPiotr Jasiukajtis * supposed value of j(n,x).
5125c28e83SPiotr Jasiukajtis *
5225c28e83SPiotr Jasiukajtis * yn(n,x) is similar in all respects, except
5325c28e83SPiotr Jasiukajtis * that forward recursion is used for all
5425c28e83SPiotr Jasiukajtis * values of n>1.
5525c28e83SPiotr Jasiukajtis *
5625c28e83SPiotr Jasiukajtis */
5725c28e83SPiotr Jasiukajtis
5825c28e83SPiotr Jasiukajtis #include "libm.h"
5925c28e83SPiotr Jasiukajtis #include <float.h> /* DBL_MIN */
6025c28e83SPiotr Jasiukajtis #include <values.h> /* X_TLOSS */
6125c28e83SPiotr Jasiukajtis #include "xpg6.h" /* __xpg6 */
6225c28e83SPiotr Jasiukajtis
6325c28e83SPiotr Jasiukajtis #define GENERIC double
6425c28e83SPiotr Jasiukajtis
6525c28e83SPiotr Jasiukajtis static const GENERIC
6625c28e83SPiotr Jasiukajtis invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
6725c28e83SPiotr Jasiukajtis two = 2.0,
6825c28e83SPiotr Jasiukajtis zero = 0.0,
6925c28e83SPiotr Jasiukajtis one = 1.0;
7025c28e83SPiotr Jasiukajtis
7125c28e83SPiotr Jasiukajtis GENERIC
jn(int n,GENERIC x)72*685c1a21SRichard Lowe jn(int n, GENERIC x)
73*685c1a21SRichard Lowe {
7425c28e83SPiotr Jasiukajtis int i, sgn;
7525c28e83SPiotr Jasiukajtis GENERIC a, b, temp = 0;
7625c28e83SPiotr Jasiukajtis GENERIC z, w, ox, on;
7725c28e83SPiotr Jasiukajtis
7825c28e83SPiotr Jasiukajtis /*
7925c28e83SPiotr Jasiukajtis * J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
8025c28e83SPiotr Jasiukajtis * Thus, J(-n,x) = J(n,-x)
8125c28e83SPiotr Jasiukajtis */
82*685c1a21SRichard Lowe ox = x;
83*685c1a21SRichard Lowe on = (GENERIC)n;
84*685c1a21SRichard Lowe
8525c28e83SPiotr Jasiukajtis if (n < 0) {
8625c28e83SPiotr Jasiukajtis n = -n;
8725c28e83SPiotr Jasiukajtis x = -x;
8825c28e83SPiotr Jasiukajtis }
8925c28e83SPiotr Jasiukajtis if (isnan(x))
9025c28e83SPiotr Jasiukajtis return (x*x); /* + -> * for Cheetah */
91*685c1a21SRichard Lowe if (!((int)_lib_version == libm_ieee ||
92*685c1a21SRichard Lowe (__xpg6 & _C99SUSv3_math_errexcept) != 0)) {
93*685c1a21SRichard Lowe if (fabs(x) > X_TLOSS)
9425c28e83SPiotr Jasiukajtis return (_SVID_libm_err(on, ox, 38));
9525c28e83SPiotr Jasiukajtis }
9625c28e83SPiotr Jasiukajtis if (n == 0)
9725c28e83SPiotr Jasiukajtis return (j0(x));
9825c28e83SPiotr Jasiukajtis if (n == 1)
9925c28e83SPiotr Jasiukajtis return (j1(x));
10025c28e83SPiotr Jasiukajtis if ((n&1) == 0)
101*685c1a21SRichard Lowe sgn = 0; /* even n */
10225c28e83SPiotr Jasiukajtis else
10325c28e83SPiotr Jasiukajtis sgn = signbit(x); /* old n */
10425c28e83SPiotr Jasiukajtis x = fabs(x);
10525c28e83SPiotr Jasiukajtis if (x == zero||!finite(x)) b = zero;
10625c28e83SPiotr Jasiukajtis else if ((GENERIC)n <= x) {
10725c28e83SPiotr Jasiukajtis /*
10825c28e83SPiotr Jasiukajtis * Safe to use
10925c28e83SPiotr Jasiukajtis * J(n+1,x)=2n/x *J(n,x)-J(n-1,x)
11025c28e83SPiotr Jasiukajtis */
111*685c1a21SRichard Lowe if (x > 1.0e91) {
11225c28e83SPiotr Jasiukajtis /*
11325c28e83SPiotr Jasiukajtis * x >> n**2
11425c28e83SPiotr Jasiukajtis * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
11525c28e83SPiotr Jasiukajtis * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
11625c28e83SPiotr Jasiukajtis * Let s=sin(x), c=cos(x),
11725c28e83SPiotr Jasiukajtis * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
11825c28e83SPiotr Jasiukajtis *
11925c28e83SPiotr Jasiukajtis * n sin(xn)*sqt2 cos(xn)*sqt2
12025c28e83SPiotr Jasiukajtis * ----------------------------------
12125c28e83SPiotr Jasiukajtis * 0 s-c c+s
122*685c1a21SRichard Lowe * 1 -s-c -c+s
12325c28e83SPiotr Jasiukajtis * 2 -s+c -c-s
12425c28e83SPiotr Jasiukajtis * 3 s+c c-s
12525c28e83SPiotr Jasiukajtis */
126*685c1a21SRichard Lowe switch (n&3) {
127*685c1a21SRichard Lowe case 0:
128*685c1a21SRichard Lowe temp = cos(x)+sin(x);
129*685c1a21SRichard Lowe break;
130*685c1a21SRichard Lowe case 1:
131*685c1a21SRichard Lowe temp = -cos(x)+sin(x);
132*685c1a21SRichard Lowe break;
133*685c1a21SRichard Lowe case 2:
134*685c1a21SRichard Lowe temp = -cos(x)-sin(x);
135*685c1a21SRichard Lowe break;
136*685c1a21SRichard Lowe case 3:
137*685c1a21SRichard Lowe temp = cos(x)-sin(x);
138*685c1a21SRichard Lowe break;
139*685c1a21SRichard Lowe }
140*685c1a21SRichard Lowe b = invsqrtpi*temp/sqrt(x);
141*685c1a21SRichard Lowe } else {
14225c28e83SPiotr Jasiukajtis a = j0(x);
14325c28e83SPiotr Jasiukajtis b = j1(x);
14425c28e83SPiotr Jasiukajtis for (i = 1; i < n; i++) {
145*685c1a21SRichard Lowe temp = b;
146*685c1a21SRichard Lowe /* avoid underflow */
147*685c1a21SRichard Lowe b = b*((GENERIC)(i+i)/x) - a;
148*685c1a21SRichard Lowe a = temp;
14925c28e83SPiotr Jasiukajtis }
15025c28e83SPiotr Jasiukajtis }
151*685c1a21SRichard Lowe } else {
152*685c1a21SRichard Lowe if (x < 1e-9) { /* use J(n,x) = 1/n!*(x/2)^n */
153*685c1a21SRichard Lowe b = pow(0.5*x, (GENERIC) n);
154*685c1a21SRichard Lowe if (b != zero) {
155*685c1a21SRichard Lowe for (a = one, i = 1; i <= n; i++)
156*685c1a21SRichard Lowe a *= (GENERIC)i;
157*685c1a21SRichard Lowe b = b/a;
158*685c1a21SRichard Lowe }
15925c28e83SPiotr Jasiukajtis } else {
160*685c1a21SRichard Lowe /*
161*685c1a21SRichard Lowe * use backward recurrence
162*685c1a21SRichard Lowe * x x^2 x^2
163*685c1a21SRichard Lowe * J(n,x)/J(n-1,x) = ---- ------ ------ .....
164*685c1a21SRichard Lowe * 2n - 2(n+1) - 2(n+2)
165*685c1a21SRichard Lowe *
166*685c1a21SRichard Lowe * 1 1 1
167*685c1a21SRichard Lowe * (for large x) = ---- ------ ------ .....
168*685c1a21SRichard Lowe * 2n 2(n+1) 2(n+2)
169*685c1a21SRichard Lowe * -- - ------ - ------ -
170*685c1a21SRichard Lowe * x x x
171*685c1a21SRichard Lowe *
172*685c1a21SRichard Lowe * Let w = 2n/x and h = 2/x, then the above quotient
173*685c1a21SRichard Lowe * is equal to the continued fraction:
174*685c1a21SRichard Lowe * 1
175*685c1a21SRichard Lowe * = -----------------------
176*685c1a21SRichard Lowe * 1
177*685c1a21SRichard Lowe * w - -----------------
178*685c1a21SRichard Lowe * 1
179*685c1a21SRichard Lowe * w+h - ---------
180*685c1a21SRichard Lowe * w+2h - ...
181*685c1a21SRichard Lowe *
182*685c1a21SRichard Lowe * To determine how many terms needed, let
183*685c1a21SRichard Lowe * Q(0) = w, Q(1) = w(w+h) - 1,
184*685c1a21SRichard Lowe * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
185*685c1a21SRichard Lowe * When Q(k) > 1e4 good for single
186*685c1a21SRichard Lowe * When Q(k) > 1e9 good for double
187*685c1a21SRichard Lowe * When Q(k) > 1e17 good for quaduple
188*685c1a21SRichard Lowe */
189*685c1a21SRichard Lowe /* determine k */
190*685c1a21SRichard Lowe GENERIC t, v;
191*685c1a21SRichard Lowe double q0, q1, h, tmp;
192*685c1a21SRichard Lowe int k, m;
193*685c1a21SRichard Lowe w = (n+n)/(double)x;
194*685c1a21SRichard Lowe h = 2.0/(double)x;
195*685c1a21SRichard Lowe q0 = w;
196*685c1a21SRichard Lowe z = w + h;
197*685c1a21SRichard Lowe q1 = w*z - 1.0;
198*685c1a21SRichard Lowe k = 1;
199*685c1a21SRichard Lowe
200*685c1a21SRichard Lowe while (q1 < 1.0e9) {
201*685c1a21SRichard Lowe k += 1;
202*685c1a21SRichard Lowe z += h;
203*685c1a21SRichard Lowe tmp = z*q1 - q0;
204*685c1a21SRichard Lowe q0 = q1;
205*685c1a21SRichard Lowe q1 = tmp;
206*685c1a21SRichard Lowe }
207*685c1a21SRichard Lowe m = n+n;
208*685c1a21SRichard Lowe for (t = zero, i = 2*(n+k); i >= m; i -= 2)
209*685c1a21SRichard Lowe t = one/(i/x-t);
210*685c1a21SRichard Lowe a = t;
211*685c1a21SRichard Lowe b = one;
212*685c1a21SRichard Lowe /*
213*685c1a21SRichard Lowe * estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
214*685c1a21SRichard Lowe * hence, if n*(log(2n/x)) > ...
215*685c1a21SRichard Lowe * single:
216*685c1a21SRichard Lowe * 8.8722839355e+01
217*685c1a21SRichard Lowe * double:
218*685c1a21SRichard Lowe * 7.09782712893383973096e+02
219*685c1a21SRichard Lowe * long double:
220*685c1a21SRichard Lowe * 1.1356523406294143949491931077970765006170e+04
221*685c1a21SRichard Lowe * then recurrent value may overflow and the result is
222*685c1a21SRichard Lowe * likely underflow to zero
223*685c1a21SRichard Lowe */
224*685c1a21SRichard Lowe tmp = n;
225*685c1a21SRichard Lowe v = two/x;
226*685c1a21SRichard Lowe tmp = tmp*log(fabs(v*tmp));
227*685c1a21SRichard Lowe if (tmp < 7.09782712893383973096e+02) {
228*685c1a21SRichard Lowe for (i = n-1; i > 0; i--) {
229*685c1a21SRichard Lowe temp = b;
230*685c1a21SRichard Lowe b = ((i+i)/x)*b - a;
231*685c1a21SRichard Lowe a = temp;
232*685c1a21SRichard Lowe }
233*685c1a21SRichard Lowe } else {
23425c28e83SPiotr Jasiukajtis for (i = n-1; i > 0; i--) {
235*685c1a21SRichard Lowe temp = b;
236*685c1a21SRichard Lowe b = ((i+i)/x)*b - a;
237*685c1a21SRichard Lowe a = temp;
23825c28e83SPiotr Jasiukajtis if (b > 1e100) {
23925c28e83SPiotr Jasiukajtis a /= b;
24025c28e83SPiotr Jasiukajtis t /= b;
24125c28e83SPiotr Jasiukajtis b = 1.0;
24225c28e83SPiotr Jasiukajtis }
24325c28e83SPiotr Jasiukajtis }
244*685c1a21SRichard Lowe }
24525c28e83SPiotr Jasiukajtis b = (t*j0(x)/b);
246*685c1a21SRichard Lowe }
24725c28e83SPiotr Jasiukajtis }
248*685c1a21SRichard Lowe if (sgn != 0)
24925c28e83SPiotr Jasiukajtis return (-b);
25025c28e83SPiotr Jasiukajtis else
25125c28e83SPiotr Jasiukajtis return (b);
25225c28e83SPiotr Jasiukajtis }
25325c28e83SPiotr Jasiukajtis
25425c28e83SPiotr Jasiukajtis GENERIC
yn(int n,GENERIC x)255*685c1a21SRichard Lowe yn(int n, GENERIC x)
256*685c1a21SRichard Lowe {
25725c28e83SPiotr Jasiukajtis int i;
25825c28e83SPiotr Jasiukajtis int sign;
25925c28e83SPiotr Jasiukajtis GENERIC a, b, temp = 0, ox, on;
26025c28e83SPiotr Jasiukajtis
261*685c1a21SRichard Lowe ox = x;
262*685c1a21SRichard Lowe on = (GENERIC)n;
26325c28e83SPiotr Jasiukajtis if (isnan(x))
26425c28e83SPiotr Jasiukajtis return (x*x); /* + -> * for Cheetah */
26525c28e83SPiotr Jasiukajtis if (x <= zero) {
26625c28e83SPiotr Jasiukajtis if (x == zero) {
26725c28e83SPiotr Jasiukajtis /* return -one/zero; */
26825c28e83SPiotr Jasiukajtis return (_SVID_libm_err((GENERIC)n, x, 12));
26925c28e83SPiotr Jasiukajtis } else {
27025c28e83SPiotr Jasiukajtis /* return zero/zero; */
27125c28e83SPiotr Jasiukajtis return (_SVID_libm_err((GENERIC)n, x, 13));
27225c28e83SPiotr Jasiukajtis }
27325c28e83SPiotr Jasiukajtis }
274*685c1a21SRichard Lowe if (!((int)_lib_version == libm_ieee ||
275*685c1a21SRichard Lowe (__xpg6 & _C99SUSv3_math_errexcept) != 0)) {
276*685c1a21SRichard Lowe if (x > X_TLOSS)
27725c28e83SPiotr Jasiukajtis return (_SVID_libm_err(on, ox, 39));
27825c28e83SPiotr Jasiukajtis }
27925c28e83SPiotr Jasiukajtis sign = 1;
28025c28e83SPiotr Jasiukajtis if (n < 0) {
28125c28e83SPiotr Jasiukajtis n = -n;
28225c28e83SPiotr Jasiukajtis if ((n&1) == 1) sign = -1;
28325c28e83SPiotr Jasiukajtis }
28425c28e83SPiotr Jasiukajtis if (n == 0)
28525c28e83SPiotr Jasiukajtis return (y0(x));
28625c28e83SPiotr Jasiukajtis if (n == 1)
28725c28e83SPiotr Jasiukajtis return (sign*y1(x));
28825c28e83SPiotr Jasiukajtis if (!finite(x))
28925c28e83SPiotr Jasiukajtis return (zero);
29025c28e83SPiotr Jasiukajtis
29125c28e83SPiotr Jasiukajtis if (x > 1.0e91) {
29225c28e83SPiotr Jasiukajtis /*
29325c28e83SPiotr Jasiukajtis * x >> n**2
29425c28e83SPiotr Jasiukajtis * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
29525c28e83SPiotr Jasiukajtis * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
29625c28e83SPiotr Jasiukajtis * Let s = sin(x), c = cos(x),
29725c28e83SPiotr Jasiukajtis * xn = x-(2n+1)*pi/4, sqt2 = sqrt(2), then
29825c28e83SPiotr Jasiukajtis *
29925c28e83SPiotr Jasiukajtis * n sin(xn)*sqt2 cos(xn)*sqt2
30025c28e83SPiotr Jasiukajtis * ----------------------------------
30125c28e83SPiotr Jasiukajtis * 0 s-c c+s
302*685c1a21SRichard Lowe * 1 -s-c -c+s
30325c28e83SPiotr Jasiukajtis * 2 -s+c -c-s
30425c28e83SPiotr Jasiukajtis * 3 s+c c-s
30525c28e83SPiotr Jasiukajtis */
30625c28e83SPiotr Jasiukajtis switch (n&3) {
307*685c1a21SRichard Lowe case 0:
308*685c1a21SRichard Lowe temp = sin(x)-cos(x);
309*685c1a21SRichard Lowe break;
310*685c1a21SRichard Lowe case 1:
311*685c1a21SRichard Lowe temp = -sin(x)-cos(x);
312*685c1a21SRichard Lowe break;
313*685c1a21SRichard Lowe case 2:
314*685c1a21SRichard Lowe temp = -sin(x)+cos(x);
315*685c1a21SRichard Lowe break;
316*685c1a21SRichard Lowe case 3:
317*685c1a21SRichard Lowe temp = sin(x)+cos(x);
318*685c1a21SRichard Lowe break;
31925c28e83SPiotr Jasiukajtis }
32025c28e83SPiotr Jasiukajtis b = invsqrtpi*temp/sqrt(x);
32125c28e83SPiotr Jasiukajtis } else {
32225c28e83SPiotr Jasiukajtis a = y0(x);
32325c28e83SPiotr Jasiukajtis b = y1(x);
32425c28e83SPiotr Jasiukajtis /*
32525c28e83SPiotr Jasiukajtis * fix 1262058 and take care of non-default rounding
32625c28e83SPiotr Jasiukajtis */
32725c28e83SPiotr Jasiukajtis for (i = 1; i < n; i++) {
32825c28e83SPiotr Jasiukajtis temp = b;
32925c28e83SPiotr Jasiukajtis b *= (GENERIC) (i + i) / x;
33025c28e83SPiotr Jasiukajtis if (b <= -DBL_MAX)
33125c28e83SPiotr Jasiukajtis break;
33225c28e83SPiotr Jasiukajtis b -= a;
33325c28e83SPiotr Jasiukajtis a = temp;
33425c28e83SPiotr Jasiukajtis }
33525c28e83SPiotr Jasiukajtis }
33625c28e83SPiotr Jasiukajtis if (sign > 0)
33725c28e83SPiotr Jasiukajtis return (b);
33825c28e83SPiotr Jasiukajtis else
33925c28e83SPiotr Jasiukajtis return (-b);
34025c28e83SPiotr Jasiukajtis }
341