41 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
42 * Note 2. About jn(n,x), yn(n,x)
43 * For n=0, j0(x) is called,
44 * for n=1, j1(x) is called,
45 * for n<x, forward recursion us used starting
46 * from values of j0(x) and j1(x).
47 * for n>x, a continued fraction approximation to
48 * j(n,x)/j(n-1,x) is evaluated and then backward
49 * recursion is used starting from a supposed value
50 * for j(n,x). The resulting value of j(0,x) is
51 * compared with the actual value to correct the
52 * supposed value of j(n,x).
53 *
54 * yn(n,x) is similar in all respects, except
55 * that forward recursion is used for all
56 * values of n>1.
57 *
58 */
59
60 #include "libm.h"
61 #include <float.h> /* LDBL_MAX */
62
63 #define GENERIC long double
64
65 static const GENERIC
66 invsqrtpi= 5.641895835477562869480794515607725858441e-0001L,
67 two = 2.0L,
68 zero = 0.0L,
69 one = 1.0L;
70
71 GENERIC
72 jnl(n,x) int n; GENERIC x;{
73 int i, sgn;
74 GENERIC a, b, temp, z, w;
75
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)
78 */
79 if(n<0){
80 n = -n;
81 x = -x;
82 }
83 if(n==0) return(j0l(x));
84 if(n==1) return(j1l(x));
85 if(x!=x) return x+x;
86 if((n&1)==0)
87 sgn=0; /* even n */
88 else
89 sgn = signbitl(x); /* old n */
90 x = fabsl(x);
91 if(x == zero||!finitel(x)) b = zero;
92 else if((GENERIC)n<=x) { /* Safe to use
93 J(n+1,x)=2n/x *J(n,x)-J(n-1,x)
94 */
194 temp = b;
195 b = ((i+i)/x)*b - a;
196 a = temp;
197 if(b>1e1000L) {
198 a /= b;
199 t /= b;
200 b = 1.0;
201 }
202 }
203 }
204 b = (t*j0l(x)/b);
205 }
206 }
207 if(sgn==1) return -b; else return b;
208 }
209
210 GENERIC ynl(n,x)
211 int n; GENERIC x;{
212 int i;
213 int sign;
214 GENERIC a, b, temp;
215
216 if(x!=x) return x+x;
217 if (x <= zero)
218 if(x==zero)
219 return -one/zero;
220 else
221 return zero/zero;
222 sign = 1;
223 if(n<0){
224 n = -n;
225 if((n&1) == 1) sign = -1;
226 }
227 if(n==0) return(y0l(x));
228 if(n==1) return(sign*y1l(x));
229 if(!finitel(x)) return zero;
230
231 if(x>1.0e91L) { /* 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)
234 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
|
41 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
42 * Note 2. About jn(n,x), yn(n,x)
43 * For n=0, j0(x) is called,
44 * for n=1, j1(x) is called,
45 * for n<x, forward recursion us used starting
46 * from values of j0(x) and j1(x).
47 * for n>x, a continued fraction approximation to
48 * j(n,x)/j(n-1,x) is evaluated and then backward
49 * recursion is used starting from a supposed value
50 * for j(n,x). The resulting value of j(0,x) is
51 * compared with the actual value to correct the
52 * supposed value of j(n,x).
53 *
54 * yn(n,x) is similar in all respects, except
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
|