/** Estimate 1/sqrt[x] */
estimateInvSqrt[x, digits=getPrecision[]] :=
{
   origPrec = getPrecision[]
   try
   {
      setPrecision[16]
      ax = abs[x]
      if ax < 1e307 && ax > 1e-307
         xn = 1/sqrt[x]   // Can use built-in sqrt operator
      else
         xn = 2^(-approxLog2[x]/2) // Maybe < 1/3 of a decimal digit accuracy!
      
      setPrecision[digits + 1]
      // println[xn]
      lastx = xn
      do
      {
         lastx = xn
         xn = xn/2 (3 - x xn^2)
         // println["$xn"]
      } until lastx == xn

      return 1. xn
   }
   finally
      setPrecision[origPrec]
}