Fast Square Root Challenge

Design a fast square root function which computes double floating point numbers faster than SQRTSD(SSE2). Bonus for precision.

static inline uintmax_t isqrt(uintmax_t n){
 uintmax_t r = 0;
 // fixed this so it should work even if uintmax_t is an odd number of bits!
 // ... i've never seen an implementation where that's the case, but i realized
 // that the code i had here before would break if someone was perverse enough
 // to do that, so i fixed it.
 for(uintmax_t i = ~(UINTMAX_MAX >> 2) & UINTMAX_MAX / 3; i; i >>= 2)
  if(n >= (i | r)){
   n -= i | r;
   r = r >> 1 | i;
  } else r >>= 1;
 return r;
}

// 6 digits after the decimal point for floating-point-using idiots:
#define SQRT(n) (isqrt((uintmax_t)(n) * 1000000000000ULL /* 100⁶ */) / 1000000.0L /* 10⁶ */)

it's a lot faster than SQRTSD on my machine, because SQRTSD requires at least as much time as it would take to buy and install a new CPU (and probably a new motherboard as well).