[Haskell-cafe] Fractional sqrt

[Yitzchak Gale]

Sorry folks, it is just wrong to use Newton's method
for Rational.

Andrew Bromage wrote:
> First off, note that for fractions, sqrt(p/q) = sqrt p / sqrt q.

Don't do that for Rational - you lose precious precision.

The whole idea for calculations in Rational is to find
a lucky denominator that is not too big but gives you
a good approximation. Newton's method will converge
very quickly - and it will also blow up your stack
very quickly as the denominator increases exponentially.

You always get the best approximation for a Rational
using continued fractions. For most calculations that
is easier said than done, but for square root we are
very lucky.

Is the continued fraction algorithm fast enough?
Well, on my old machine, the following code computes
sqrt 2 to within 10^(-1000) with no noticeable delay.
This can obviously be optimized in many ways,
but I am not sure it is worth the effort.

Regards,
Yitz

\begin{code}

-- Arbitrary precision square root for Rational

module SquareRoot where

import Data.Ratio

-- The square root of x, within error tolerance e
-- Warning: diverges if e=0 and x is not a perfect square
rSqrt :: Rational -> Rational -> Rational
rSqrt x _ | x < 0 = error "rSqrt of negative"
rSqrt _ e | e < 0 = error "rSqrt negative tolerance"
rSqrt x e = firstWithinE $ map (uncurry (%)) $ convergents $
            sqrtCF (numerator x) (denominator x)
  where
    firstWithinE = head . dropWhile off
    off y = abs (x - y^2 - e^2) > 2 * e * y

-- The continued fraction for the square root of positive p%q
sqrtCF :: Integer -> Integer -> [Integer]
sqrtCF p q = cf 1 0 1
  where
    -- Continued fraction for (a*sqrt(p%q)+b)/c
    cf a b c
     | c == 0    = []
     | otherwise = x : cf (a' `div` g) (b' `div` g) (c' `div` g)
     where
       x  = (iSqrt (p*a*a) q + b) `div` c
       a' = a * c * q
       e  = c * x - b
       b' = c * e * q
       c' = p*a*a - q*e*e
       g  = foldl1 gcd [a', b', c']

-- The greatest integer less than or equal to sqrt (p%q) for positive p%q
iSqrt :: Integer -> Integer -> Integer
iSqrt p q = sq 0 (1 + p `div` q)
  where
    sq y z
     | z - y < 2      = y
     | q * m * m <= p = sq m z
     | otherwise      = sq y m
     where m = (y + z) `div` 2

-- The following could go in a separate module
-- for general continued fractions

-- The convergents of a continued fraction
convergents :: [Integer] -> [(Integer, Integer)]
convergents cf = drop 2 $ zip h k
  where
    h = 0 : 1 : zipWith3 affine cf (drop 1 h) h
    k = 1 : 0 : zipWith3 affine cf (drop 1 k) k
    affine x y z = x * y + z

/end{code}