0

I am currently trying to implement the Sieve of Atkin in Haskell

In step 3 on the Wikipedia article on the Sieve of Atkin I need to find the number of Integer solutions to multiple equations.

However my solution to the first of these equations (4x² + y² = n, x > 0, y > 0 with n being a entry in a list of positive Integers) produces an infinite loop upon a query with any n.

This is my code for this part of the problem so far:

eq1 :: Integer -> Integer
eq1 n = eq1_ n []
eq1_ :: Integer -> [(Integer, Integer)] -> Integer
eq1_ n list     | (x > 0) && (y > 0) && (n == 4*(x^2) + (y^2)) && (notElem ((x,y)) list) = eq1_ n ([(x, y)] ++ list)
                | otherwise = toInteger (length list)
                where
                    x = floor (sqrt (fromIntegral ((n - y^2) `div` 4)))
                    y = floor (sqrt (fromIntegral (n - 4*(x^2))))

It is loaded just fine by WinGHCi, but when I query e.g. eq1 0 it just stays in an infinite loop and has to be interrupted before producing an answer. I suspect it goes in a loop between the two assignments of x and y.

How can I prevent this? Is this even possible?

Edit: Realised where the infinite loop must be.

D'Deletham i-Ra'tleihfi
  • 1
  • 1

1 Answers1

1

I'm going to start by reformatting your code a tad to make it more readable. Line breaks are helpful! Also, the order of operations can reduce the weight of parentheses. Side note:

f x | e1 && e2 && e3 = e4

can also be written

f x | e1
    , e2
    , e3
    = e4

which may be easier on the eyes.

eq1 :: Integer -> Integer
eq1 n = eq1_ n []

eq1_ :: Integer -> [(Integer, Integer)] -> Integer
eq1_ n list
  | x > 0 &&
    y > 0 &&
    n == 4*x^2 + y^2 &&
    notElem (x,y) list
  = eq1_ n ([(x, y)] ++ list)
  | otherwise
  = toInteger (length list)
  where
    isqrt = floor . sqrt . fromIntegral
    x = isqrt $ (n - y^2) `div` 4
    y = isqrt $ n - 4*(x^2)

Now I can immediately see that the logic is wonky. Given n, you calculate x and y. Then you either stop or call the function recursively. On the recursive call, however, you're guaranteed to stop! So even if you were otherwise right, you'd definitely have a semantic problem, always returning 0 or 1.

But as you've seen, that's not the only problem. You're also defining x in terms of y and y in terms of x. Now there are important situations where such mutual recursion is useful. But when the mutually recursive values are "atomic" things like integers, you're sure to get an infinite loop. Haskell won't solve the equations for you; that's your job!

Here's my suggestion:

Start with a brute force list comprehension solution:

sols n
  = [(x,y)
    |x <- takeWhile (\p -> 4 * p^2 < n) [1..]
    ,y <- takeWhile (\q -> f x y <= n) [1..]
    ,f x y = n]
  where
    f x y = 4*x^2+y^2

Next, you can use an approximate integer square root to narrow the search space for y:

sols n
  = [(x,y)
    |x <- takeWhile (\p -> 4 * p^2 < n) [1..]
    ,y <- takeWhile
            (\q -> f x y <= n)
          [floor(sqrt(fromIntegral(n-4*x^2)))..]
    ,f x y = n]
  where
    f x y = 4*x^2+y^2
dfeuer
  • 48,079
  • 5
  • 63
  • 167