[Haskell-cafe] Fractional sqrt
ajb at spamcop.net
ajb at spamcop.net
Fri Jan 19 10:48:26 EST 2007
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
More information about the Haskell-Cafe
mailing list