- else {
- /*
- * if allow error up to 2 ulp, simply return
- * -1.0 / (x+r) here
- */
- /* compute -1.0 / (x+r) accurately */
- double a, t;
- z = w;
- SET_LOW_WORD(z,0);
- v = r - (z - x); /* z+v = r+x */
- t = a = -1.0 / w; /* a = -1.0/w */
- SET_LOW_WORD(t,0);
- s = 1.0 + t * z;
- return t + a * (s + t * v);
- }
+ /* -1.0/(x+r) has up to 2ulp error, so compute it accurately */
+ w0 = w;
+ SET_LOW_WORD(w0, 0);
+ v = r - (w0 - x); /* w0+v = r+x */
+ a0 = a = -1.0 / w;
+ SET_LOW_WORD(a0, 0);
+ return a0 + a*(1.0 + a0*w0 + a0*v);