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) We first deal with the (possibly multiple) factor 2, then go through the odd numbers as divisors only.

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

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 result res ++ [z] else result res
       where 
         res = divisions (DivIter { dividend = o1,
                                    divisor = 3,  
                                    bound = intsqrt(o1),
                                    result = o2})
         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, snd m ++ [2] )

divstep :: DivIter -> DivIter 
divstep y |ximod > 0 = y { divisor=(divisor y) + 2 } 
          |otherwise = y {dividend = xidiv, 
                          bound = intsqrt xidiv, 
                          result = result y ++ [divisor y]}
    where (xidiv, ximod) = divMod (dividend y) (divisor y)

divisions :: DivIter -> DivIter
divisions y |divisor y <= bound y = divisions (divstep y) 
            |otherwise            = 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". Indeed the function "divisions" is separated from "divstep" here not for readability, but for possible replacement. Compiled with this version:

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 20 minutes (with a constant footprint of 2MB) on my system, which still is not very good compared with the two hours python needs.