Difference between revisions of "Talk:99 questions/31 to 41"
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)) |
|||
(15 intermediate revisions by one other user not shown) | |||
Line 1: | Line 1: | ||
Hi, |
Hi, |
||
− | there are several problems with the solution to 35. |
+ | 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^ |
+ | 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: |
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. |
||
⚫ | |||
+ | <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): |
||
⚫ | |||
− | |||
⚫ | |||
− | |||
⚫ | |||
⚫ | |||
⚫ | |||
⚫ | |||
⚫ | |||
⚫ | |||
intsqrt m = floor (sqrt $ fromInteger m) |
intsqrt m = floor (sqrt $ fromInteger m) |
||
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 |
|
+ | (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, |
+ | | 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 } |
− | + | |otherwise = y {dividend = xidiv, |
|
+ | bound = intsqrt xidiv, |
||
⚫ | |||
− | + | result = result 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 |
||
+ | 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: |
||
+ | <haskell>divisions = do |
||
+ | y <- get |
||
+ | if divisor y <= bound y then do |
||
+ | put ( divstep y ) |
||
+ | divisions |
||
+ | else |
||
+ | return y</haskell> |
||
+ | called from primefactors by |
||
− | The reader should notice, that "twosect" is an easy model for the (still easy) function "divstep". |
||
+ | <haskell>res = execState divisions (DivIter { dividend = o1, |
||
− | This proposal assumes, that named fields really are syntactic sugar and will be completely factored out by compilation. Otherwise a quadruple should be used. |
||
⚫ | |||
+ | bound = intsqrt(o1), |
||
+ | result = o2 })</haskell> |
||
+ | it computes |
||
− | "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. |
||
+ | 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. |
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.