- return one/(x+x*x);
- if ((ix|lx) == 0)
- return -one/zero;
- if (hx < 0)
- return zero/zero;
- if (ix >= 0x40000000) { /* |x| >= 2.0 */
- /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
- * where x0 = x-pi/4
- * Better formula:
- * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
- * = 1/sqrt(2) * (sin(x) + cos(x))
- * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
- * = 1/sqrt(2) * (sin(x) - cos(x))
- * To avoid cancellation, use
- * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
- * to compute the worse one.
- */
- s = sin(x);
- c = cos(x);
- ss = s-c;
- cc = s+c;
- /*
- * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
- * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
- */
- if (ix < 0x7fe00000) { /* make sure x+x does not overflow */
- z = -cos(x+x);
- if (s*c < zero)
- cc = z/ss;
- else
- ss = z/cc;
- }
- if (ix > 0x48000000)
- z = (invsqrtpi*ss)/sqrt(x);
- else {
- u = pzero(x);
- v = qzero(x);
- z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
- }
- return z;
+ return 1/x;
+
+ if (ix >= 0x40000000) { /* x >= 2 */
+ /* large ulp errors near zeros: 3.958, 7.086,.. */
+ return common(ix,x,1);