- return 1.0f/(x+x*x);
- if (ix == 0)
- return -1.0f/0.0f;
- if (hx < 0)
- return 0.0f/0.0f;
- if (ix >= 0x40000000) { /* |x| >= 2.0 */
- s = sinf(x);
- c = cosf(x);
- ss = -s-c;
- cc = s-c;
- if (ix < 0x7f000000) { /* make sure x+x not overflow */
- z = cosf(x+x);
- if (s*c > 0.0f)
- cc = z/ss;
- else
- ss = z/cc;
- }
- /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
- * where x0 = x-3pi/4
- * Better formula:
- * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
- * = 1/sqrt(2) * (sin(x) - cos(x))
- * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
- * = -1/sqrt(2) * (cos(x) + sin(x))
- * To avoid cancellation, use
- * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
- * to compute the worse one.
- */
- if (ix > 0x48000000)
- z = (invsqrtpi*ss)/sqrtf(x);
- else {
- u = ponef(x);
- v = qonef(x);
- z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
- }
- return z;
- }
- if (ix <= 0x24800000) /* x < 2**-54 */
+ return 1/x;
+ if (ix >= 0x40000000) /* |x| >= 2.0 */
+ return common(ix,x,1,0);
+ if (ix < 0x32000000) /* x < 2**-27 */