G'day all. Quoting Henning Thielemann <[EMAIL PROTECTED]>:
> 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 _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe