201 if(b>1e1000L) {
202 a /= b;
203 t /= b;
204 b = 1.0;
205 }
206 }
207 }
208 b = (t*j0l(x)/b);
209 }
210 }
211 if(sgn==1) return -b; else return b;
212 }
213
214 GENERIC ynl(n,x)
215 int n; GENERIC x;{
216 int i;
217 int sign;
218 GENERIC a, b, temp;
219
220 if(x!=x) return x+x;
221 if (x <= zero)
222 if(x==zero)
223 return -one/zero;
224 else
225 return zero/zero;
226 sign = 1;
227 if(n<0){
228 n = -n;
229 if((n&1) == 1) sign = -1;
230 }
231 if(n==0) return(y0l(x));
232 if(n==1) return(sign*y1l(x));
233 if(!finitel(x)) return zero;
234
235 if(x>1.0e91L) { /* x >> n**2
236 Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
237 Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
238 Let s=sin(x), c=cos(x),
239 xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
240
241 n sin(xn)*sqt2 cos(xn)*sqt2
242 ----------------------------------
243 0 s-c c+s
244 1 -s-c -c+s
245 2 -s+c -c-s
|
201 if(b>1e1000L) {
202 a /= b;
203 t /= b;
204 b = 1.0;
205 }
206 }
207 }
208 b = (t*j0l(x)/b);
209 }
210 }
211 if(sgn==1) return -b; else return b;
212 }
213
214 GENERIC ynl(n,x)
215 int n; GENERIC x;{
216 int i;
217 int sign;
218 GENERIC a, b, temp;
219
220 if(x!=x) return x+x;
221 if (x <= zero) {
222 if(x==zero)
223 return -one/zero;
224 else
225 return zero/zero;
226 }
227 sign = 1;
228 if(n<0){
229 n = -n;
230 if((n&1) == 1) sign = -1;
231 }
232 if(n==0) return(y0l(x));
233 if(n==1) return(sign*y1l(x));
234 if(!finitel(x)) return zero;
235
236 if(x>1.0e91L) { /* x >> n**2
237 Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
238 Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
239 Let s=sin(x), c=cos(x),
240 xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
241
242 n sin(xn)*sqt2 cos(xn)*sqt2
243 ----------------------------------
244 0 s-c c+s
245 1 -s-c -c+s
246 2 -s+c -c-s
|