[Haskell-cafe] puzzle: prove this floorSqrt correct

[Ronny Wichers Schreur]

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<x then help y else x
 >
 > {-
 > 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)