[Haskell-cafe] Fractional sqrt

[ajb at spamcop.net]

G'day all.

Quoting Henning Thielemann <lemming at henning-thielemann.de>:

> 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