Christian Sievers writes (in the Haskell Cafe):
sqrtInt n = help n where help x = let y = ((x + (n `div` x)) `div` 2) in if y
{- following p. 38f of Henri Cohen, A Course in Computational Algebraic Number Theory, where a proof and a suggestion for improvement (choose a better start value) is given.
For large values of n there's a faster Karatsuba variant by Paul Zimmermann, see http://www.inria.fr/rrrt/rr-3805.html.
As far as I know, ghc uses gmp, so I wonder if there is access to functions like mpz_sqrt or mpz_perfect_square_p. -}
The sqrt functions in GMP use Paul Zimmerman's algorithm. The original poster asked for a proof, and interestingly enough the C-implementation in GMP was checked with the proof assistant Coq, see http://www.inria.fr/rrrt/rr-4475.html. This is quite a feat, because GMP's implementation doesn't allocate any intermediate integers. The formal proof has to verify that all array accesses are legal and that the intermediate values don't overlap. The proof is 10000 lines of Coq. The report also contains a functional version of the algorithm with a proof in Coq. This proof is (of course!) much shorter with 500 lines. We can directly convert the functional Coq source from the article (page 14) to Haskell by simply ignoring all proof certificates:
bitsPerLimb = 32 -- bitSize (0 :: Int) -- should be even beta = 2^bitsPerLimb -- beta is the base decompose n = let l = n `div` 2 in (n - l, l) fix f = f (fix f) sqrtCoq = fix (flip sqrt_F) sqrt_F h sqrt n | h <= 1 = normalised_base_case h n | otherwise = let (h', l) = decompose h (n', n'') = divMod n (beta^(2*l)) (n1, n0) = divMod n'' (beta^l) (s', r') = sqrt h' n' (q, r'') = divMod (r' * beta^l + n1) (2 * s') in if (0 <= r'' * beta^l + n0 - q^2) then (s' * beta^l + q, r'' * beta^l + n0 - q^2) else (s' * beta^l + q - 1, r'' * beta^l + n0 - q^2 + 2 * (s' * beta^l + q) - 1)
Zimmermann's algorithm needs an auxiliary procedure to compute sqrts of small integers, let's use Christian's version:
normalised_base_case _ n = let s = sqrtInt n in (s, n-s*s)
and it must do some preconditioning (see the reports for details):
sqrtRem :: Integer -> (Integer, Integer) sqrtRem n | n < 0 = error "sqrtRem neg" | n == 0 = (0, 0) | otherwise = let (h, msl) = limbs 0 n mul0 = normalise 1 msl mul = (if odd h then beta else 1) * mul0 n' = n*mul*mul (s0, r0) = sqrtCoq (fst (decompose h)) n' r1 = s0 `mod` mul in (s0 `div` mul, (r0+2*s0*r1) `div` (mul*mul)) where -- count limbs and return most significant limb limbs :: Integer -> Integer -> (Integer, Integer) limbs size n | n < beta = (size, n) | otherwise = limbs (size+1) (n `div` beta) -- shift left until either of the top two bits is set normalise :: Integer -> Integer -> Integer normalise m msl | msl < beta `div` 4 = normalise (m*2) (msl*4) | otherwise = m
main = print $ sqrtRem (2^100000)