1:
2: 3: 4: 5: 6: 7: 8: 9: 10: 11:
12:
13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: 32: 33: 34: 35: 36: 37:
38:
39: #include "math.h"
40: #include "math_private.h"
41:
42: static const double
43: invsqrtpi= 5.64189583547756279280e-01,
44: two = 2.00000000000000000000e+00,
45: one = 1.00000000000000000000e+00;
46:
47: static const double zero = 0.00000000000000000000e+00;
48:
49: double
50: jn(int n, double x)
51: {
52: int32_t i,hx,ix,lx, sgn;
53: double a, b, temp, di;
54: double z, w;
55:
56: 57: 58:
59: EXTRACT_WORDS(hx,lx,x);
60: ix = 0x7fffffff&hx;
61:
62: if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
63: if(n<0){
64: n = -n;
65: x = -x;
66: hx ^= 0x80000000;
67: }
68: if(n==0) return(j0(x));
69: if(n==1) return(j1(x));
70: sgn = (n&1)&(hx>>31);
71: x = fabs(x);
72: if((ix|lx)==0||ix>=0x7ff00000)
73: b = zero;
74: else if((double)n<=x) {
75:
76: if(ix>=0x52D00000) {
77: 78: 79: 80: 81: 82: 83: 84: 85: 86: 87: 88: 89:
90: switch(n&3) {
91: case 0: temp = cos(x)+sin(x); break;
92: case 1: temp = -cos(x)+sin(x); break;
93: case 2: temp = -cos(x)-sin(x); break;
94: case 3: temp = cos(x)-sin(x); break;
95: }
96: b = invsqrtpi*temp/sqrt(x);
97: } else {
98: a = j0(x);
99: b = j1(x);
100: for(i=1;i<n;i++){
101: temp = b;
102: b = b*((double)(i+i)/x) - a;
103: a = temp;
104: }
105: }
106: } else {
107: if(ix<0x3e100000) {
108: 109: 110:
111: if(n>33)
112: b = zero;
113: else {
114: temp = x*0.5; b = temp;
115: for (a=one,i=2;i<=n;i++) {
116: a *= (double)i;
117: b *= temp;
118: }
119: b = b/a;
120: }
121: } else {
122:
123: 124: 125: 126: 127: 128: 129: 130: 131: 132: 133: 134: 135: 136: 137: 138: 139: 140: 141: 142: 143: 144: 145: 146: 147: 148: 149:
150:
151: double t,v;
152: double q0,q1,h,tmp; int32_t k,m;
153: w = (n+n)/(double)x; h = 2.0/(double)x;
154: q0 = w; z = w+h; q1 = w*z - 1.0; k=1;
155: while(q1<1.0e9) {
156: k += 1; z += h;
157: tmp = z*q1 - q0;
158: q0 = q1;
159: q1 = tmp;
160: }
161: m = n+n;
162: for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
163: a = t;
164: b = one;
165: 166: 167: 168: 169: 170: 171: 172:
173: tmp = n;
174: v = two/x;
175: tmp = tmp*log(fabs(v*tmp));
176: if(tmp<7.09782712893383973096e+02) {
177: for(i=n-1,di=(double)(i+i);i>0;i--){
178: temp = b;
179: b *= di;
180: b = b/x - a;
181: a = temp;
182: di -= two;
183: }
184: } else {
185: for(i=n-1,di=(double)(i+i);i>0;i--){
186: temp = b;
187: b *= di;
188: b = b/x - a;
189: a = temp;
190: di -= two;
191:
192: if(b>1e100) {
193: a /= b;
194: t /= b;
195: b = one;
196: }
197: }
198: }
199: b = (t*j0(x)/b);
200: }
201: }
202: if(sgn==1) return -b; else return b;
203: }
204:
205: double
206: yn(int n, double x)
207: {
208: int32_t i,hx,ix,lx;
209: int32_t sign;
210: double a, b, temp;
211:
212: EXTRACT_WORDS(hx,lx,x);
213: ix = 0x7fffffff&hx;
214:
215: if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
216: if((ix|lx)==0) return -one/zero;
217: if(hx<0) return zero/zero;
218: sign = 1;
219: if(n<0){
220: n = -n;
221: sign = 1 - ((n&1)<<1);
222: }
223: if(n==0) return(y0(x));
224: if(n==1) return(sign*y1(x));
225: if(ix==0x7ff00000) return zero;
226: if(ix>=0x52D00000) {
227: 228: 229: 230: 231: 232: 233: 234: 235: 236: 237: 238: 239:
240: switch(n&3) {
241: case 0: temp = sin(x)-cos(x); break;
242: case 1: temp = -sin(x)-cos(x); break;
243: case 2: temp = -sin(x)+cos(x); break;
244: case 3: temp = sin(x)+cos(x); break;
245: }
246: b = invsqrtpi*temp/sqrt(x);
247: } else {
248: u_int32_t high;
249: a = y0(x);
250: b = y1(x);
251:
252: GET_HIGH_WORD(high,b);
253: for(i=1;i<n&&high!=0xfff00000;i++){
254: temp = b;
255: b = ((double)(i+i)/x)*b - a;
256: GET_HIGH_WORD(high,b);
257: a = temp;
258: }
259: }
260: if(sign>0) return b; else return -b;
261: }