在 Haskell 中寻找完美数

Question

因此，使用 Haskell 找到前 4 个完美数字（6、28、496、8192），并且有大量的实现。但目前我正在尝试设计一个有效的解决方案来找到从 1 -> 1,000,000,000,000 的所有完美数字。

我当前的代码如下所示：

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

isPerfect :: Integral a => a -> Bool
isPerfect n = n == sumDivisors n

upUntil :: Integral a => a -> [a]
upUntil n = filter isPerfect [1..n]

它可以相当快地找到 10000 个，但我需要能够在有效的时间空间中进一步查找。有什么建议吗？

Answer 1

您需要更好的算法，因此请在互联网上查找。与计算机程序员相比，数学家在如何有效地找到完美数字方面是更好的专家。 https://web.stanford.edu/class/cs106b/assignments/1-cpp/perfect是我发现的一种算法的描述，该算法似乎至少能够找到比穷举搜索算法更完美的数字。

Answer 2

在 Haskell 中，通过结合 Euclid-Euler 定理和 Miller-Rabin 素性测试，您可以快速找到非常大的完美数。如果没有米勒·拉宾，你可能会在运行时在第 8 个和第 12 个完美数之间碰壁。

这个要点解释了为什么欧几里得-欧拉和米勒-拉宾可以加快速度； Euclid-Euler 可让您搜索梅森素数（并将它们乘以 2 的幂）来找到完美数，而不是计算每个数字的等分总和。 请注意，如果存在奇数完美数，该算法将错过它们。

米勒·拉宾 (Miller Rabin) 走得更远，使我们能够以概率方式（高置信度）检查素性。它不是检查直到 n 的平方根的每个数字的模，而是依赖于 费马小定理 和“见证”数字，从而大大降低了复杂性。

这是我几年前编写的算法，使用来自 wiki.haskell.org 的 Miller-Rabin 测试以及我自己基于欧几里得-欧拉定理的代码：

{-| 
  Euclid-Based Solution with Miller-Rabin Primality Testing
  --
  Below, till marked end: miller-rabin funcs, from wiki.haskell.org
 -}
-- (eq. to) find2km (2^k * n) = (k,n)
find2km :: Integral a => a -> (a,a)
find2km n = f 0 n
    where 
        f k m
            | r == 1 = (k,m)
            | otherwise = f (k+1) q
            where (q,r) = quotRem m 2        

-- n is the number to test; a is the (presumably randomly chosen) witness
millerRabinPrimality :: Integer -> Integer -> Bool
millerRabinPrimality n a
    | a <= 1 || a >= n-1 = 
        error $ "millerRabinPrimality: a out of range (" 
              ++ show a ++ " for "++ show n ++ ")" 
    | n < 2 = False
    | even n = False
    | b0 == 1 || b0 == n' = True
    | otherwise = iter (tail b)
    where
        n' = n-1
        (k,m) = find2km n'
        b0 = powMod n a m
        b = take (fromIntegral k) $ iterate (squareMod n) b0
        iter [] = False
        iter (x:xs)
            | x == 1 = False
            | x == n' = True
            | otherwise = iter xs

-- (eq. to) pow' (*) (^2) n k = n^k
pow' :: (Num a, Integral b) => (a->a->a) -> (a->a) -> a -> b -> a
pow' _ _ _ 0 = 1
pow' mul sq x' n' = f x' n' 1
    where 
        f x n y
            | n == 1 = x `mul` y
            | r == 0 = f x2 q y
            | otherwise = f x2 q (x `mul` y)
            where
                (q,r) = quotRem n 2
                x2 = sq x
 
mulMod :: Integral a => a -> a -> a -> a
mulMod a b c = (b * c) `mod` a
squareMod :: Integral a => a -> a -> a
squareMod a b = (b * b) `rem` a

-- (eq. to) powMod m n k = n^k `mod` m
powMod :: Integral a => a -> a -> a -> a
powMod m = pow' (mulMod m) (squareMod m)

{-| End Miller-Rabin funcs
    Start Perfect-Finding funcs 
 -}
genpow2 :: Integer -> [(Integer, Integer)]
genpow2 n = [(x, y) | y <- [2..n], x <- [2^y]]

genprimes :: Integer -> [(Integer, Integer)]
genprimes n = filter 
    (\x -> 
        if fst(x)-1 < 4 
        then True 
        else all (== True) (map (millerRabinPrimality (fst(x)-1)) (filter (<= snd(x)) [2,3,5,7,11,13,17,31,73]))
    ) (genpow2 n)

perfects :: Integer -> [Integer]
perfects n = [(fst(x)-1)*(2^(snd(x)-1)) | x <- genprimes n]

main = print (perfects 256)

此代码足以在现代机器上在几分钟内找到非常大的完美数（1000 个数字）（如果允许它使用所有 ram 和 CPU，可能会更快）。使用

main = print (perfects 2281)

运行时，大约需要 14 分钟才能找到第 17 个完全数（1373 位）。