- /* 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 = sinf(x);
- c = cosf(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 < 0x7f000000) { /* make sure x+x not overflow */
- z = -cosf(x+x);
- if (s*c < zero)
- cc = z/ss;
- else
- ss = z/cc;
- }
- if (ix > 0x80000000)
- z = (invsqrtpi*ss)/sqrtf(x);
- else {
- u = pzerof(x);
- v = qzerof(x);
- z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
- }
- return z;