12 /* j1(x), y1(x)
22  * 		y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
38  * Method -- y1(x):
39  *	1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
42  *		y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
43  *	   therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
44  *	   We use the following function to approximate y1,
45  *		y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
49  *	   Note: For tiny x, 1/x dominate y1 and hence
50  *		y1(tiny) = -2/pi/tiny, (choose tiny<2**-54)
52  * 		y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
103 	 * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)  in j1()
139 y1(double x)  in y1()  function
147 	 * y1(NaN) = NaN.  in y1()
148 	 * y1(Inf) = 0.  in y1()
149 	 * y1(-Inf) = NaN and raise invalid exception.  in y1()
152 	/* y1(+-0) = -inf and raise divide-by-zero exception. */  in y1()
154 	/* y1(x<0) = NaN and raise invalid exception. */  in y1()
165         /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))  in y1()