Looking for suggestions to improve my algorithm

Hello Haskellers, I have been trying to learn a bit about Haskell by solving Project Euler problems. For those of you who don't know what Project Euler is, see http://projecteuler.net After solving problem 21, which is related to amicable pairs, I decided to attempt problem 95 since it is an extension of the same topic. The full description of problem 95 is here: http://projecteuler.net/index.php?section=problems&id=95 This is the summary: "Find the smallest member of the longest amicable chain with no element exceeding one million." I have attempted to solve this problem, but my solution is too resource hungry to search the full set of 1000000 numbers. I am hoping someone reading this list can suggest : - How I can improve my algorithm - An alternative algorithm that will be more efficient - Ways to profile a Haskell program to help me understand why my implementation is performing poorly. In addition to the question above, I would also be interested in comments on the style of the code. If there is a more idiomatic way to write portions of the code, I would like to see what it is. My solution is at the bottom of this e-mail. The program will probably run obscenely slow or die from stack overflow unless you replace [1..999999] with [1..somethingSmaller] in main. Thanks, David Frey --- BEGIN Main.hs --- module Main where import ProjectEuler (takeUntil, divisors) import qualified Data.Map as M import qualified Data.IntSet as I main = let initialContext = Context (I.fromList []) 0 0 in print $ cycleStart $ foldl checkForChain initialContext [1..999999] {- Idea: * Put all the numbers that have been visited into a map regardless of whether they are a part of a chain or not. * Store the min element in the cycle and the number of elements in the cycle * As you process, from 1->n if the stopping conditions for a sumOfDivisors result are: * has already been seen before * number is less than the start of this chain attempt * >= 1,000,000 -} data Context = Context {seenNum::I.IntSet, cycleStart::Int, cycleLength::Int} hasBeenSeen :: Int -> Context -> Bool hasBeenSeen n context = I.member n (seenNum context) markSeen :: Int -> Context -> Context markSeen n context = context { seenNum = (I.insert n (seenNum context)) } deleteFromSeen :: Int -> Context -> Context deleteFromSeen n context = context { seenNum = (I.delete n (seenNum context)) } {- - Examines the context to see if the input has potential to be a chain or not - based on whether the input number has been visited before. If it could be a - chain, an attempt is made to build the chain. -} checkForChain :: Context -> Int -> Context checkForChain context n = if hasBeenSeen n context then deleteFromSeen n context else buildChain context (sum $ divisors n) n [n] {- - Builds a chain until ones of the 3 stopping conditions are met or a chain is - found. If a chain is found the context will be updated with the new chain if - it is the longest. - - Stopping Conditions: - * Number has already been seen before - * Number is less than the start of this chain attempt - * Number >= 1,000,000 -} buildChain :: Context -> Int -> Int -> [Int] -> Context buildChain context candidate first cycleList = if elem candidate cycleList then foundChain (takeUntil ((==) candidate) cycleList) context else if candidate < first || candidate >= 1000000 || hasBeenSeen candidate context then context else buildChain (markSeen candidate context) (sum $ divisors candidate) first (candidate : cycleList) {- - Updates the context with the new longest chain and the start value if the - chain input parameter is longer than the one currently in the context. -} foundChain :: [Int] -> Context -> Context foundChain ls context = let l = length ls m = minimum ls in if l > (cycleLength context) then context { cycleLength = l, cycleStart = m } else if l == (cycleLength context) then let m = minimum ls in if m < (cycleStart context) then context { cycleStart = m } else context else context --- END Main.hs --- I put a bunch of common functions in ProjectEuler.hs. Here are the relevant functions for this problem: {- - Gets all of the proper divisors of a number. That is all divisors starting - from 1, but not including itself. -} divisors :: (Integral a) => a -> [a] divisors n = let p1 = [x | x <- [1 .. floor $ sqrt $ fromIntegral n], n `mod` x == 0] p2 = map (div n) (tail p1) in p1 `concatNoDupe` (reverse p2) where {- - Concatenate two lists. If the last element in the first list and the - first element in the second list are ==, produce only the value from - the first list in the output. -} concatNoDupe :: (Eq a) => [a] -> [a] -> [a] concatNoDupe [] ys = ys concatNoDupe [x] (y:ys) = if x == y then (y : ys) else (x : y : ys) concatNoDupe (x:xs) ys = x : (concatNoDupe xs ys) {- - Similar to takeWhile, but also takes the first element that fails the - predicate. -} takeUntil :: (a -> Bool) -> [a] -> [a] takeUntil pred (x:xs) = (x : if pred x then [] else takeUntil pred xs) takeUntil _ [] = []

On 8/29/07, David Frey
Hello Haskellers,
I have been trying to learn a bit about Haskell by solving Project Euler problems. For those of you who don't know what Project Euler is, see http://projecteuler.net
After solving problem 21, which is related to amicable pairs, I decided to attempt problem 95 since it is an extension of the same topic.
The full description of problem 95 is here: http://projecteuler.net/index.php?section=problems&id=95
This is the summary: "Find the smallest member of the longest amicable chain with no element exceeding one million."
I have attempted to solve this problem, but my solution is too resource hungry to search the full set of 1000000 numbers.
I am hoping someone reading this list can suggest : - How I can improve my algorithm - An alternative algorithm that will be more efficient - Ways to profile a Haskell program to help me understand why my implementation is performing poorly.
In addition to the question above, I would also be interested in comments on the style of the code. If there is a more idiomatic way to write portions of the code, I would like to see what it is.
My solution is at the bottom of this e-mail. The program will probably run obscenely slow or die from stack overflow unless you replace [1..999999] with [1..somethingSmaller] in main.
Thanks, David Frey
Hi David, Project Euler is a good way to learn some Haskell, although it does tend to give you a crash course in understanding laziness, efficiency and such in Haskell (whether that's good or bad depends on your point of view!). All you need to do to stop your program from overflowing the stack is to change foldl to foldl' in the definition of main (foldl' is in Data.List so you'll also have to import that). The problem with foldl is that since Haskell is lazy, instead of evaluating things as it goes along, it builds up a huge string of thunks (=unevaluated expressions) and doesn't even start evaluating anything until it reaches the end of the list, which can easily blow the stack. foldl' is a strict version of foldl which forces evaluation to take place as it goes along. For an excellent explanation of all of this, I suggest you read http://haskell.org/haskellwiki/Stack_overflow. As soon as I replaced foldl with foldl' in your code, it no longer overflows the stack -- however, it's still quite slow! I didn't wait long enough for it to finish when I ran it up to a million. I don't have any advice on that front at the moment, perhaps someone else will come along with a suggestion or two. In the meantime, you can poke around http://haskell.org/haskellwiki/Performance to see if it gives you any ideas. Hope this helps, -Brent

Am Donnerstag, 30. August 2007 01:08 schrieb Brent Yorgey:
Hi David,
Project Euler is a good way to learn some Haskell, although it does tend to give you a crash course in understanding laziness, efficiency and such in Haskell (whether that's good or bad depends on your point of view!).
All you need to do to stop your program from overflowing the stack is to change foldl to foldl' in the definition of main (foldl' is in Data.List so you'll also have to import that).
I didn't even need to do that. Worked fine for me as it is (I compiled with -O2, that probably helped).
The problem with foldl is that since Haskell is lazy, instead of evaluating things as it goes along, it builds up a huge string of thunks (=unevaluated expressions) and doesn't even start evaluating anything until it reaches the end of the list, which can easily blow the stack. foldl' is a strict version of foldl which forces evaluation to take place as it goes along. For an excellent explanation of all of this, I suggest you read http://haskell.org/haskellwiki/Stack_overflow.
As soon as I replaced foldl with foldl' in your code, it no longer overflows the stack -- however, it's still quite slow! I didn't wait long enough for it to finish when I ran it up to a million.
45 seconds on my 1200MHz thing, not incredibly fast, but not obscenely slow either.
I don't have any advice on that front at the moment, perhaps someone else will come along with a suggestion or two. In the meantime, you can poke around http://haskell.org/haskellwiki/Performance to see if it gives you any ideas.
Hope this helps, -Brent
Cheers, Daniel

Am Mittwoch, 29. August 2007 23:12 schrieb David Frey:
Hello Haskellers,
I have been trying to learn a bit about Haskell by solving Project Euler problems.
Good :)
For those of you who don't know what Project Euler is, see http://projecteuler.net
After solving problem 21, which is related to amicable pairs, I decided to attempt problem 95 since it is an extension of the same topic.
The full description of problem 95 is here: http://projecteuler.net/index.php?section=problems&id=95
This is the summary: "Find the smallest member of the longest amicable chain with no element exceeding one million."
I have attempted to solve this problem, but my solution is too resource hungry to search the full set of 1000000 numbers.
I haven't looked closely, but what stands out is the obscenely slow divisors function. This could be sped up if you factored the numbers using the list of primes less than 1000. But since you need the factorization of all numbers <= 1000000 (that ought to be included, the wording is 'not exceeding' - but it doesn't appear in the chain, so you won't get a false answer by omitting it), you can do even better. My suggestion: 1. create an array or a map containing the sum of divisors of each number from 1 to 1000000 1.1. 'smallest prime factor' is useful for that 2. look for chains in that array/map
I am hoping someone reading this list can suggest : - How I can improve my algorithm - An alternative algorithm that will be more efficient - Ways to profile a Haskell program to help me understand why my implementation is performing poorly.
compile with -prof -auto-all and run with euler95 +RTS -P then read euler95.prof and take a look at the ghc users guide for more profiling options
In addition to the question above, I would also be interested in comments on the style of the code. If there is a more idiomatic way to write portions of the code, I would like to see what it is.
I'll take a look and see. But don't hold your breath, 'idiomatic' isn't one of my strengths.
My solution is at the bottom of this e-mail. The program will probably run obscenely slow or die from stack overflow unless you replace [1..999999] with [1..somethingSmaller] in main.
Thanks, David Frey
Cheers, Daniel

The algorithm I would use: Treat the problem as a directed graph of 1 million vertices where each vertex has at most 1 outgoing edge. Each vertex is labeled with a number N(v). For each vertex v, if sum(divisors(N(v))) <= 1 million, that vertex has an outgoing edge to vertex v', where N(v') == sum(divisors(N(v))). Call this graph G0. 1) No vertex with 0 incoming edges can be part of an amicable chain. 2) No vertex with 0 outgoing edges can be part of an amicable chain. If there are any such vertices in the graph, you can remove them. Call this new graph G1. Lemma 3) G1 is a directed graph where every vertex has at most one outgoing edge. Proof) G0 has this property and we have only deleted vertices from G0, which can only remove the You can repeat this process until no such vertices remain. Call the final resultant graph Gn. Lemma 4) Each vertex in Gn has exactly 1 incoming edge and exactly 1 outgoing edge. Proof: 1) No vertex has 0 incoming edges or else we would have removed it above. 2) Pidgeonhole principle: There are exactly as many edges as vertices. If any vertex had 2 incoming edges, there must be at least one vertex with 0 incoming vertices. Leamma 5) Each vertex in Gn is part of exactly one cycle. Proof) Follow the outgoing vertex chain from a vertex. Since Gn is finite and each vertex has exactly 1 incoming & 1 outgoing edge, eventually you will return to that vertex. This forms a cycle and there were no other possible edges to follow. Therefore: Algorithm: 1) Generate G0 2) Iteratively remove vertices from G0 to create G1, G2, etc. until there are no vertices with 0 incoming or outgoing edges. Let this graph be Gn 3) Split the vertices of Gn into equivalence classes where each class represents a cycle 4) Find the largest equivalence class. 5) Find the smallest label in that equivalence class. As to how to profile:
ghc --make -auto-all -prof -O2 main.hs main +RTS -p will generate a profile report in main.prof. I've found that functions that do a large amount of allocation tend to be the easiest to optimize. Adding additional strictness annotations (seq) sometimes helps.
For constructing graphs, you may want to look at http://www.haskell.org/haskellwiki/Tying_the_Knot as an alternative to the "IsSeen" maps that you might use to iteratively construct a graph. I'd probably use something like this: type Graph a = Map a Vertex data Vertex a = Vertex a [Vertex a] instance Eq a => Eq (Vertex a) where Vertex a _ == Vertex b _ = a == b then: constructGraph :: Ord a => (a -> Graph a -> [Vertex a]) -> [a] -> Graph a constructGraph mkEdges labels = graph where graph = Map.fromList [(label, Vertex label (mkEdges label graph)) | label <- labels] Notice that graph is used in its own definition! You can picture this as first creating a map like this: 1 -> Vertex 1 _ 2 -> Vertex 2 _ 3 -> Vertex 3 _ ... Then when you query the map for the edges of Vertex 2, only then does the mkEdges get called, with the constructed graph as input, updating the "_" with pointers to the actual vertices. -- ryan

I managed it in 7 seconds (on 1500 MHz) with an idea close to yours (but I used IntSet, not IntMap), Daniel Fisher gave you some good ideas to achieve it, the real snail in this problem is the sumDivisors function. -- Jedaï

2007/8/30, Chaddaï Fouché
I managed it in 7 seconds (on 1500 MHz) with an idea close to yours (but I used IntSet, not IntMap), Daniel Fisher gave you some good ideas to achieve it, the real snail in this problem is the sumDivisors function.
I put my final solution on the wiki, it get it done in 6s now (on a Pentium M 1.73Mhz). -- Jedaï

I played around with this for a while based on the same sort of algorithm and ended up with a similar solution too. It turns out the operations saved by keeping track of already visited nodes are more than outweighed by the cost of doing so. (As you can see, I still have the hook in my code where amCn returns the visited list that I'm just disregarding. Ii had initially kept a giant array of all values to calculate, but the cost of that was unmanageable. After that, my big roadblock was that I couldn't come up with a good sumDivisors function to save my life. I tried a number of "optimized" methods that combined trial division with reduction based on prime factorizations, but i had to either reduce the lists by checking divisibility again somewhere along the way, or calling nub, and strictness and memoization still didn't let me produce a factorization in reasonable time. In the end, I ended up lifting yours. The only problem is that I've been staring at it for a bit and am not really sure how it works. I'd love an explanation. In any case, the code of my solution follows: amCn maxval n = amCn' (propDSum n) [] where amCn' cur visitedlist = if cur > maxval || cur < n then (0,visitedlist) else if elem cur visitedlist then (0,visitedlist) else if (cur == n) then ((length visitedlist) + 1,visitedlist) else (amCn' $! (propDSum cur)) $! (cur:visitedlist) longestAmTo maxval = longestAm' 2 (0,0) where longestAm' n bestFit@(chainLen,minVal) = if n > maxval then bestFit else longestAm' (n+1) $! bestFit' where (count, visited) = amCn maxval n bestFit' = if chainLen > count then bestFit else (count,n) properDivisorsSum :: UArray Int Int properDivisorsSum = accumArray (+) 1 (0,1000000) $ (0,-1):[(k,factor)| factor<-[2..1000000 `div` 2] , k<-[2*factor,2*factor+factor..1000000] ] propDSum n = properDivisorsSum ! n --S On Aug 30, 2007, at 11:33 AM, Chaddaï Fouché wrote:
2007/8/30, Chaddaï Fouché
: I managed it in 7 seconds (on 1500 MHz) with an idea close to yours (but I used IntSet, not IntMap), Daniel Fisher gave you some good ideas to achieve it, the real snail in this problem is the sumDivisors function.
I put my final solution on the wiki, it get it done in 6s now (on a Pentium M 1.73Mhz).
-- Jedaï _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe

Right, your program is 2 times faster than mine on my machine... I wonder if there is a better structure to do this bookkeeping than IntSet (maybe Sequence slightly remanied ?), anyway it goes to show how sometimes the bookkeeping can be more expensive than the operations it's meant to prevent ! As for my sumOfProperDivisors function it's dead simple, I would even say stupid (but it's faster than anything else I tried for now). An accumArray use the function you give it to update a cell, ok ? Here it's just (+) and all cells began their life as 1 since 1 is a proper divisors of all numbers (except 1, thus the (1,-1) for correctness). The following list just associate each proper divisor of a num with it, so the final value of a cell is the sum of those proper divisors. To achieve that we make factor take the value of all proper divisors possible for numbers from 1 to 1000000, in other words [2..1000000 `div` 2] (the `div` 2 is ok since we're speaking about proper divisors here), and then we go on to associate this divisor with all the numbers he divides properly, which are [factor * 2, factor * 3,...1000000]. Is it clear now ? -- Jedaï
participants (6)
-
Brent Yorgey
-
Chaddaï Fouché
-
Daniel Fischer
-
David Frey
-
Ryan Ingram
-
Sterling Clover