55 * that forward recursion is used for all
56 * values of n>1.
57 *
58 */
59
60 #include "libm.h"
61 #include "longdouble.h"
62 #include <float.h> /* LDBL_MAX */
63
64 #define GENERIC long double
65
66 static const GENERIC
67 invsqrtpi= 5.641895835477562869480794515607725858441e-0001L,
68 two = 2.0L,
69 zero = 0.0L,
70 one = 1.0L;
71
72 GENERIC
73 jnl(n,x) int n; GENERIC x;{
74 int i, sgn;
75 GENERIC a, b, temp, z, w;
76
77 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
78 * Thus, J(-n,x) = J(n,-x)
79 */
80 if(n<0){
81 n = -n;
82 x = -x;
83 }
84 if(n==0) return(j0l(x));
85 if(n==1) return(j1l(x));
86 if(x!=x) return x+x;
87 if((n&1)==0)
88 sgn=0; /* even n */
89 else
90 sgn = signbitl(x); /* old n */
91 x = fabsl(x);
92 if(x == zero||!finitel(x)) b = zero;
93 else if((GENERIC)n<=x) { /* Safe to use
94 J(n+1,x)=2n/x *J(n,x)-J(n-1,x)
95 */
195 temp = b;
196 b = ((i+i)/x)*b - a;
197 a = temp;
198 if(b>1e1000L) {
199 a /= b;
200 t /= b;
201 b = 1.0;
202 }
203 }
204 }
205 b = (t*j0l(x)/b);
206 }
207 }
208 if(sgn==1) return -b; else return b;
209 }
210
211 GENERIC ynl(n,x)
212 int n; GENERIC x;{
213 int i;
214 int sign;
215 GENERIC a, b, temp;
216
217 if(x!=x) return x+x;
218 if (x <= zero)
219 if(x==zero)
220 return -one/zero;
221 else
222 return zero/zero;
223 sign = 1;
224 if(n<0){
225 n = -n;
226 if((n&1) == 1) sign = -1;
227 }
228 if(n==0) return(y0l(x));
229 if(n==1) return(sign*y1l(x));
230 if(!finitel(x)) return zero;
231
232 if(x>1.0e91L) { /* x >> n**2
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)
235 Let s=sin(x), c=cos(x),
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
|
55 * that forward recursion is used for all
56 * values of n>1.
57 *
58 */
59
60 #include "libm.h"
61 #include "longdouble.h"
62 #include <float.h> /* LDBL_MAX */
63
64 #define GENERIC long double
65
66 static const GENERIC
67 invsqrtpi= 5.641895835477562869480794515607725858441e-0001L,
68 two = 2.0L,
69 zero = 0.0L,
70 one = 1.0L;
71
72 GENERIC
73 jnl(n,x) int n; GENERIC x;{
74 int i, sgn;
75 GENERIC a, b, temp = 0, z, w;
76
77 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
78 * Thus, J(-n,x) = J(n,-x)
79 */
80 if(n<0){
81 n = -n;
82 x = -x;
83 }
84 if(n==0) return(j0l(x));
85 if(n==1) return(j1l(x));
86 if(x!=x) return x+x;
87 if((n&1)==0)
88 sgn=0; /* even n */
89 else
90 sgn = signbitl(x); /* old n */
91 x = fabsl(x);
92 if(x == zero||!finitel(x)) b = zero;
93 else if((GENERIC)n<=x) { /* Safe to use
94 J(n+1,x)=2n/x *J(n,x)-J(n-1,x)
95 */
195 temp = b;
196 b = ((i+i)/x)*b - a;
197 a = temp;
198 if(b>1e1000L) {
199 a /= b;
200 t /= b;
201 b = 1.0;
202 }
203 }
204 }
205 b = (t*j0l(x)/b);
206 }
207 }
208 if(sgn==1) return -b; else return b;
209 }
210
211 GENERIC ynl(n,x)
212 int n; GENERIC x;{
213 int i;
214 int sign;
215 GENERIC a, b, temp = 0;
216
217 if(x!=x)
218 return x+x;
219 if (x <= zero) {
220 if(x==zero)
221 return -one/zero;
222 else
223 return zero/zero;
224 }
225 sign = 1;
226 if(n<0){
227 n = -n;
228 if((n&1) == 1) sign = -1;
229 }
230 if(n==0) return(y0l(x));
231 if(n==1) return(sign*y1l(x));
232 if(!finitel(x)) return zero;
233
234 if(x>1.0e91L) { /* x >> n**2
235 Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
236 Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
237 Let s=sin(x), c=cos(x),
238 xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
239
240 n sin(xn)*sqt2 cos(xn)*sqt2
241 ----------------------------------
242 0 s-c c+s
243 1 -s-c -c+s
244 2 -s+c -c-s
|