如何用Haskell高效查找1到1000000000000范围内的所有完全数

核心原理优化（降维打击级优化） #

根据欧几里得-欧拉定理，所有偶完全数都可以表示为 2^(p-1) * (2^p - 1) 的形式，其中 2^p - 1 必须是梅森素数（且p本身为素数）。目前已证明1e12范围内不存在奇完全数，因此不需要遍历所有整数，只需要枚举符合条件的素数p生成候选即可，复杂度直接从O(n√n)降到O(1)（因为符合条件的p数量极少）。

注：你提到的第4个完全数有误，8192不是完全数，正确值为8128。

原代码问题 #

1. 真约数和计算bug：当输入值为完全平方数时，平方根因子会被重复相加，比如a=4时，你的函数计算得到的约数和是1+2+2=5，而实际真约数和为1+2=3，会导致判断错误。
2. 遍历逻辑效率极低：逐个数检查的逻辑在n超过1e6后就会出现明显性能问题，到1e12完全不可能在合理时间内跑完。

优化实现代码 #

-- 朴素素数判断，因为p上限极低，不需要复杂筛法
isPrime :: Integral a => a -> Bool
isPrime n
  | n <= 1 = False
  | n == 2 = True
  | even n = False
  | otherwise = all (\i -> n `mod` i /= 0) [3,5..floor (sqrt (fromIntegral n))]

-- 生成指定范围内的所有完全数
perfectNumbersUpTo :: Integer -> [Integer]
perfectNumbersUpTo limit =
  [ perfect | p <- filter isPrime [2..maxP]
            , let mersenne = 2^p - 1
            , isPrime mersenne
            , let perfect = 2^(p-1) * mersenne
            , perfect <= limit ]
  where
    -- 计算p的最大可能值，满足2^(p-1)*(2^p-1) <= limit
    maxP = floor $ logBase 2 (sqrt (fromIntegral (limit * 2 + 1)) + 0.5)

-- 调用示例：查找1e12以内的所有完全数
main :: IO ()
main = print $ perfectNumbersUpTo 1000000000000

运行结果 #

上述代码输出为 [6,28,496,8128,33550336,8589869056,137438691328]，即为1e12以内的所有完全数，运行时间小于1毫秒。

可选小规模优化（如果坚持用原遍历逻辑） #

如果一定要保留逐个数检查的逻辑，可以做以下优化：

• 修复sumDivisors的平方数问题：

  sumDivisors a =
    foldr (\n -> let (q,r) = a `quotRem` n in
        if r==0 then (+ if n == q then n else n + q) else id) 1
        [2..(floor . sqrt . fromIntegral $ a)]
  

• 只遍历偶数：因为所有已知完全数都是偶数，且小于1e12的奇数不可能是完全数，直接过滤掉所有奇数，速度提升一倍。
• 提前剪枝：sumDivisors计算过程中如果和已经超过n，可以直接提前终止计算，不用遍历到平方根。

内容的提问来源于stack exchange，提问作者VanGendt