Talk:99 questions/31 to 41

Hi,

there are several problems with the solution to 35.

First: The signature of primefactors is Integer -> [Integer]. That means, the function is designed to work on integers larger than 2^64. (a reasonable task now, when ghc outperforms Fortran). But then we cannot make a list of the primes up to sqrt n, which exceeds 2^32. 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 it is done completely. Thus we must go through the divisors only once for the whole process (not start by dividing by 2 again after every "successful" division). This reduces the number of computations at least to the half.

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)n*n is computed at each step. We should compute a bound ("intsqrt") 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 {x      :: Integer, 
             divi   :: 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 {x    = o1, 
                                                      divi = 3, 
                                                      bound = intsqrt(o1), 
                                                      result = o2})
                            revres = reverse (result res)
                            z = x res  -- yes, sometimes equal to 1
                            (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 { divi=(divi y) +2 } else
    y {x = xidiv, bound = intsqrt xidiv, result = (divi y):(result y) }
    where
        (xidiv, ximod) = divMod (x y) (divi y)

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

The reader should notice, that "twosect" is an easy model for the (still easy) function "divstep".

This proposal assumes, that named fields really are syntactic sugar and will be completely factored out by compilation. Otherwise a quadruple should be used.

"divisions" is separated from "divstep" here not for readability, but for possible replacement. We MUST ensure tail recursion here, or use a monadic version of "divisions" with destructive updates - but that's a more advanced topic.