如何使用 Haskell（有效地）找到数字的最大素因数？

Question

我正在尝试通过解决 Project Euler 上的一些任务来练习 Haskell。在问题 3 中，我们必须找到数字 的最大质因数600851475143，几年前我在 Java 中就已经做到了。

我想出了以下几点：

primes :: [Int]
primes = sieve [2..]
    where sieve (p:xs) = p : sieve (filter (\x -> x `rem` p /= 0) xs)

biggestPrimeFactor :: Int -> Int
biggestPrimeFactor 1 = 0
biggestPrimeFactor x =
    if x `elem` takeWhile (< x + 1) primes
        then x
        else last (filter (\y -> x `rem` y == 0) (takeWhile (< x `div` 2) primes))

这对于较小的数字非常有效，但效率非常低，因此对于我给出的数字来说效果不佳。这似乎很明显，因为程序会迭代所有小于除以 2 的数的素数（如果它本身不是素数），但我不确定该怎么办。理想情况下，我能够进一步限制可能的检查，但我不知道如何实现这一点。

请注意，我并不是在寻找“最佳解决方案”，而是在寻找一种对于更大的数字至少具有适度效率的解决方案，并且易于理解和实现，因为我仍然是 Haskell 的初学者。

Answer 1

这里有两个缓慢的主要原因。更容易解决的是 中的边界条件biggestPrimeFactor。检查到p > x `div` 2渐进地比检查到 更糟糕p^2 > x。但当一个数字有很多因素时，即使这样也不是最理想的。最大的因子可能远小于sqrt x。如果您在发现因素时不断减少目标数量，则可以考虑这一点并大大加快随机输入的处理速度。

这是一个例子，包括 Daniel Wagner 的评论注释：

-- Naive trial division against a list of primes. Doesn't do anything
-- intelligent when asked to factor a number less than 2.
factorsNaive :: [Integer] -> Integer -> [Integer]
factorsNaive primes@(p : ps) x
    | p * p > x = [x]
    | otherwise = case x `quotRem` p of
        (q, 0) -> p : factorsNaive primes q
        _ -> factorsNaive ps x

一些注意事项：

1. 我决定传入primes列表。这与下一节相关，但它也允许我在没有帮助者的情况下编写此内容。
2. 我专门使用 toInteger而不是，Int因为我想向它扔大量数字而不关心它maxBound :: Int是什么。这比较慢，但我决定首先默认正确性。
3. 我删除了对输入列表的遍历。一次性完成这件事效率更高一些，但大多数情况下它更干净。
4. 严格来说，即使输入列表包含非素数，只要列表从 2 开始，单调非递减，并最终包含每个素数，这也是正确的。
5. 请注意，当它递归时，它要么丢弃一个素数，要么产生一个素数。它永远不会同时做这两件事。这是确保不会错过重复因素的简单方法。

我命名它factorsNaive只是为了表明它并没有用数论做任何聪明的事情。有很多事情可以做，比这复杂得多，但这是一个很好的停止点，可以理解相对较小的数字的因式分解......

或者至少只要你有一个方便的质数列表，就可以进行因式分解。事实证明，这是代码速度变慢的第二个主要原因。随着素数列表变长，生成速度会很慢。

您的定义primes本质上是在输入列表上堆叠一堆过滤器。产生的每个素数都必须对之前的每个素数进行过滤测试。这听起来可能很熟悉 - 生成第一个素数至少需要 O(n^2) 的工作n。（实际上更多，因为随着数字变大，除法的成本更高，但现在让我们忽略这一点。）这是一个已知的（对于数学家来说，我必须查一下才能确定）结果，质数的数量小于或等于随着变大而n接近。随着变大，它接近线性，因此生成的列表随着变大而接近 O(n^2) 。n/ln nnnprimesnn

（是的，这个论证是一团糟。它的正式版本在梅丽莎·奥尼尔的论文“埃拉托斯特尼的真正筛子”中提出。请参阅它以获取对结果的更严格的论证。）

可以写出更有效的定义，primes它具有更好的常数因子和更好的渐近性。由于这就是上面括号中提到的论文的全部要点，因此我不会过多讨论细节。我只想指出第一个可能的优化：

-- trial division. let's work in Integer for predictable correctness
-- on positive numbers
trialPrimes :: [Integer]
trialPrimes = 2 : sieve [3, 5 ..]
  where
    sieve (p : ps) = p : sieve (filter (\x -> x `rem` p /= 0) ps)

这比你想象的要少。它不会使速度加倍，因为性能改进最终会被前面提到的过滤器堆栈所抵消。此版本仅从该堆栈中删除一个过滤器，但至少它是在初始版本中拒绝最多输入的过滤器。

在 ghci 中（没有编译或优化，这些确实可以产生影响），这足够快，可以在几秒钟内分解两个五位数素数的乘积。

ghci> :set +s
ghci> factorsNaive trialPrimes $ 84761 * 60821
[60821,84761]
(5.98 secs, 4,103,321,840 bytes)

具有几个小因子的数字处理速度要快得多。另请注意，由于素数列表是顶级绑定，因此计算结果会被缓存。再次运行计算现在已预先计算出素数列表。

ghci> factorsNaive trialPrimes $ 84761 * 60821
[60821,84761]
(0.01 secs, 6,934,688 bytes)

这也表明运行时间绝对由生成素数列表决定。当素数列表已经在内存中时，在这种规模下，朴素因式分解几乎是即时的。

但您不应该真正信任解释代码的性能。

main :: IO ()
main = print (factorsNaive trialPrimes $ 84761 * 60821)

给出

carl@DESKTOP:~/hask/2023$ ghc -O2 -rtsopts factor.hs
[1 of 2] Compiling Main             ( factor.hs, factor.o )
[2 of 2] Linking factor
carl@DESKTOP:~/hask/2023$ ./factor +RTS -s
[60821,84761]
   1,884,787,896 bytes allocated in the heap
      32,303,080 bytes copied during GC
          89,072 bytes maximum residency (2 sample(s))
          29,400 bytes maximum slop
               7 MiB total memory in use (0 MB lost due to fragmentation)

                                     Tot time (elapsed)  Avg pause  Max pause
  Gen  0       326 colls,     0 par    0.021s   0.021s     0.0001s    0.0002s
  Gen  1         2 colls,     0 par    0.000s   0.000s     0.0002s    0.0004s

  INIT    time    0.000s  (  0.000s elapsed)
  MUT     time    0.523s  (  0.522s elapsed)
  GC      time    0.021s  (  0.022s elapsed)
  EXIT    time    0.000s  (  0.007s elapsed)
  Total   time    0.545s  (  0.550s elapsed)

  %GC     time       0.0%  (0.0% elapsed)

  Alloc rate    3,603,678,988 bytes per MUT second

  Productivity  96.0% of total user, 94.8% of total elapsed

这将运行时间从六秒缩短到半秒。（是的，+RTS -s这相当冗长，但它快速且简单。）我认为这是停止初学者级代码的合理位置。

如果您想研究更有效的素数生成，hackage 上的 primes 包包含 O'Neill 论文中算法的实现以及与此处等效的朴素因式分解的实现。