G'day all.
Quoting Henning Thielemann
Newton method for sqrt is very fast. It converges quadratically, that is in each iteration the number of correct digits doubles. The problem is to find a good starting approximation.
Yup. So how might we go about doing this? First off, note that for fractions, sqrt(p/q) = sqrt p / sqrt q. Secondly, note that for floating point numbers, sqrt(m * b^e) = sqrt m * b^(e/2). So an approach to get an initial approximation to p/q might be: 1. Compute approximate square roots for p and q, then divide. 2. For each of p and q, express as a floating-point number m * b^(2e), compute an approximate square root for m (which will be less than 2b), then multiply by b^e. We've now reduced the range of approximate square roots that we want to numbers in the range 0 to 2b, where b is a parameter that you can use to tune the algorithm. OK, so let's pick b=16 out of the air. Our code might look something like this: sqrtApprox :: Rational -> Rational sqrtApprox x = sqrtApprox' (numerator x) / sqrtApprox' (denominator x) sqrtApprox' :: Integer -> Rational sqrtApprox' n | n < 0 = error "sqrtApprox'" | otherwise = approx n 1 where approx n acc | n < 256 = (acc%1) * approxSmallSqrt (fromIntegral n) | otherwise = approx (n `div` 256) (acc*16) approxSmallSqrt :: Int -> Rational approxSmallSqrt n = {- detail omitted -} The approxSmallSqrt function is a good candidate for the memoising CAFs approach: http://haskell.org/hawiki/MemoisingCafs You can populate the memo table by using a fixed number of iterations of the Newton-Raphson algorithm on a suitably sane initial estimate. You can get a slightly better estimate by being a bit more careful about rounding. Suppose that n = 2 + 256 * n', then a good estimate for sqrt n is 16 * sqrt n'. But if n = 255 + 256 * n', the algorithm above would also estimate sqrt n as 16 * sqrt n'. It's only slightly more expensive (and occasionally better) to round up instead, and use 16 * sqrt (n'+1) as the estimate. So you might want to do this instead: sqrtApprox' :: Integer -> Rational sqrtApprox' n | n < 0 = error "sqrtApprox'" | otherwise = approx n 1 where approx n acc | n < 272 = (acc%1) * approxSmallSqrt (fromIntegral n) | r < 16 = approx q (acc*16) | otherwise = approx (q+1) (acc*16) where (q,r) = n `divMod` 256 approxSmallSqrt :: Int -> Rational approxSmallSqrt n = {- detail omitted -} I've also extended the range for approxSmallSqrt here from (0,255) to (0,271). It is left as an exercise as to why this might be a good idea. (Hint: 272 is approximately 16.5*16.5.) Cheers, Andrew Bromage