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 __jnl = jnl
31ddc0e0b5SRichard Lowe #pragma weak __ynl = ynl
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 "longdouble.h"
6025c28e83SPiotr Jasiukajtis #include <float.h> /* LDBL_MAX */
6125c28e83SPiotr Jasiukajtis
6225c28e83SPiotr Jasiukajtis #define GENERIC long double
6325c28e83SPiotr Jasiukajtis
6425c28e83SPiotr Jasiukajtis static const GENERIC
6525c28e83SPiotr Jasiukajtis invsqrtpi = 5.641895835477562869480794515607725858441e-0001L,
6625c28e83SPiotr Jasiukajtis two = 2.0L,
6725c28e83SPiotr Jasiukajtis zero = 0.0L,
6825c28e83SPiotr Jasiukajtis one = 1.0L;
6925c28e83SPiotr Jasiukajtis
7025c28e83SPiotr Jasiukajtis GENERIC
jnl(int n,GENERIC x)71*685c1a21SRichard Lowe jnl(int n, GENERIC x)
72*685c1a21SRichard Lowe {
7325c28e83SPiotr Jasiukajtis int i, sgn;
7425c28e83SPiotr Jasiukajtis GENERIC a, b, temp = 0, z, w;
7525c28e83SPiotr Jasiukajtis
7625c28e83SPiotr Jasiukajtis /*
7725c28e83SPiotr Jasiukajtis * J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
7825c28e83SPiotr Jasiukajtis * Thus, J(-n,x) = J(n,-x)
7925c28e83SPiotr Jasiukajtis */
8025c28e83SPiotr Jasiukajtis if (n < 0) {
8125c28e83SPiotr Jasiukajtis n = -n;
8225c28e83SPiotr Jasiukajtis x = -x;
8325c28e83SPiotr Jasiukajtis }
84*685c1a21SRichard Lowe if (n == 0)
85*685c1a21SRichard Lowe return (j0l(x));
86*685c1a21SRichard Lowe if (n == 1)
87*685c1a21SRichard Lowe return (j1l(x));
88*685c1a21SRichard Lowe if (x != x)
89*685c1a21SRichard Lowe return (x+x);
9025c28e83SPiotr Jasiukajtis if ((n&1) == 0)
91*685c1a21SRichard Lowe sgn = 0; /* even n */
9225c28e83SPiotr Jasiukajtis else
9325c28e83SPiotr Jasiukajtis sgn = signbitl(x); /* old n */
9425c28e83SPiotr Jasiukajtis x = fabsl(x);
9525c28e83SPiotr Jasiukajtis if (x == zero || !finitel(x)) b = zero;
9625c28e83SPiotr Jasiukajtis else if ((GENERIC)n <= x) {
9725c28e83SPiotr Jasiukajtis /*
9825c28e83SPiotr Jasiukajtis * Safe to use
9925c28e83SPiotr Jasiukajtis * J(n+1,x)=2n/x *J(n,x)-J(n-1,x)
10025c28e83SPiotr Jasiukajtis */
101*685c1a21SRichard Lowe if (x > 1.0e91L) {
10225c28e83SPiotr Jasiukajtis /*
10325c28e83SPiotr Jasiukajtis * x >> n**2
10425c28e83SPiotr Jasiukajtis * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
10525c28e83SPiotr Jasiukajtis * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
10625c28e83SPiotr Jasiukajtis * Let s=sin(x), c=cos(x),
10725c28e83SPiotr Jasiukajtis * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
10825c28e83SPiotr Jasiukajtis *
10925c28e83SPiotr Jasiukajtis * n sin(xn)*sqt2 cos(xn)*sqt2
11025c28e83SPiotr Jasiukajtis * ----------------------------------
11125c28e83SPiotr Jasiukajtis * 0 s-c c+s
112*685c1a21SRichard Lowe * 1 -s-c -c+s
11325c28e83SPiotr Jasiukajtis * 2 -s+c -c-s
11425c28e83SPiotr Jasiukajtis * 3 s+c c-s
11525c28e83SPiotr Jasiukajtis */
116*685c1a21SRichard Lowe switch (n&3) {
117*685c1a21SRichard Lowe case 0:
118*685c1a21SRichard Lowe temp = cosl(x)+sinl(x);
119*685c1a21SRichard Lowe break;
120*685c1a21SRichard Lowe case 1:
121*685c1a21SRichard Lowe temp = -cosl(x)+sinl(x);
122*685c1a21SRichard Lowe break;
123*685c1a21SRichard Lowe case 2:
124*685c1a21SRichard Lowe temp = -cosl(x)-sinl(x);
125*685c1a21SRichard Lowe break;
126*685c1a21SRichard Lowe case 3:
127*685c1a21SRichard Lowe temp = cosl(x)-sinl(x);
128*685c1a21SRichard Lowe break;
129*685c1a21SRichard Lowe }
130*685c1a21SRichard Lowe b = invsqrtpi*temp/sqrtl(x);
131*685c1a21SRichard Lowe } else {
13225c28e83SPiotr Jasiukajtis a = j0l(x);
13325c28e83SPiotr Jasiukajtis b = j1l(x);
13425c28e83SPiotr Jasiukajtis for (i = 1; i < n; i++) {
135*685c1a21SRichard Lowe temp = b;
136*685c1a21SRichard Lowe /* avoid underflow */
137*685c1a21SRichard Lowe b = b*((GENERIC)(i+i)/x) - a;
138*685c1a21SRichard Lowe a = temp;
13925c28e83SPiotr Jasiukajtis }
14025c28e83SPiotr Jasiukajtis }
141*685c1a21SRichard Lowe } else {
142*685c1a21SRichard Lowe if (x < 1e-17L) { /* use J(n,x) = 1/n!*(x/2)^n */
143*685c1a21SRichard Lowe b = powl(0.5L*x, (GENERIC)n);
144*685c1a21SRichard Lowe if (b != zero) {
145*685c1a21SRichard Lowe for (a = one, i = 1; i <= n; i++)
146*685c1a21SRichard Lowe a *= (GENERIC)i;
147*685c1a21SRichard Lowe b = b/a;
148*685c1a21SRichard Lowe }
149*685c1a21SRichard Lowe } else {
150*685c1a21SRichard Lowe /* BEGIN CSTYLED */
151*685c1a21SRichard Lowe /*
152*685c1a21SRichard Lowe * use backward recurrence
153*685c1a21SRichard Lowe * x x^2 x^2
154*685c1a21SRichard Lowe * J(n,x)/J(n-1,x) = ---- ------ ------ .....
155*685c1a21SRichard Lowe * 2n - 2(n+1) - 2(n+2)
156*685c1a21SRichard Lowe *
157*685c1a21SRichard Lowe * 1 1 1
158*685c1a21SRichard Lowe * (for large x) = ---- ------ ------ .....
159*685c1a21SRichard Lowe * 2n 2(n+1) 2(n+2)
160*685c1a21SRichard Lowe * -- - ------ - ------ -
161*685c1a21SRichard Lowe * x x x
162*685c1a21SRichard Lowe *
163*685c1a21SRichard Lowe * Let w = 2n/x and h=2/x, then the above quotient
164*685c1a21SRichard Lowe * is equal to the continued fraction:
165*685c1a21SRichard Lowe * 1
166*685c1a21SRichard Lowe * = -----------------------
167*685c1a21SRichard Lowe * 1
168*685c1a21SRichard Lowe * w - -----------------
169*685c1a21SRichard Lowe * 1
170*685c1a21SRichard Lowe * w+h - ---------
171*685c1a21SRichard Lowe * w+2h - ...
172*685c1a21SRichard Lowe *
173*685c1a21SRichard Lowe * To determine how many terms needed, let
174*685c1a21SRichard Lowe * Q(0) = w, Q(1) = w(w+h) - 1,
175*685c1a21SRichard Lowe * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
176*685c1a21SRichard Lowe * When Q(k) > 1e4 good for single
177*685c1a21SRichard Lowe * When Q(k) > 1e9 good for double
178*685c1a21SRichard Lowe * When Q(k) > 1e17 good for quaduple
179*685c1a21SRichard Lowe */
180*685c1a21SRichard Lowe /* END CSTYLED */
181*685c1a21SRichard Lowe /* determine k */
182*685c1a21SRichard Lowe GENERIC t, v;
183*685c1a21SRichard Lowe double q0, q1, h, tmp;
184*685c1a21SRichard Lowe int k, m;
185*685c1a21SRichard Lowe w = (n+n)/(double)x;
186*685c1a21SRichard Lowe h = 2.0/(double)x;
187*685c1a21SRichard Lowe q0 = w;
188*685c1a21SRichard Lowe z = w+h;
189*685c1a21SRichard Lowe q1 = w*z - 1.0;
190*685c1a21SRichard Lowe k = 1;
191*685c1a21SRichard Lowe while (q1 < 1.0e17) {
192*685c1a21SRichard Lowe k += 1;
193*685c1a21SRichard Lowe z += h;
194*685c1a21SRichard Lowe tmp = z*q1 - q0;
195*685c1a21SRichard Lowe q0 = q1;
196*685c1a21SRichard Lowe q1 = tmp;
197*685c1a21SRichard Lowe }
198*685c1a21SRichard Lowe m = n+n;
199*685c1a21SRichard Lowe for (t = zero, i = 2*(n+k); i >= m; i -= 2)
200*685c1a21SRichard Lowe t = one/(i/x-t);
201*685c1a21SRichard Lowe a = t;
202*685c1a21SRichard Lowe b = one;
20325c28e83SPiotr Jasiukajtis /*
20425c28e83SPiotr Jasiukajtis * Estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
20525c28e83SPiotr Jasiukajtis * hence, if n*(log(2n/x)) > ...
206*685c1a21SRichard Lowe * single:
207*685c1a21SRichard Lowe * 8.8722839355e+01
208*685c1a21SRichard Lowe * double:
209*685c1a21SRichard Lowe * 7.09782712893383973096e+02
210*685c1a21SRichard Lowe * long double:
211*685c1a21SRichard Lowe * 1.1356523406294143949491931077970765006170e+04
21225c28e83SPiotr Jasiukajtis * then recurrent value may overflow and the result is
21325c28e83SPiotr Jasiukajtis * likely underflow to zero.
21425c28e83SPiotr Jasiukajtis */
215*685c1a21SRichard Lowe tmp = n;
216*685c1a21SRichard Lowe v = two/x;
217*685c1a21SRichard Lowe tmp = tmp*logl(fabsl(v*tmp));
218*685c1a21SRichard Lowe if (tmp < 1.1356523406294143949491931077970765e+04L) {
21925c28e83SPiotr Jasiukajtis for (i = n-1; i > 0; i--) {
220*685c1a21SRichard Lowe temp = b;
221*685c1a21SRichard Lowe b = ((i+i)/x)*b - a;
222*685c1a21SRichard Lowe a = temp;
22325c28e83SPiotr Jasiukajtis }
224*685c1a21SRichard Lowe } else {
22525c28e83SPiotr Jasiukajtis for (i = n-1; i > 0; i--) {
226*685c1a21SRichard Lowe temp = b;
227*685c1a21SRichard Lowe b = ((i+i)/x)*b - a;
228*685c1a21SRichard Lowe a = temp;
229*685c1a21SRichard Lowe if (b > 1e1000L) {
23025c28e83SPiotr Jasiukajtis a /= b;
23125c28e83SPiotr Jasiukajtis t /= b;
23225c28e83SPiotr Jasiukajtis b = 1.0;
23325c28e83SPiotr Jasiukajtis }
23425c28e83SPiotr Jasiukajtis }
235*685c1a21SRichard Lowe }
23625c28e83SPiotr Jasiukajtis b = (t*j0l(x)/b);
237*685c1a21SRichard Lowe }
23825c28e83SPiotr Jasiukajtis }
239*685c1a21SRichard Lowe if (sgn != 0)
240*685c1a21SRichard Lowe return (-b);
24125c28e83SPiotr Jasiukajtis else
242*685c1a21SRichard Lowe return (b);
24325c28e83SPiotr Jasiukajtis }
24425c28e83SPiotr Jasiukajtis
24525c28e83SPiotr Jasiukajtis GENERIC
ynl(int n,GENERIC x)246*685c1a21SRichard Lowe ynl(int n, GENERIC x)
247*685c1a21SRichard Lowe {
24825c28e83SPiotr Jasiukajtis int i;
24925c28e83SPiotr Jasiukajtis int sign;
25025c28e83SPiotr Jasiukajtis GENERIC a, b, temp = 0;
25125c28e83SPiotr Jasiukajtis
25225c28e83SPiotr Jasiukajtis if (x != x)
253*685c1a21SRichard Lowe return (x+x);
25425c28e83SPiotr Jasiukajtis if (x <= zero) {
25525c28e83SPiotr Jasiukajtis if (x == zero)
256*685c1a21SRichard Lowe return (-one/zero);
25725c28e83SPiotr Jasiukajtis else
258*685c1a21SRichard Lowe return (zero/zero);
25925c28e83SPiotr Jasiukajtis }
26025c28e83SPiotr Jasiukajtis sign = 1;
26125c28e83SPiotr Jasiukajtis if (n < 0) {
26225c28e83SPiotr Jasiukajtis n = -n;
263*685c1a21SRichard Lowe if ((n&1) == 1)
264*685c1a21SRichard Lowe sign = -1;
26525c28e83SPiotr Jasiukajtis }
266*685c1a21SRichard Lowe if (n == 0)
267*685c1a21SRichard Lowe return (y0l(x));
268*685c1a21SRichard Lowe if (n == 1)
269*685c1a21SRichard Lowe return (sign*y1l(x));
270*685c1a21SRichard Lowe if (!finitel(x))
271*685c1a21SRichard Lowe return (zero);
27225c28e83SPiotr Jasiukajtis
27325c28e83SPiotr Jasiukajtis if (x > 1.0e91L) {
27425c28e83SPiotr Jasiukajtis /*
27525c28e83SPiotr Jasiukajtis * x >> n**2
27625c28e83SPiotr Jasiukajtis * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
27725c28e83SPiotr Jasiukajtis * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
27825c28e83SPiotr Jasiukajtis * Let s=sin(x), c=cos(x),
27925c28e83SPiotr Jasiukajtis * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
28025c28e83SPiotr Jasiukajtis *
28125c28e83SPiotr Jasiukajtis * n sin(xn)*sqt2 cos(xn)*sqt2
28225c28e83SPiotr Jasiukajtis * ----------------------------------
283*685c1a21SRichard Lowe * 0 s-c c+s
284*685c1a21SRichard Lowe * 1 -s-c -c+s
285*685c1a21SRichard Lowe * 2 -s+c -c-s
28625c28e83SPiotr Jasiukajtis * 3 s+c c-s
28725c28e83SPiotr Jasiukajtis */
28825c28e83SPiotr Jasiukajtis switch (n&3) {
289*685c1a21SRichard Lowe case 0:
290*685c1a21SRichard Lowe temp = sinl(x)-cosl(x);
291*685c1a21SRichard Lowe break;
292*685c1a21SRichard Lowe case 1:
293*685c1a21SRichard Lowe temp = -sinl(x)-cosl(x);
294*685c1a21SRichard Lowe break;
295*685c1a21SRichard Lowe case 2:
296*685c1a21SRichard Lowe temp = -sinl(x)+cosl(x);
297*685c1a21SRichard Lowe break;
298*685c1a21SRichard Lowe case 3:
299*685c1a21SRichard Lowe temp = sinl(x)+cosl(x);
300*685c1a21SRichard Lowe break;
30125c28e83SPiotr Jasiukajtis }
30225c28e83SPiotr Jasiukajtis b = invsqrtpi*temp/sqrtl(x);
30325c28e83SPiotr Jasiukajtis } else {
30425c28e83SPiotr Jasiukajtis a = y0l(x);
30525c28e83SPiotr Jasiukajtis b = y1l(x);
30625c28e83SPiotr Jasiukajtis /*
30725c28e83SPiotr Jasiukajtis * fix 1262058 and take care of non-default rounding
30825c28e83SPiotr Jasiukajtis */
30925c28e83SPiotr Jasiukajtis for (i = 1; i < n; i++) {
31025c28e83SPiotr Jasiukajtis temp = b;
31125c28e83SPiotr Jasiukajtis b *= (GENERIC) (i + i) / x;
31225c28e83SPiotr Jasiukajtis if (b <= -LDBL_MAX)
31325c28e83SPiotr Jasiukajtis break;
31425c28e83SPiotr Jasiukajtis b -= a;
31525c28e83SPiotr Jasiukajtis a = temp;
31625c28e83SPiotr Jasiukajtis }
31725c28e83SPiotr Jasiukajtis }
31825c28e83SPiotr Jasiukajtis if (sign > 0)
319*685c1a21SRichard Lowe return (b);
32025c28e83SPiotr Jasiukajtis else
321*685c1a21SRichard Lowe return (-b);
32225c28e83SPiotr Jasiukajtis }
323