Difference between revisions of "Talk:99 questions/31 to 41"
m |
m (shorten code lines to make it readable on lower resolutions (<haskell> texts get scroll bars ON TOP of code, absolutely UNREADABLE, please someone FIX IT)) |
||
(8 intermediate revisions by one other user not shown) | |||
Line 7: | Line 7: | ||
Second: |
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. |
||
⚫ | |||
+ | <haskell> |
||
⚫ | |||
⚫ | |||
− | |||
− | 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): |
||
− | |||
− | |||
⚫ | |||
DivIter { dividend :: Integer, |
DivIter { dividend :: Integer, |
||
divisor :: Integer, |
divisor :: Integer, |
||
Line 27: | Line 22: | ||
primefactors :: Integer -> [Integer] |
primefactors :: Integer -> [Integer] |
||
− | primefactors n |
+ | primefactors n |
− | + | | n < 2 = [] |
|
− | + | | even n = o2 ++ (primefactors o1) |
|
− | + | | otherwise = if z /=1 then result res ++ [z] else result res |
|
⚫ | |||
⚫ | |||
− | + | res = divisions (DivIter { dividend = o1, |
|
− | + | divisor = 3, |
|
− | + | bound = intsqrt(o1), |
|
− | + | result = o2}) |
|
+ | z = dividend res |
||
--indeed gets 1 in certain situations |
--indeed gets 1 in certain situations |
||
− | + | (o1, o2) = twosect (n, []) |
|
twosect :: (Integer, [Integer]) -> (Integer, [Integer]) |
twosect :: (Integer, [Integer]) -> (Integer, [Integer]) |
||
− | twosect m | odd (fst m) = m |
+ | twosect m | odd (fst m) = m |
| even (fst m) = twosect ( div (fst m) 2, snd m ++ [2] ) |
| even (fst m) = twosect ( div (fst m) 2, snd m ++ [2] ) |
||
divstep :: DivIter -> DivIter |
divstep :: DivIter -> DivIter |
||
− | divstep y |
+ | divstep y |ximod > 0 = y { divisor=(divisor y) + 2 } |
− | y { |
+ | |otherwise = y {dividend = xidiv, |
⚫ | |||
⚫ | |||
− | + | result = result y ++ [divisor y]} |
|
− | bound = intsqrt xidiv, |
||
− | result = result y ++ [divisor y] } |
||
where (xidiv, ximod) = divMod (dividend y) (divisor y) |
where (xidiv, ximod) = divMod (dividend y) (divisor y) |
||
divisions :: DivIter -> DivIter |
divisions :: DivIter -> DivIter |
||
− | divisions y |
+ | divisions y |divisor y <= bound y = divisions (divstep y) |
− | + | |otherwise = y</haskell> |
|
+ | |||
This computes |
This computes |
||
− | primefactors (2^90+1) = |
+ | 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". |
+ | 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: |
<haskell>divisions = do |
<haskell>divisions = do |
||
Line 75: | Line 71: | ||
it computes |
it computes |
||
− | primefactors (2^120+1) = |
+ | primefactors (2^120+1) = |
+ | [97,257,673,394783681,4278255361,46908728641] |
||
− | in about |
+ | 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. |
Latest revision as of 14:28, 15 July 2011
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.