Talk:99 questions/31 to 41

Hi,

there are several problems with the solution to 35, given from "Another possibility is .." on.

First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^32. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n. And there are doubts anyway concerning the price for building such a list first.

Second:

a) If found a 'matching' divisor, we must factor it out completely, that means divide by it, as often as possible. Then the dividend for the following divisions gets possibly smaller and - what's more - the square root of it too.

b) The division by a current divisor is done twice in the solution, once for getting the quotient, a second time for getting the remainder. We must use divMod. This reduces the number of the expensive divisions to the half again.

c) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only. Another reduction by 50%.

d) factor*factor is computed at each step. We should compute a bound ("intsqrt" below) only once for each "successful" division.

All in all these proposals should give improvement in performance in a magnitude of 1000% to the solution in the text, with code still naive and easily readable for beginners in Haskell (as me):

data DivIter = 
    DivIter { dividend :: Integer, 
              divisor  :: Integer,
              bound    :: Integer,
              result   :: [Integer] }

intsqrt m = floor (sqrt $ fromInteger m)

primefactors :: Integer -> [Integer]
primefactors n | n < 2 = []
               | even n = o2 ++ (primefactors o1)
               | otherwise = if z /=1 then revres ++ [z] else revres
                    where 
                        res = divisions (DivIter { dividend = o1,
                                                   divisor = 3,  
                                                   bound = intsqrt(o1),
                                                   result = o2})
                        revres = reverse (result res)
                        z = dividend res  
                               --indeed gets 1 in certain situations
                        (o1, o2) = twosect (n, [])

twosect :: (Integer, [Integer]) -> (Integer, [Integer])
twosect m  | odd (fst m) = m
           | even (fst m) = twosect ( div (fst m) 2,  2:(snd m) )

divstep :: DivIter -> DivIter 
divstep y = if ximod>0 then 
            y { divisor=(divisor y) + 2 } 
        else
            y {dividend = xidiv, 
               bound = intsqrt xidiv, 
               result = (divisor y):(result y) }
    where (xidiv, ximod) = divMod (dividend y) (divisor y)

divisions :: DivIter -> DivIter
divisions y = if divisor y <= bound y then divisions (divstep y) 
              else y

This computes

  primefactors (2^90+1) = [5,5,13,37,41,61,109,181,1321,54001,29247661]

in a moment, as an equivalent python script does, but fails to beat that script for 2^120+1 - due to the recursion in "divisions". Compiled with this version for "divisions" (indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement):

divisions = do
    y <- get
    if divisor y <= bound y then do
        put ( divstep y )
        divisions
        else 
            return y

called from primefactors by

res = execState divisions (DivIter { dividend = o1, 
                                     divisor = 3, 
                                     bound = intsqrt(o1),
                                     result = o2 })

it computes

   primefactors (2^120+1) = [97,257,673,394783681,4278255361,46908728641]

in about 30 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.