List comprehension takes too much memory

Question

I'm a beginner to Haskell and used it to solve some 50 problems of Project Euler but now I'm stuck at problem 66. The problem is that the compiled code (ghc -O2 --make problem66.hs) takes all my machine's free memory after 10-20 seconds. My code looks like this:

-- Project Euler, problem 66

diophantine x y d = x^2 - d*y^2 == 1

minimalsolution d = take 1 [(x, y, d) | y <- [2..],
                            let x = round $ sqrt $ fromIntegral (d*y^2+1),
                            diophantine x y d]

issquare x = (round $ sqrt $ fromIntegral x)^2 == x

main = do
    print (map minimalsolution (filter (not . issquare) [1..1000]))

I have a hunch that the problem lies in the infinite list inside the list comprehension for minimalsolution.

I actually thought that due to lazyness, Haskell would evaluate the list only until it finds one element (because of take 1) and on the way discard everything for which diophantine evaluates to False. Am I wrong there?

Interestingly, I did not see this behaviour in ghci. Is it because processing inside ghci is so much slower that I just would have to wait until I see the memory consumption explode - or is it something else?

No spoilers, please. All I want to know is where the extreme memory consumption comes from and how I can fix it.

Answer 1

I haven't profiled before, so stone throwers are welcome.

Haskell determines that [2..] is a constant and is reused for every element of the list, despite take 1 using only one element of that list; so it memoizes the list for computing future elements of the same list. You get stuck computing value for d=61.

---

Edit:

What's interesting, this one terminates for [1..1000]:

minimalsolution d = take 1 [(x, y, d) | y <- [2..] :: [Int],
                            let x = round $ sqrt $ fromIntegral (d*y^2+1),
                            diophantine x y d]

Just added :: [Int]. Memory use looks stable at 1MB. Using Int64 reproduces the problem.

minimalsolution d = take 1 [(x, y, d) | y <- [2..] :: [Int64],
                            let x = round $ sqrt $ fromIntegral (d*y^2+1),
                            diophantine x y d]

---

Edit:

Well, as has been suggested, the difference is caused by overflow. The solution to d=61 is reported as (5983,20568,61), but 5983^2 is nowhere near 61*20568^2.

Answer 2

Inside of the comprehension creating unnecessary Double instances on each value of y.

I couldn't find a solution using list comprehensions that didn't have the space blowup. But rewriting using recursion yields a stable memory profile.

diophantine :: Int -> Int -> Int -> Bool
diophantine x y d = x^2 - d*y^2 == 1

minimalsolution ::  Int -> (Int, Int, Int)
minimalsolution d = go 2
  where
    d0 = fromIntegral d
    go a =
      let y = fromIntegral a
          x = round $ sqrt $ (d0*y^2+1) in
      if diophantine x y d then
        (x, y, d)
      else
        go (y+1)