On Thu, 12 Aug 2004, Remi Turk wrote:

On Thu, Aug 12, 2004 at 06:59:26PM +0200, Christian Sievers wrote:

...

blaetterrascheln@web.de wrote:

...

-- Here's the discrete version of Newton's method for finding -- the square root. Does it always work? Any literature?

I recently used, without range check,

sqrtInt n = help n where help x = let y = ((x + (n `div` x)) `div` 2) in if y

If I urgently need factors of an integer I check "factor*factor > n" rather than "factor > intSqrt n". :-]

On Thu, 12 Aug 2004, Remi Turk wrote:

On Thu, Aug 12, 2004 at 09:01:03PM +0200, Henning Thielemann wrote:

...

On Thu, 12 Aug 2004, Remi Turk wrote:

...

I usually (each time I urgently need to calculate primes ;)) use a simple intSqrt = floor . sqrt . fromIntegral (which will indeed give wrong answers for big numbers)

If I urgently need factors of an integer I check "factor*factor > n" rather than "factor > intSqrt n". :-]

but you'll have to (^2) once every "iteration". The following code only has to sqrt once.

Right. Hm, but if you get a 'div' for free for every 'mod' (e.g. divMod) you could also check "factor > n `div` factor"

Oh, the joy of premature optimization. Or even premature coding ;)

-- Will lie when given stupid questions isPrime 1 = False isPrime 2 = True isPrime n = not $ any (n`isDivBy`) (2:[3,5..intSqrt n]) where n `isDivBy` d = n `rem` d == 0 intSqrt = floor . sqrt . fromIntegral

If the numbers become too big, will 'fromIntegral' round down or up? If it rounds down, then intSqrt might lead to a too small result.

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)