4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
30 #pragma weak __jnl = jnl
31 #pragma weak __ynl = ynl
34 * floating point Bessel's function of the 1st and 2nd kind
35 * of order n: jn(n,x),yn(n,x);
38 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
39 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
40 * Note 2. About jn(n,x), yn(n,x)
41 * For n=0, j0(x) is called,
42 * for n=1, j1(x) is called,
43 * for n<x, forward recursion us used starting
44 * from values of j0(x) and j1(x).
45 * for n>x, a continued fraction approximation to
46 * j(n,x)/j(n-1,x) is evaluated and then backward
47 * recursion is used starting from a supposed value
48 * for j(n,x). The resulting value of j(0,x) is
49 * compared with the actual value to correct the
50 * supposed value of j(n,x).
52 * yn(n,x) is similar in all respects, except
53 * that forward recursion is used for all
59 #include "longdouble.h"
60 #include <float.h> /* LDBL_MAX */
62 #define GENERIC long double
65 invsqrtpi
= 5.641895835477562869480794515607725858441e-0001L,
71 jnl(n
, x
) int n
; GENERIC x
; {
73 GENERIC a
, b
, temp
= 0, z
, w
;
76 * J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
77 * Thus, J(-n,x) = J(n,-x)
83 if (n
== 0) return (j0l(x
));
84 if (n
== 1) return (j1l(x
));
85 if (x
!= x
) return x
+x
;
89 sgn
= signbitl(x
); /* old n */
91 if (x
== zero
|| !finitel(x
)) b
= zero
;
92 else if ((GENERIC
)n
<= x
) {
95 * J(n+1,x)=2n/x *J(n,x)-J(n-1,x)
100 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
101 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
102 * Let s=sin(x), c=cos(x),
103 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
105 * n sin(xn)*sqt2 cos(xn)*sqt2
106 * ----------------------------------
113 case 0: temp
= cosl(x
)+sinl(x
); break;
114 case 1: temp
= -cosl(x
)+sinl(x
); break;
115 case 2: temp
= -cosl(x
)-sinl(x
); break;
116 case 3: temp
= cosl(x
)-sinl(x
); break;
118 b
= invsqrtpi
*temp
/sqrtl(x
);
122 for (i
= 1; i
< n
; i
++) {
124 b
= b
*((GENERIC
)(i
+i
)/x
) - a
; /* avoid underflow */
129 if (x
< 1e-17L) { /* use J(n,x) = 1/n!*(x/2)^n */
130 b
= powl(0.5L*x
, (GENERIC
) n
);
132 for (a
= one
, i
= 1; i
<= n
; i
++) a
*= (GENERIC
)i
;
137 * use backward recurrence
139 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
140 * 2n - 2(n+1) - 2(n+2)
143 * (for large x) = ---- ------ ------ .....
145 * -- - ------ - ------ -
148 * Let w = 2n/x and h=2/x, then the above quotient
149 * is equal to the continued fraction:
151 * = -----------------------
153 * w - -----------------
158 * To determine how many terms needed, let
159 * Q(0) = w, Q(1) = w(w+h) - 1,
160 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
161 * When Q(k) > 1e4 good for single
162 * When Q(k) > 1e9 good for double
163 * When Q(k) > 1e17 good for quaduple
167 double q0
, q1
, h
, tmp
; int k
, m
;
168 w
= (n
+n
)/(double)x
; h
= 2.0/(double)x
;
169 q0
= w
; z
= w
+h
; q1
= w
*z
- 1.0; k
= 1;
170 while (q1
< 1.0e17
) {
177 for (t
= zero
, i
= 2*(n
+k
); i
>= m
; i
-= 2) t
= one
/(i
/x
-t
);
181 * Estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
182 * hence, if n*(log(2n/x)) > ...
183 * single 8.8722839355e+01
184 * double 7.09782712893383973096e+02
185 * long double 1.1356523406294143949491931077970765006170e+04
186 * then recurrent value may overflow and the result is
187 * likely underflow to zero.
191 tmp
= tmp
*logl(fabsl(v
*tmp
));
192 if (tmp
< 1.1356523406294143949491931077970765e+04L) {
193 for (i
= n
-1; i
> 0; i
--) {
199 for (i
= n
-1; i
> 0; i
--) {
220 ynl(n
, x
) int n
; GENERIC x
; {
223 GENERIC a
, b
, temp
= 0;
236 if ((n
&1) == 1) sign
= -1;
238 if (n
== 0) return (y0l(x
));
239 if (n
== 1) return (sign
*y1l(x
));
240 if (!finitel(x
)) return zero
;
245 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
246 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
247 * Let s=sin(x), c=cos(x),
248 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
250 * n sin(xn)*sqt2 cos(xn)*sqt2
251 * ----------------------------------
258 case 0: temp
= sinl(x
)-cosl(x
); break;
259 case 1: temp
= -sinl(x
)-cosl(x
); break;
260 case 2: temp
= -sinl(x
)+cosl(x
); break;
261 case 3: temp
= sinl(x
)+cosl(x
); break;
263 b
= invsqrtpi
*temp
/sqrtl(x
);
268 * fix 1262058 and take care of non-default rounding
270 for (i
= 1; i
< n
; i
++) {
272 b
*= (GENERIC
) (i
+ i
) / x
;