# Talk:99 questions/31 to 41

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

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.