Hi, since about three weeks i am learning Haskell now. One of my first excercises is to decompose an Integer into its primefactors. I already posted discussion on the solution to the problem 35 in "99 excercises". My simple algorithm uses a datatype DivIter of 4 named fields together with the core iteration divstep :: DivIter -> DivIter divstep x | divisor x > bound x = x | ximod > 0 = x { divisor = (divisor x) +2 } | otherwise = x {dividend=xidiv, bound=intsqrt(xidiv), result = result x ++ [divisor x] } where (xidiv, ximod) = divMod (dividend x) (divisor x) (dividend x is already odd, when this is called). The problem to solve for really large Integers is how to call divstep iterated without not accumulating billions of stack frames. Here is one possibility: divisions = do y <- get if divisor y <= bound y then do put ( divstep y ) divisions else return y (this for a version of divstep without the first guard) called from res = execState divisions (DivIter { dividend = o1, divisor = 3, bound = intsqrt(o1), result = o2 }) ( where o1 "the odd part" of the number to decompose, o2 a list of its "contained" twos). This computes the primefactors of 2^120+1 in 17 minutes after all. But i cannot help feeling that this is an abuse of the State monad. The MonadFix has a function fix (a -> a) -> a and i added the first guard in divstep above for making this a fixpoint problem. For me the signature looks, as if using fix doesn't afford to create explicitely a variable of a MonadFix instance and a simple twoliner for divisions could do the job. What i do not understand at all from the documentation of fix is: "fix f is the least fixed point of the function f, i.e. the least defined x such that f x = x." What does "least" mean here ? There is nothing said about x being a variable of an instance of Ord. And why fix has not the type a -> (a -> a) -> a, means: How can i provide a starting point of the iteration x ==> f x ==> f (f x) ==> ...?
Am Dienstag, 18. Dezember 2007 17:26 schrieb Joost Behrends:
Hi,
since about three weeks i am learning Haskell now. One of my first excercises is to decompose an Integer into its primefactors. I already posted discussion on the solution to the problem 35 in "99 excercises".
My simple algorithm uses a datatype DivIter of 4 named fields together with the core iteration
But a simple recursion that returns the list of primefactors lazily would also solve the stack frame problem, wouldn't it? sort of factor 0 = error "0 has no factorisation" factor 1 = [] factor n | n < 0 = (-1):factor (-n) | even n = 2:factor (n `div` 2) | otherwise = oddFactors 3 n oddFactors k n | k*k > n = [n] | r == 0 = k:oddFactors k q | otherwise = oddFactors (k+2) n where (q,r) = n `divMod` k you can then start consuming the prime factors as they come, before the factorisation is complete.
divstep :: DivIter -> DivIter divstep x | divisor x > bound x = x
| ximod > 0 = x { divisor = (divisor x) +2 } | otherwise = x {dividend=xidiv,
bound=intsqrt(xidiv), result = result x ++ [divisor x] } where (xidiv, ximod) = divMod (dividend x) (divisor x)
(dividend x is already odd, when this is called).
The problem to solve for really large Integers is how to call divstep iterated without not accumulating billions of stack frames. Here is one possibility:
divisions = do y <- get if divisor y <= bound y then do put ( divstep y ) divisions else return y
(this for a version of divstep without the first guard) called from
res = execState divisions (DivIter { dividend = o1, divisor = 3, bound = intsqrt(o1), result = o2 })
( where o1 "the odd part" of the number to decompose, o2 a list of its "contained" twos). This computes the primefactors of 2^120+1 in 17 minutes after all. But i cannot help feeling that this is an abuse of the State monad. The MonadFix has a function fix (a -> a) -> a and i added the first guard in divstep above for making this a fixpoint problem.
For me the signature looks, as if using fix doesn't afford to create explicitely a variable of a MonadFix instance and a simple twoliner for divisions could do the job. What i do not understand at all from the documentation of fix is:
"fix f is the least fixed point of the function f, i.e. the least defined x such that f x = x."
What does "least" mean here ? There is nothing said about x being a variable of an instance of Ord. And why fix has not the type a -> (a -> a) -> a, means: How can i provide a starting point of the iteration x ==> f x ==> f (f x) ==> ...?
It's quite another thing, fix is not a fixed point iteration as used in calculus, least here means 'least defined'. The least defined of all values is 'undefined' (aka bottom, often denoted by _|_). For simple data types like Int or Bool, a value is either completely defined or undefined, and both True and False are more defined than bottom, but neither is more or less defined than the other. For a recursive data type like lists, you have a more interesting hierarrchy of definedness: _|_ is less defined than _|_:_|_ is less defined than _|_:_|_:_|_ is less defined than _|_:_|_:_|_:_|_ ... And definedness of list elements is also interesting, _|_:_|_:_|_ is less defined than 1:_|_:_|_ is less defined than 1:2:_|_ is less defined than 1:2:3:_|_ is less defined than ... You cannot use fix on a strict function (a function is strict iff f _|_ = _|_), as by its implementation, fix f = let x = fix f in f x IIRC, it's calculated without knowing what x is when f is called. fix f is basically lim xn, when x0 = undefined, x(n+1) = f xn And since f x is always at least as defined as x, xn is a monotonic sequence (monotonic with respect to definedness), so the limit exists - it's _|_ for strict functions, and if we follow a few steps of the easy example fix (1:), we see x0 = _|_ x1 = 1:_|_ x2 = 1:1:_|_ x3 = 1:1:1:_|_, hence fix (1:) == repeat 1. HTH, Daniel
Daniel Fischer
Am Dienstag, 18. Dezember 2007 17:26 schrieb Joost Behrends:
Hi,
since about three weeks i am learning Haskell now. One of my first excercises is to decompose an Integer into its primefactors. I already posted discussion on the solution to the problem 35 in "99 excercises".
My simple algorithm uses a datatype DivIter of 4 named fields together with the core iteration
But a simple recursion that returns the list of primefactors lazily would also solve the stack frame problem, wouldn't it? sort of factor 0 = error "0 has no factorisation" factor 1 = [] factor n | n < 0 = (-1):factor (-n) | even n = 2:factor (n `div` 2) | otherwise = oddFactors 3 n
oddFactors k n | k*k > n = [n] | r == 0 = k:oddFactors k q | otherwise = oddFactors (k+2) n where (q,r) = n `divMod` k
you can then start consuming the prime factors as they come, before the factorisation is complete.
Hi and thanks for your answers, @Daniel: no, this doesn't solve the stack problem. These are the primefactors of 2^120+1: [97,257,673,394783681,4278255361,46908728641]. oddFactors k n | otherwise = oddFactors (k+2) n could eventually push 394783681-673 function calls onto the stack before finding the factor 394783681. And indeed a recursive version of divisions trying to compute this ran more than two hours on my machine, before i stopped it (this is the time a python script needs for the computation). And there were peaks of memory use > 300 MB ! While the version with the State monad seems to work tail recursive - it has an absolutely constant memory use, slightly different per run, i watched 2044k and 2056k. And it takes around 17 minutes on my machine getting the result. Thus it's vital here, how the recursion is done. Because of that i separated divisions and divstep - for experimenting with divisions. If a factor is done, it should leave as little traces as possible on the machine. Both of your instructions for fix are well readable - thanks again. I'll spend some time studying them, but it seems, fix doesn't fix the problem.
Joost Behrends wrote:
@Daniel: no, this doesn't solve the stack problem. These are the primefactors of 2^120+1: [97,257,673,394783681,4278255361,46908728641].
oddFactors k n | otherwise = oddFactors (k+2) n
could eventually push 394783681-673 function calls onto the stack before finding the factor 394783681. And indeed a recursive version of divisions trying to compute this ran more than two hours on my machine, before i stopped it (this is the time a python script needs for the computation). And there were peaks of memory use > 300 MB ! While the version with the State monad seems to work tail recursive - it has an absolutely constant memory use, slightly different per run, i watched 2044k and 2056k. And it takes around 17 minutes on my machine getting the result.
Theoretically the recursions in oddFactors k n | otherwise = oddFactors (k+2) n and (*) divisions y |divisor y <= bound y = divisions (divstep y) do not cost stack space. They are tail recursions too! In general similar tail recursions do not cost stack space. Some programs using them use a lot of stack space, but that is because of non-strictness, not because of tail recursion itself. And such non-strictness is absent here due to the way the parameters have to be evaluated for all sorts of decisions before further recursive calls. Empirically my own experiments agree with the theory. I use GHC 6.8.1, ghc -O2, linux x86, P4 2 GHz, 1 GB RAM. First I just try 2^61+1 which is more tractable and still challenging enough. I think it's a great reference for calibration: it is significant enough as a test, and not too long to wait for. 2^61+1 = 3 * 768614336404564651 Your code on http://www.haskell.org/haskellwiki/?title=Talk:99_questions/31_to_41&oldid=17649 using (*) for divisions takes under 3 minutes and 1.4 MB memory throughout. Daniel's code takes under 6 minutes and 1.4 MB memory throughout. Your code using Control.Monad.State (which is Control.Monad.State.Lazy) for divisions takes under 11 minutes, and memory slowly grows from 1 MB to 6 MB. Your code using Control.Monad.State.Strict for divisions is similar to Control.Monad.State. Now I try 2^120+1. This probably is six times longer than 2^61+1, so I only try your code using (*) and Daniel's code. Your code using (*) takes under 19 minutes and 1.4 MB memory throughout. Daniel's code takes 27 minutes and 1.4 MB memory throughout. I have also tried turning off code optimizations and even deleting compiled code and running interpreted code in ghci for a while. I still see constant memory usage. My finding agrees with my theory and disagrees with your finding and probably your theory.
Albert Y. C. Lai
Theoretically the recursions in
oddFactors k n | otherwise = oddFactors (k+2) n
and
(*) divisions y |divisor y <= bound y = divisions (divstep y)
do not cost stack space. They are tail recursions too!
In general similar tail recursions do not cost stack space. Some programs using them use a lot of stack space, but that is because of non-strictness, not because of tail recursion itself. And such non-strictness is absent here due to the way the parameters have to be evaluated for all sorts of decisions before further recursive calls.
Thanks for all that benchwork (and especially your exponent 61). I must admit, that i ran the (*) version in ghci only (not having compiled it to a standalone version). Maybe ghci is configured wrongly on my system. As i remember, i tried (*) twice, coming near to memory exhaustion and not awaiting the result. I would really like a non-monadic version. What is interesting also is how near your 19 minutes came to my 17 (Windows XP, 2.2GHZ, 512MB). And the comparations to Daniels code seem to imply, that my named fields in DivIter are not very expensive, if at all. Is there documentation of tail recursion anywhere ? I searched (googling with site:www.haskell.org) and didn't find anything else than entries in mailing lists. Cheers, Joost
Am Dienstag, 18. Dezember 2007 schrieb Joost Behrends: <snip>
"fix f is the least fixed point of the function f, i.e. the least defined x such that f x = x."
What does "least" mean here ? There is nothing said about x being a variable of an instance of Ord. And why fix has not the type a -> (a -> a) -> a, means: How can i provide a starting point of the iteration x ==> f x ==> f (f x) ==> ...?
<snip> the starting point is "undefined". the ordering of functions is is_subset_of. a more detailed explanation: a function "A -> B" is a subset of the cartesian product "A x B", where for each element in A there is not more than one element in B. subsets are partially ordered. the empty set (the function "const undefined" or simply "undefined") is the lowest subset, and AxB is the largest (but in most cases not a function). the function "f0 (_::a) = (undefined::b)" is the lowest subset. the function "f1 ('x'::a) = (5+fromEnum 'x'::b)" is larger than f0. the function "f1' ('y'::a) = (7::b)" is larger than f0, and not equal (neither equal, nor lower, nor larger) to f1. the function "f2 (c::a) | isUpper c = (5 + fromEnum c::b)" is larger than f1, and not equal (neither equal, nor lower, nor larger) to f1'. the function "fn (c::a) | True = (5 + fromEnum c::b)" is one maximal defined function: it is defined on every input parameter. now, the "fix" function takes a function "construct_f::(a->a)->(a->a)" and calculates first "(construct_f undefined) :: (a->a)". "undefined :: (a->a)" equals f0, the lowest function/element, but it is not a fixpoint. "construct_f undefined" is a bit more "defined".... "construct_f . construct_f . construct_f . construct_f . ... (oo times) $ undefined" is the largest thing you can get this way, it does not need to be defined everywhere, but it is a fixpoint. there may be larger fixpoints, but no lower. example: fix construct_f where construct_f f = \x -> (if x==0 then 42 else f (x-1)) look at "construct_f undefined": it constructs a function which is defined on the input x==0; otherwise it tries to evaluate "undefined (x-1)", which is completely undefined. look at "construct_f $ construct_f undefined": it constructs a function which is defined on the input x==0 and x==1. "fix cf = cf (fix cf)" is the fixpoint function, and with this... "fix construct_f" constructs a function which is defined on all inputs x>=0, but not on inputs x<0. this function is one fixpoint (the least one) of construct_f. another function is a fixpoint of construct_f: "\x->42". but this is a larger function than the above fixpoint, so this is not the LEAST FIXPOINT; the above one is. you can test, whether it is a fixpoint: "construct_f (\x->42) == (\x->if x==0 then 42 else (\x->42)(x-1)) == (\x->if x==0 then 42 else 42) == (\x->42)" exercise1: "construct_f (\x->if x>=0 then 42 else 23) == ...?" exercise2: "construct_f (\x->if x>=0 then 42 else undefined) == ...?" another example: lists. fix (\fibs->1:1:zipWith (+) fibs (tail fibs)) i hope to have helped. - marc
Joost Behrends wrote:
since about three weeks i am learning Haskell now. One of my first exercises is to decompose an Integer into its primefactors.
How about separating the candidate prime numbers from the recursion factorize :: Integer -> [Integer] factorize = f primes' where primes' = 2:[3,5..] f (p:ps) n | r == 0 = p : f (p:ps) q | p*p > n = [n] | otherwise = f ps n where (q,r) = n `divMod` p For a faster factorization, just plug in a better primes' . Regards, apfelmus
apfelmus
How about separating the candidate prime numbers from the recursion
factorize :: Integer -> [Integer] factorize = f primes' where primes' = 2:[3,5..] f (p:ps) n | r == 0 = p : f (p:ps) q | p*p > n = [n] | otherwise = f ps n where (q,r) = n `divMod` p
Providing effectively primes' for that is simply impossible (besides: p < intsqrt n must stay, otherwise you have the expensive p*p at every step) talking about really big numbers as i did in my post. There are no fast generators iterating just through the primes firstly, and these lists get much too big also (for 2^120 you cannot even begin to use a list of the primes up to 2^(any appropriate x) ). What can be done is to iterate through odd numbers meeting as many primes as possible. We could do this: iterdivisors x | x == 0 = 3 | x == 1 = 5 | otherwise x = iterdivisors (x-1) + ((cycle [2,4]) !! x) This gives 7,11,13,17,19,23,25,29,31,35,37,41,43,47,49,53,55,59,61,63,67 ... i.e. exactly all primes and odds with greater primefactors as 3. We can improve that cycle avoiding the multiples of 5: ... | otherwise x = iterdivisors (x-1) + ((cycle [2,4,2,4,2,4,6,2,6] !! x) and we can do better by avoiding the multiples of 7 and so on (the length of these lists grows fast - it gets multiplied by every new avoided prime -, but we could provide that lists programmatically). And we must be sure, that cycle doesn't eat up memory for each new pass through the list. And we should use a more efficient representaion for the list of summands than a list. This is the first front. But the title of my post and much more interesting topic for learning Haskell is, how to avoid memory exhaustion by recursion. THIS was my intention and the reason why i erroneously brought MonadFix into the game. The recursion i described as follows
divisions = do y <- get if divisor y <= bound y then do put ( divstep y ) divisions else return y
makes a DESTRUCTIVE UPDATE of the DivIters (by put) and this kind of recursion seems not to "remember" itself (as i have understood, that is achieved by "tail recursion"). I just didn't like making DivIters to States. It's kind of lying code. However it worked and improved performance by a factor around 10 (or oo - perhaps a normal recursion exhausts 512MB memory for 2^120+1, as it will do for much smaller Integers, if they are prime) not to talk about footprint. Compiled for running standalone, it took 17 minutes, an equivalent python script 2 hours. This factor near 7 is not fully satisfactory. There are more workarounds beside the State monad for having destructive updates with Haskell. There are IORefs, there are updatable arrays and more. THAT is my question: Is this (only basically) the most efficient way to get them here ? Or is there still a way of getting a tail recursive Haskell function for iterating through the DivIters (outside any monads) ?? I would not get stroke dead by surprise if yes, but i couldn't find any. Cheers, Joost
On 12/20/07, Joost Behrends
makes a DESTRUCTIVE UPDATE of the DivIters (by put) and this kind of recursion seems not to "remember" itself (as i have understood, that is achieved by "tail recursion"). I just didn't like making DivIters to States. It's kind of lying code.
I just want to point out that this isn't true; "put" in the State monad doesn't do any destructive update at all (unlike, for example, IORef). You can tell this for yourself by looking at the type of "State s a" in Control.Monad.State: http://haskell.org/ghc/docs/latest/html/libraries/mtl/src/Control-Monad-Stat... newtype State s a = State { runState :: s -> (a,s) } That is, your "divisions" function of type divisions :: State DivIter DivIter is equivalent to the type runState divisions :: DivIter -> (DivIter, DivIter) and the code is the same as if you'd just passed the DivIter directly along. -- ryan
Ryan Ingram
I just want to point out that this isn't true; "put" in the State monad doesn't do any destructive update at all (unlike, for example, IORef). newtype State s a = State { runState :: s -> (a,s) } That is, your "divisions" function of type divisions :: State DivIter DivIter is equivalent to the type runState divisions :: DivIter -> (DivIter, DivIter) and the code is the same as if you'd just passed the DivIter directly along. -- ryan
Oops ? sry, i do not understand nothing. The syntax with the block in "newtype State s a = State { runState :: s -> (a,s) }" is not in the tutorials i read. And i do not understand " ... passed the DivIter directly along ". "Passed along" ?? But you'd probably better try not to explain this. Continued work in Haskell by itself will let me understand that some weeks from now. Thanks, Joost
On 12/20/07, Joost Behrends
The syntax with the block in
"newtype State s a = State { runState :: s -> (a,s) }"
is not in the tutorials i read.
newtype creates a new type which is treated exactly the same as an existing type at runtime, but which is distinct to the typechecker. For example, this code:
newtype X = MkX [Int] -- MkX :: [Int] -> X unX :: X -> [Int] unX (MkX val) = val
Now, you cannot call "sum (MkX [1,2,3])" because sum takes a list and you're passing an "X". But you can call "sum (unX (MkX [1,2,3]))". Since this is a newtype, the work done in "unX" and "MkX" is entirely erased at runtime; no allocation or pointer-following takes place. The syntax given is slightly more complicated:
newtype State s a = State { runState :: s -> (a,s) }
This is using the "record" syntax for data constructors, but it's basically the same as the following:
newtype State s a = State (s -> (a,s)) runState (State f) = f
So, runState is just unwrapping the function from the newtype. Since this is a newtype, there's no actual pointer traversal taking place; the distinction is only made during typechecking. The following two functions should compile to exactly the same machine code:
test1 :: Int -> ((), Int) test1 = \x -> ((), x+1)
test2 :: Int -> ((), Int) test2 = runState (State (\x -> ((), x+1)))
And i do not understand " ... passed the DivIter directly along ". "Passed along" ??
Recall your original code:
divisions :: State DivIter DivIter divisions = do y <- get if divisor y <= bound y then do put ( divstep y ) divisions else return y
Lets de-sugar that:
divisions = get >>= \y -> if divisor y <= bound y then (put (divstep y) >> divisions) else return y
Inlining get, put, and return:
divisions = (State (\s -> (s,s))) >>= \y -> if divisor y <= bound y then (State (\s -> ((), divstep y)) >> divisions) else State (\s -> (y, s))
After inlining (>>=) and (>>), and running a few simplification passes (in particular, runState (State x) = x, and beta-reduction), you end up with this:
divisions = State (\y -> if divisor y <= bound y then runState divisions (divstep y) else (y,y))
You can remove the State monad from this pretty trivially:
divisions :: DivIter -> (DivIter, DivIter) divisions y = if divisor y <= bound y then divisions (divstep y) else (y,y)
or,
divisions y | divisor y <= bound y = divisions (divstep y) | otherwise = (y,y)
This version of the code threads the "DivIter" state manually through each call, but it's exactly the same as what the State version is doing. No destructive updates anywhere to be seen! -- ryan
Joost Behrends wrote:
apfelmus writes:
How about separating the candidate prime numbers from the recursion
factorize :: Integer -> [Integer] factorize = f primes' where primes' = 2:[3,5..] f (p:ps) n | r == 0 = p : f (p:ps) q | p*p > n = [n] | otherwise = f ps n where (q,r) = n `divMod` p
(besides: p < intsqrt n must stay, otherwise you have the expensive p*p at every step)
Huh? p < intsqrt n is evaluated just as often as p*p > n , with changing n . Why would that be less expensive? Btw, the code above test for r==0 first, which means that the following p*p > n is tested exactly once for every prime candidate p .
Providing effectively primes' for that is simply impossible talking about really big numbers as i did in my post. There are no fast generators iterating just through the primes firstly
Sure, that's why I called it primes' . It's indented to be an easily computable list of prime candidates and as you write below, we can do better than using all odd numbers for that.
and these lists get much too big also (for 2^120 you cannot even begin to use a list of the primes up to 2^(any appropriate x) ).
Thanks to lazy evaluation, the list will be generated on demand and its elements are garbage collect once used. So, the code above will run in constant space. The list is more like a suspended loop.
What can be done is to iterate through odd numbers meeting as many primes as possible. We could do this:
iterdivisors x | x == 0 = 3 | x == 1 = 5 | otherwise x = iterdivisors (x-1) + ((cycle [2,4]) !! x)
This gives 7,11,13,17,19,23,25,29,31,35,37,41,43,47,49,53,55,59,61,63,67 ...
i.e. exactly all primes and odds with greater primefactors as 3. We can improve that cycle avoiding the multiples of 5:
... | otherwise x = iterdivisors (x-1) + ((cycle [2,4,2,4,2,4,6,2,6] !! x)
and we can do better by avoiding the multiples of 7 and so on
(the length of these lists grows fast - it gets multiplied by every new avoided prime -, but we could provide that lists programmatically). And we must be sure, that cycle doesn't eat up memory for each new pass through the list. And we should use a more efficient representaion for the list of summands than a list.
Huh, this looks very expensive. I'd avoid indices like x altogether and use a plain list instead, we don't need random access to the prime candidates, after all.
But the title of my post and much more interesting topic for learning Haskell is, how to avoid memory exhaustion by recursion.
Maybe you stumbled over common examples for a stack overflow like foldr (+) 0 foldl (+) 0 whereas foldl' (+) 0 runs without. See also http://en.wikibooks.org/wiki/Haskell/Performance_Introduction http://www.haskell.org/haskellwiki/Stack_overflow
THIS was my intention and the reason why i erroneously brought MonadFix into the game. The recursion i described as follows
divisions = do y <- get if divisor y <= bound y then do put ( divstep y ) divisions else return y
makes a DESTRUCTIVE UPDATE of the DivIters (by put)
Huh? The State monad doesn't do destructive updates, to the contrary. (neither do IORefs or STRefs, only STArrays or something do).
and this kind of recursion seems not to "remember" itself (as i have understood, that is achieved by "tail recursion"). I just didn't like making DivIters to States. It's kind of lying code.
Because of lazy evaluation, tail recursion is not what it seems to be in Haskell.
However it worked and improved performance by a factor around 10 (or oo - perhaps a normal recursion exhausts 512MB memory for 2^120+1, as it will do for much smaller Integers, if they are prime) not to talk about footprint. Compiled for running standalone, it took 17 minutes, an equivalent python script 2 hours. This factor near 7 is not fully satisfactory.
Using the State monad introduces unnecessary overhead. Also, I assume that you ran the compiler with the -O2 flag to enable optimizations?
Or is there still a way of getting a tail recursive Haskell function for iterating through the DivIters (outside any monads) ?? I would not get stroke dead by surprise if yes, but i couldn't find any.
I don't understand why you use a complicated DivIters data structure. Passing two simple parameters does the job just fine. Regards, apfelmus
@apfelmus, please read my code. I introduced DivIter to separate divstep from divisions. But it stores intsqrt dividend also. Thus the sqrt is only recomputed, when a new factor is found. Concerning primes': With the sieve of Eratosthenes we cannot make a lazy list, we need the whole list at any point of computation. And beyond the end we fall back to x -> x+2. This is the difference to the list of summands i proposed. And we have the arbitrary choice of [2,4], [2,4,2,4,2,4,2,6,2,6] or any longer. Concerning State: There is useful instruction by Ryan Ingram too here. I will study that and yours. Thank you ! Cheers, Joost
Am Freitag, 21. Dezember 2007 11:33 schrieb apfelmus:
Joost Behrends wrote:
apfelmus writes:
How about separating the candidate prime numbers from the recursion
factorize :: Integer -> [Integer] factorize = f primes' where primes' = 2:[3,5..] f (p:ps) n
| r == 0 = p : f (p:ps) q | p*p > n = [n] | otherwise = f ps n
where (q,r) = n `divMod` p
(besides: p < intsqrt n must stay, otherwise you have the expensive p*p at every step)
Huh? p < intsqrt n is evaluated just as often as p*p > n , with changing n . Why would that be less expensive? Btw, the code above test for r==0 first, which means that the following p*p > n is tested exactly once for every prime candidate p .
However, when you do the sensible thing (which Joost did) and have the intsqrt a parameter of the function, like in factorize :: Integer -> [Integer] factorize n = f n (intsqrt n) primes' where primes' = something more or less clever here f m sr (p:ps) | r == 0 = p:f q (intsqrt q) (p:ps) | p > sr = if m == 1 then [] else [m] | otherwise = f m sr ps where (q,r) = m `quotRem` p , then you only have the expensive intsqrt for each prime factor, and the test for each candidate is only one comparison instead of one multiplication and one comparison. Since usually the number of prime factors is minuscule compared to the number of candidates to be tested, that's indeed faster for most inputs (in my tests ~28.5% less run time). Some speed is also gained by using quotRem instead of divMod (pace, Henning, I agree divMod is preferable, but quotRem is a primitive). But of course, the most could be gained from a better algorithm.
Regards, apfelmus
Cheers, Daniel
Daniel Fischer wrote:
apfelmus writes:
| r == 0 = p : f (p:ps) q | p*p > n = [n] | otherwise = f ps n
However, when you do the sensible thing (which Joost did) and have the intsqrt a parameter of the function, like in
factorize :: Integer -> [Integer] factorize n = f n (intsqrt n) primes' where primes' = something more or less clever here f m sr (p:ps) | r == 0 = p:f q (intsqrt q) (p:ps) | p > sr = if m == 1 then [] else [m] | otherwise = f m sr ps where (q,r) = m `quotRem` p
, then you only have the expensive intsqrt for each prime factor, and the test for each candidate is only one comparison instead of one multiplication and one comparison.
Ah thanks, sorry for not seeing it earlier. My thinking was that intsqrt q is calculated multiple times (for q , q/p, q/p^2, ...) per prime candidate p when n is divisible by a large power of p . But I didn't see that, in return, intsqrt q stays the same for candidates that aren't factors of n . Compared to that, p*p is only calculated once per candidate, but then for every candidate. The former is clearly preferable to the latter. In fact, thanks to lazy evaluation, the above code (test r==0 before p > sr) evaluates intsqrt q at most once per prime candidate and thus combines both advantages. Regards, apfelmus
@apfelmus, please read my code. I introduced DivIter to separate divstep from divisions. But it stores intsqrt dividend also. Thus the sqrt is only recomputed, when a new factor is found. Concerning primes': With the sieve of Eratosthenes we cannot make a lazy list, we need the whole list at any point of computation. And beyond the end we fall back to x -> x+2. This is the difference to the list of summands i proposed. And we have the arbitrary choice of [2,4], [2,4,2,4,2,4,2,6,2,6] or any longer. Concerning State: There is useful instruction by Ryan Ingram too here. I will study that and yours. Thank you ! Cheers, Joost
apfelmus
Huh? p < intsqrt n is evaluated just as often as p*p > n , with changing n . Why would that be less expensive? Btw, the code above test for r==0 first, which means that the following p*p > n is tested exactly once for every prime candidate p .
No. One point in the introduction of DivIter is, that intsqrt dividend is stored there and only recomputed, when a new factor is found. And concerning my cycled lists of summands as [2,4] or [2,4,2,4,2,4,2,6,2,6]: There simply is no function easily yielding primes for your list primes'. If we use the sieve of Eratosthenes, we must have the whole list of found primes up to a certain point in memory for proceeding beyond that certain point. We cannot gain anything by lazy evaluation. Not with the sieve of Eratosthenes - and there is no other reliable mechanism. What is more - if we have a list of primes, possibly up to 1,000,000 - what shall we do for efficiently yielding divisors beyond 1,000,000 ? We would have to fall back to x -> x+2. Thus an easily computable function stepping through all primes can only be a function, which yields primes plus some more odd numbers. This is, what i tried. Alternating addition of 2 and 4 to the current divisor can be continued beyond any bound. And i am not forced to use any of the longer list of summands - indeed, which of these lists to choose should be adapted to the size of the number to decompose. Concerning the State monad: Yes, i was wrong here completely. Ryan Ingram gave detailed instructions why. And Albert Y.C. Lai pointed out, that normal recursion for division works tail recursive too. It didn't on my computer - perhaps i misconfigured ghci (where i tried it). Let us close discussion of this point. Cheers, Joost
Joost Behrends wrote:
apfelmus writes:
Huh? p < intsqrt n is evaluated just as often as p*p > n , with changing n . Why would that be less expensive? Btw, the code above test for r==0 first, which means that the following p*p > n is tested exactly once for every prime candidate p .
No. One point in the introduction of DivIter is, that intsqrt dividend is stored there and only recomputed, when a new factor is found.
Yes, I'm sorry, it didn't occur to me that recomputing per actual prime factor (with multiplicities, i.e. p^5 counted 5 times) is better than recomputing per candidate factor (without multiplicities, i.e. p^5 counted only once).
And concerning my cycled lists of summands as [2,4] or [2,4,2,4,2,4,2,6,2,6]:
an easily computable function stepping through all primes can only be a function, which yields primes plus some more odd numbers. This is, what i tried.
Yes, this scheme was my intention, too. The list primes' doesn't need to (and indeed shouldn't) be a list of actual primes, just a good guess like primes' = 2:[3,5] primes' = 2:3:scanl (+) 5 (cycle [2,4]) or something with [2,4,2,4,2,4,2,6,2,6]. So, it's an infinite list of numbers that qualify as candidates for being a prime factor of n (which I called "prime candidates". Not a good name, since they don't need to be actual prime numbers.) What I want to say is that using such an infinite is a nice way to separate the generate-prime-factor-candidates-logic from the test-all-those-candidates-loop. It's not necessary to hard-code it with a predicate like
iterdivisors x | x == 0 = 3 | x == 1 = 5 | otherwise x = iterdivisors (x-1) + ((cycle [2,4]) !! x)
(which, in this particular form, is hopelessly inefficient) or special step functions like
d2 :: DivIter -> DivIter d2 x |dividend x `mod` divisor x > 0 = x { divisor = divisor x + 2} |otherwise = found x d4 :: DivIter -> DivIter d4 x |dividend x `mod` divisor x > 0 = x { divisor = divisor x + 4} |otherwise = found x d6 :: DivIter -> DivIter d6 x |dividend x `mod` divisor x > 0 = x { divisor = divisor x + 6} |otherwise = found x
divisions :: DivIter -> DivIter divisions y |or[divisor y == 3, divisor y == 5] = divisions (d2 y) |divisor y <= bound y = divisions (d6$d2$d6$d4$d2$d4$d2$d4 y) |otherwise = y
It's not even necessary to treat 2 in a special way like
twosect :: (Integer,[Integer]) -> (Integer,[Integer]) twosect m |odd (fst m) = m |even (fst m) = twosect (div (fst m) 2, snd m ++ [2])
does, the primes' list handles it all. Regards apfelmus
participants (7)
-
Albert Y. C. Lai -
apfelmus -
Daniel Fischer -
Joost Behrends -
Joost Behrends -
Marc A. Ziegert -
Ryan Ingram