54 * values of n>1.
55 *
56 */
57
58 #include "libm.h"
59 #include <float.h> /* DBL_MIN */
60 #include <values.h> /* X_TLOSS */
61 #include "xpg6.h" /* __xpg6 */
62
63 #define GENERIC double
64
65 static const GENERIC
66 invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
67 two = 2.0,
68 zero = 0.0,
69 one = 1.0;
70
71 GENERIC
72 jn(int n, GENERIC x) {
73 int i, sgn;
74 GENERIC a, b, temp;
75 GENERIC z, w, ox, on;
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 ox = x; on = (GENERIC)n;
81 if(n<0){
82 n = -n;
83 x = -x;
84 }
85 if(isnan(x)) return x*x; /* + -> * for Cheetah */
86 if (!((int) _lib_version == libm_ieee ||
87 (__xpg6 & _C99SUSv3_math_errexcept) != 0)) {
88 if(fabs(x) > X_TLOSS) return _SVID_libm_err(on,ox,38);
89 }
90 if(n==0) return(j0(x));
91 if(n==1) return(j1(x));
92 if((n&1)==0)
93 sgn=0; /* even n */
94 else
200 temp = b;
201 b = ((i+i)/x)*b - a;
202 a = temp;
203 if(b>1e100) {
204 a /= b;
205 t /= b;
206 b = 1.0;
207 }
208 }
209 }
210 b = (t*j0(x)/b);
211 }
212 }
213 if(sgn==1) return -b; else return b;
214 }
215
216 GENERIC
217 yn(int n, GENERIC x) {
218 int i;
219 int sign;
220 GENERIC a, b, temp, ox, on;
221
222 ox = x; on = (GENERIC)n;
223 if(isnan(x)) return x*x; /* + -> * for Cheetah */
224 if (x <= zero)
225 if(x==zero)
226 /* return -one/zero; */
227 return _SVID_libm_err((GENERIC)n,x,12);
228 else
229 /* return zero/zero; */
230 return _SVID_libm_err((GENERIC)n,x,13);
231 if (!((int) _lib_version == libm_ieee ||
232 (__xpg6 & _C99SUSv3_math_errexcept) != 0)) {
233 if(x > X_TLOSS) return _SVID_libm_err(on,ox,39);
234 }
235 sign = 1;
236 if(n<0){
237 n = -n;
238 if((n&1) == 1) sign = -1;
239 }
240 if(n==0) return(y0(x));
241 if(n==1) return(sign*y1(x));
242 if(!finite(x)) return zero;
243
244 if(x>1.0e91) { /* x >> n**2
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
249
250 n sin(xn)*sqt2 cos(xn)*sqt2
|
54 * values of n>1.
55 *
56 */
57
58 #include "libm.h"
59 #include <float.h> /* DBL_MIN */
60 #include <values.h> /* X_TLOSS */
61 #include "xpg6.h" /* __xpg6 */
62
63 #define GENERIC double
64
65 static const GENERIC
66 invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
67 two = 2.0,
68 zero = 0.0,
69 one = 1.0;
70
71 GENERIC
72 jn(int n, GENERIC x) {
73 int i, sgn;
74 GENERIC a, b, temp = 0;
75 GENERIC z, w, ox, on;
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 ox = x; on = (GENERIC)n;
81 if(n<0){
82 n = -n;
83 x = -x;
84 }
85 if(isnan(x)) return x*x; /* + -> * for Cheetah */
86 if (!((int) _lib_version == libm_ieee ||
87 (__xpg6 & _C99SUSv3_math_errexcept) != 0)) {
88 if(fabs(x) > X_TLOSS) return _SVID_libm_err(on,ox,38);
89 }
90 if(n==0) return(j0(x));
91 if(n==1) return(j1(x));
92 if((n&1)==0)
93 sgn=0; /* even n */
94 else
200 temp = b;
201 b = ((i+i)/x)*b - a;
202 a = temp;
203 if(b>1e100) {
204 a /= b;
205 t /= b;
206 b = 1.0;
207 }
208 }
209 }
210 b = (t*j0(x)/b);
211 }
212 }
213 if(sgn==1) return -b; else return b;
214 }
215
216 GENERIC
217 yn(int n, GENERIC x) {
218 int i;
219 int sign;
220 GENERIC a, b, temp = 0, ox, on;
221
222 ox = x; on = (GENERIC)n;
223 if(isnan(x)) return x*x; /* + -> * for Cheetah */
224 if (x <= zero) {
225 if(x==zero) {
226 /* return -one/zero; */
227 return _SVID_libm_err((GENERIC)n,x,12);
228 } else {
229 /* return zero/zero; */
230 return _SVID_libm_err((GENERIC)n,x,13);
231 }
232 }
233 if (!((int) _lib_version == libm_ieee ||
234 (__xpg6 & _C99SUSv3_math_errexcept) != 0)) {
235 if(x > X_TLOSS) return _SVID_libm_err(on,ox,39);
236 }
237 sign = 1;
238 if(n<0){
239 n = -n;
240 if((n&1) == 1) sign = -1;
241 }
242 if(n==0) return(y0(x));
243 if(n==1) return(sign*y1(x));
244 if(!finite(x)) return zero;
245
246 if(x>1.0e91) { /* x >> n**2
247 Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
248 Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
249 Let s=sin(x), c=cos(x),
250 xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
251
252 n sin(xn)*sqt2 cos(xn)*sqt2
|