Haskell：如何在Atkin的Sieve中找到方程的整数解的数量？

Question

我目前正在尝试在Haskell中实现Atkin的Sieve

在Wikipedia article on the Sieve of Atkin的第3步中，我需要找到多个方程的整数解的数量。

然而，我对这些方程中的第一个（4x²+y²= n，x> 0，y> 0，其中n是正整数列表中的条目）的解决方案在使用任何n的查询时产生无限循环。

这是我到目前为止这部分问题的代码：

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))))

它被WinGHCi加载得很好，但当我查询时，例如eq1 0它只是停留在无限循环中，在产生答案之前必须被打断。我怀疑它在x和y的两个任务之间循环。

我怎么能阻止这个？这甚至可能吗？

编辑：实现无限循环必须在哪里。

Answer 1

我将首先重新格式化您的代码，使其更具可读性。换行有用！此外，操作顺序可以减少括号的重量。边注：

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

也可以写

f x | e1
    , e2
    , e3
    = e4

这对眼睛来说可能更容易。

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)

现在我可以立即看到逻辑是不稳定的。鉴于n，你计算x和y。然后，您可以递归停止或调用该函数。但是，在递归调用中，您可以保证停止！所以，即使你是对的，你肯定会有一个语义问题，总是返回0或1。

但正如你所见，这不是唯一的问题。根据x，你也用y和y来定义x。现在有一些重要的情况，这种相互递归是有用的。但是当相互递归的值是像整数一样的“原子”时，你肯定会得到一个无限循环。哈斯克尔不会为你解决方程式;那是你的工作！

---

这是我的建议：

从强力列表理解解决方案开始：

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

接下来，您可以使用近似整数平方根来缩小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