/* Compute 1.0f / a almost correctly rounded. Halley iteration with cubic convergence */
float rcpss(float in)
{
    float r, out;
    __asm__ __volatile__("rcpss %1, %0\n" :"=x"(r) :"x"(in));
    //out = r*out*(1.0 + out);
    out = (r - (r*r*in));
    out = out*(r+out); /* Halley step. */
    return r+out;
}