Deleted Added
sdiff udiff text old ( 97413 ) new ( 141296 )
full compact
e_jn.c (97413) e_jn.c (141296)
1/* @(#)e_jn.c 5.1 93/09/24 */
1
2/* @(#)e_jn.c 1.4 95/01/18 */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
3/*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Developed at SunSoft, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13#ifndef lint
10 * is preserved.
11 * ====================================================
12 */
13
14#ifndef lint
14static char rcsid[] = "$FreeBSD: head/lib/msun/src/e_jn.c 97413 2002-05-28 18:15:04Z alfred $";
15static char rcsid[] = "$FreeBSD: head/lib/msun/src/e_jn.c 141296 2005-02-04 18:26:06Z das $";
15#endif
16
17/*
18 * __ieee754_jn(n, x), __ieee754_yn(n, x)
19 * floating point Bessel's function of the 1st and 2nd kind
20 * of order n
16#endif
17
18/*
19 * __ieee754_jn(n, x), __ieee754_yn(n, x)
20 * floating point Bessel's function of the 1st and 2nd kind
21 * of order n
21 *
22 *
22 * Special cases:
23 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
24 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
25 * Note 2. About jn(n,x), yn(n,x)
26 * For n=0, j0(x) is called,
27 * for n=1, j1(x) is called,
28 * for n<x, forward recursion us used starting
29 * from values of j0(x) and j1(x).
30 * for n>x, a continued fraction approximation to
31 * j(n,x)/j(n-1,x) is evaluated and then backward
32 * recursion is used starting from a supposed value
33 * for j(n,x). The resulting value of j(0,x) is
34 * compared with the actual value to correct the
35 * supposed value of j(n,x).
36 *
37 * yn(n,x) is similar in all respects, except
38 * that forward recursion is used for all
39 * values of n>1.
23 * Special cases:
24 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
25 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
26 * Note 2. About jn(n,x), yn(n,x)
27 * For n=0, j0(x) is called,
28 * for n=1, j1(x) is called,
29 * for n<x, forward recursion us used starting
30 * from values of j0(x) and j1(x).
31 * for n>x, a continued fraction approximation to
32 * j(n,x)/j(n-1,x) is evaluated and then backward
33 * recursion is used starting from a supposed value
34 * for j(n,x). The resulting value of j(0,x) is
35 * compared with the actual value to correct the
36 * supposed value of j(n,x).
37 *
38 * yn(n,x) is similar in all respects, except
39 * that forward recursion is used for all
40 * values of n>1.
40 *
41 *
41 */
42
43#include "math.h"
44#include "math_private.h"
45
46static const double
47invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
48two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */

--- 10 unchanged lines hidden (view full) ---

59
60 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
61 * Thus, J(-n,x) = J(n,-x)
62 */
63 EXTRACT_WORDS(hx,lx,x);
64 ix = 0x7fffffff&hx;
65 /* if J(n,NaN) is NaN */
66 if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
42 */
43
44#include "math.h"
45#include "math_private.h"
46
47static const double
48invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
49two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */

--- 10 unchanged lines hidden (view full) ---

60
61 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
62 * Thus, J(-n,x) = J(n,-x)
63 */
64 EXTRACT_WORDS(hx,lx,x);
65 ix = 0x7fffffff&hx;
66 /* if J(n,NaN) is NaN */
67 if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
67 if(n<0){
68 if(n<0){
68 n = -n;
69 x = -x;
70 hx ^= 0x80000000;
71 }
72 if(n==0) return(__ieee754_j0(x));
73 if(n==1) return(__ieee754_j1(x));
74 sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
75 x = fabs(x);
76 if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */
77 b = zero;
69 n = -n;
70 x = -x;
71 hx ^= 0x80000000;
72 }
73 if(n==0) return(__ieee754_j0(x));
74 if(n==1) return(__ieee754_j1(x));
75 sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
76 x = fabs(x);
77 if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */
78 b = zero;
78 else if((double)n<=x) {
79 else if((double)n<=x) {
79 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
80 if(ix>=0x52D00000) { /* x > 2**302 */
80 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
81 if(ix>=0x52D00000) { /* x > 2**302 */
81 /* (x >> n**2)
82 /* (x >> n**2)
82 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
83 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
83 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
84 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
84 * Let s=sin(x), c=cos(x),
85 * Let s=sin(x), c=cos(x),
85 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
86 *
87 * n sin(xn)*sqt2 cos(xn)*sqt2
88 * ----------------------------------
89 * 0 s-c c+s
90 * 1 -s-c -c+s
91 * 2 -s+c -c-s
92 * 3 s+c c-s
93 */
94 switch(n&3) {
95 case 0: temp = cos(x)+sin(x); break;
96 case 1: temp = -cos(x)+sin(x); break;
97 case 2: temp = -cos(x)-sin(x); break;
98 case 3: temp = cos(x)-sin(x); break;
99 }
100 b = invsqrtpi*temp/sqrt(x);
86 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
87 *
88 * n sin(xn)*sqt2 cos(xn)*sqt2
89 * ----------------------------------
90 * 0 s-c c+s
91 * 1 -s-c -c+s
92 * 2 -s+c -c-s
93 * 3 s+c c-s
94 */
95 switch(n&3) {
96 case 0: temp = cos(x)+sin(x); break;
97 case 1: temp = -cos(x)+sin(x); break;
98 case 2: temp = -cos(x)-sin(x); break;
99 case 3: temp = cos(x)-sin(x); break;
100 }
101 b = invsqrtpi*temp/sqrt(x);
101 } else {
102 } else {
102 a = __ieee754_j0(x);
103 b = __ieee754_j1(x);
104 for(i=1;i<n;i++){
105 temp = b;
106 b = b*((double)(i+i)/x) - a; /* avoid underflow */
107 a = temp;
108 }
109 }
110 } else {
111 if(ix<0x3e100000) { /* x < 2**-29 */
103 a = __ieee754_j0(x);
104 b = __ieee754_j1(x);
105 for(i=1;i<n;i++){
106 temp = b;
107 b = b*((double)(i+i)/x) - a; /* avoid underflow */
108 a = temp;
109 }
110 }
111 } else {
112 if(ix<0x3e100000) { /* x < 2**-29 */
112 /* x is tiny, return the first Taylor expansion of J(n,x)
113 /* x is tiny, return the first Taylor expansion of J(n,x)
113 * J(n,x) = 1/n!*(x/2)^n - ...
114 */
115 if(n>33) /* underflow */
116 b = zero;
117 else {
118 temp = x*0.5; b = temp;
119 for (a=one,i=2;i<=n;i++) {
120 a *= (double)i; /* a = n! */
121 b *= temp; /* b = (x/2)^n */
122 }
123 b = b/a;
124 }
125 } else {
126 /* use backward recurrence */
114 * J(n,x) = 1/n!*(x/2)^n - ...
115 */
116 if(n>33) /* underflow */
117 b = zero;
118 else {
119 temp = x*0.5; b = temp;
120 for (a=one,i=2;i<=n;i++) {
121 a *= (double)i; /* a = n! */
122 b *= temp; /* b = (x/2)^n */
123 }
124 b = b/a;
125 }
126 } else {
127 /* use backward recurrence */
127 /* x x^2 x^2
128 /* x x^2 x^2
128 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
129 * 2n - 2(n+1) - 2(n+2)
130 *
129 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
130 * 2n - 2(n+1) - 2(n+2)
131 *
131 * 1 1 1
132 * 1 1 1
132 * (for large x) = ---- ------ ------ .....
133 * 2n 2(n+1) 2(n+2)
133 * (for large x) = ---- ------ ------ .....
134 * 2n 2(n+1) 2(n+2)
134 * -- - ------ - ------ -
135 * -- - ------ - ------ -
135 * x x x
136 *
137 * Let w = 2n/x and h=2/x, then the above quotient
138 * is equal to the continued fraction:
139 * 1
140 * = -----------------------
141 * 1
142 * w - -----------------
143 * 1
144 * w+h - ---------
145 * w+2h - ...
146 *
147 * To determine how many terms needed, let
148 * Q(0) = w, Q(1) = w(w+h) - 1,
149 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
136 * x x x
137 *
138 * Let w = 2n/x and h=2/x, then the above quotient
139 * is equal to the continued fraction:
140 * 1
141 * = -----------------------
142 * 1
143 * w - -----------------
144 * 1
145 * w+h - ---------
146 * w+2h - ...
147 *
148 * To determine how many terms needed, let
149 * Q(0) = w, Q(1) = w(w+h) - 1,
150 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
150 * When Q(k) > 1e4 good for single
151 * When Q(k) > 1e9 good for double
152 * When Q(k) > 1e17 good for quadruple
151 * When Q(k) > 1e4 good for single
152 * When Q(k) > 1e9 good for double
153 * When Q(k) > 1e17 good for quadruple
153 */
154 /* determine k */
155 double t,v;
156 double q0,q1,h,tmp; int32_t k,m;
157 w = (n+n)/(double)x; h = 2.0/(double)x;
158 q0 = w; z = w+h; q1 = w*z - 1.0; k=1;
159 while(q1<1.0e9) {
160 k += 1; z += h;

--- 62 unchanged lines hidden (view full) ---

223 if(n<0){
224 n = -n;
225 sign = 1 - ((n&1)<<1);
226 }
227 if(n==0) return(__ieee754_y0(x));
228 if(n==1) return(sign*__ieee754_y1(x));
229 if(ix==0x7ff00000) return zero;
230 if(ix>=0x52D00000) { /* x > 2**302 */
154 */
155 /* determine k */
156 double t,v;
157 double q0,q1,h,tmp; int32_t k,m;
158 w = (n+n)/(double)x; h = 2.0/(double)x;
159 q0 = w; z = w+h; q1 = w*z - 1.0; k=1;
160 while(q1<1.0e9) {
161 k += 1; z += h;

--- 62 unchanged lines hidden (view full) ---

224 if(n<0){
225 n = -n;
226 sign = 1 - ((n&1)<<1);
227 }
228 if(n==0) return(__ieee754_y0(x));
229 if(n==1) return(sign*__ieee754_y1(x));
230 if(ix==0x7ff00000) return zero;
231 if(ix>=0x52D00000) { /* x > 2**302 */
231 /* (x >> n**2)
232 /* (x >> n**2)
232 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
233 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
233 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
234 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
234 * Let s=sin(x), c=cos(x),
235 * Let s=sin(x), c=cos(x),
235 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
236 *
237 * n sin(xn)*sqt2 cos(xn)*sqt2
238 * ----------------------------------
239 * 0 s-c c+s
240 * 1 -s-c -c+s
241 * 2 -s+c -c-s
242 * 3 s+c c-s

--- 23 unchanged lines hidden --- 		236 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
237 *
238 * n sin(xn)*sqt2 cos(xn)*sqt2
239 * ----------------------------------
240 * 0 s-c c+s
241 * 1 -s-c -c+s
242 * 2 -s+c -c-s
243 * 3 s+c c-s

--- 23 unchanged lines hidden ---