18  *	y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
19  *	y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
25  *	For n>x, a continued fraction approximation to
26  *	j(n,x)/j(n-1,x) is evaluated and then backward
43 invsqrtpi=  5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
54 	double z, w;  in jn()  local
56     /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)  in jn()
57      * Thus, J(-n,x) = J(n,-x)  in jn()
62 	if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;  in jn()
64 		n = -n;  in jn()
65 		x = -x;  in jn()
70 	sgn = (n&1)&(hx>>31);	/* even n -- 0, odd n -- sign(x) */  in jn()
75 		/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */  in jn()
78      *	    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)  in jn()
79      *	    Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)  in jn()
81      *		xn=x-(2n+1)*pi/4, sqt2 = sqrt(2), then  in jn()
84      *		----------------------------------  in jn()
85      *		   0	 s-c		 c+s  in jn()
86      *		   1	-s-c 		-c+s  in jn()
87      *		   2	-s+c		-c-s  in jn()
88      *		   3	 s+c		 c-s  in jn()
93 		    case 1: temp = -c+s; break;  in jn()
94 		    case 2: temp = -c-s; break;  in jn()
95 		    case 3: temp =  c-s; break;  in jn()
103 		    b = b*((double)(i+i)/x) - a; /* avoid underflow */  in jn()
108 	    if(ix<0x3e100000) {	/* x < 2**-29 */  in jn()
110      * J(n,x) = 1/n!*(x/2)^n  - ...  in jn()
125 		 *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....  in jn()
126 		 *			2n  - 2(n+1) - 2(n+2)  in jn()
129 		 *  (for large x)   =  ----  ------   ------   .....  in jn()
131 		 *			-- - ------ - ------ -  in jn()
135 		 * is equal to the continued fraction:  in jn()
137 		 *	= -----------------------  in jn()
139 		 *	   w - -----------------  in jn()
141 		 * 	        w+h - ---------  in jn()
142 		 *		       w+2h - ...  in jn()
145 		 * Q(0) = w, Q(1) = w(w+h) - 1,  in jn()
146 		 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),  in jn()
155 		q0 = w;  z = w+h; q1 = w*z - 1.0; k=1;  in jn()
157 			k += 1; z += h;  in jn()
158 			tmp = z*q1 - q0;  in jn()
163 		for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);  in jn()
178 	    	    for(i=n-1,di=(double)(i+i);i>0;i--){  in jn()
181 			b  = b/x - a;  in jn()
183 			di -= two;  in jn()
186 	    	    for(i=n-1,di=(double)(i+i);i>0;i--){  in jn()
189 			b  = b/x - a;  in jn()
191 			di -= two;  in jn()
200 		z = j0(x);  in jn()
202 		if (fabs(z) >= fabs(w))  in jn()
203 		    b = (t*z/b);  in jn()
208 	if(sgn==1) return -b; else return b;  in jn()
221 	if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;  in yn()
222 	/* yn(n,+-0) = -inf and raise divide-by-zero exception. */  in yn()
223 	if((ix|lx)==0) return -one/vzero;  in yn()
228 		n = -n;  in yn()
229 		sign = 1 - ((n&1)<<1);  in yn()
236      *	    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)  in yn()
237      *	    Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)  in yn()
239      *		xn=x-(2n+1)*pi/4, sqt2 = sqrt(2), then  in yn()
242      *		----------------------------------  in yn()
243      *		   0	 s-c		 c+s  in yn()
244      *		   1	-s-c 		-c+s  in yn()
245      *		   2	-s+c		-c-s  in yn()
246      *		   3	 s+c		 c-s  in yn()
250 		    case 0: temp =  s-c; break;  in yn()
251 		    case 1: temp = -s-c; break;  in yn()
252 		    case 2: temp = -s+c; break;  in yn()
260 	/* quit if b is -inf */  in yn()
264 		b = ((double)(i+i)/x)*b - a;  in yn()
269 	if(sign>0) return b; else return -b;  in yn()