尝试优化此算法/代码以计算 a 和 b 之间互质整数的最大子集

Question

我正在用 Haskell 编写一个函数，它接收两个整数 a 和 b，并计算 [a,b] 的最大子集的长度，使得所有元素互质。现在，这样做的原因是因为我相信研究这些值可能会产生产生微不足道的较小界限的方法（可能一个大到足以暗示各种推测范围内的实际素数，例如连续平方）。所以很自然地，我试图对一些相当大的数字运行它。

不幸的是，下面的代码运行速度不够快，以至于我无法使用它。我认为缺陷是 Haskell 的惰性策略导致问题，但我既不知道语法也不知道我需要强制执行以防止操作累积的位置。

subsets [] = [[]]
subsets (x:xs) = subsets xs ++ map (x:) (subsets xs)

divides a b = (mod a b == 0)

coprime a b = ((gcd a b) == 1)

mutually_coprime []  = True
mutually_coprime (x:xs) | (coprime_list x xs) = mutually_coprime xs
                        | otherwise = False

coprime_list _ [] = True
coprime_list a (x:xs) | (coprime a x) = coprime_list a xs
                      | otherwise = False

coprime_subsets a b = coprime_subsets_helper (subsets [a..b])

coprime_subsets_helper [] = []
coprime_subsets_helper (x:xs) | (mutually_coprime x) = [x] ++ (coprime_subsets_helper xs)
                              | otherwise = coprime_subsets_helper xs

coprime_subset_length a b = max_element (map list_length (coprime_subsets a b))

list_length [] = 0
list_length (x:xs) = 1 + list_length xs

max_element a = max_element_helper a 0

max_element_helper [] a = a
max_element_helper (x:xs) a | (x > a) = max_element_helper xs x
                            | otherwise = max_element_helper xs a

只是为了弄清楚这是什么类型的输入，“coprime_subsets 100 120”对我来说似乎从未停止过。我实际上让它运行，起床，做了一些其他的事情，然后又回来了。 它仍在运行。我怀疑一个很大的瓶颈是立即计算所有子集。但是，我不想对生成的子集设置人为的下限。这可能会筛选掉所有互质集，让我一无所有。

到目前为止我所尝试的：

• 我用 gcd 替换了原来的互质函数。最初，这对所有整数使用模数和迭代检查。我认为 gcd 使用了类似 Euclid's Algorithm 的东西，理论上它应该运行得更快。

• 我一直在尝试想办法将子集的生成构建到互质集生成器中。到目前为止，我还没有想出任何东西。我也不确定这是否真的会有所帮助。

• 我一直在寻找 Haskell 的懒惰政策可能会伤害我的任何地方。没什么特别的，但我敢肯定。

我也知道这可能是我使用的环境 (winhugs) 的效率问题。我会说这就是问题所在；但是，当我询问如何在 Math Stack Exchange 上确定最大子集（对于一般 n 大小的范围）时，我收到的答案表明，从计算上讲，这是一件非常缓慢的事情。如果是这样，那很好。我只是希望能够跑完一些我感兴趣的范围，而不用花几乎永远的时间。

我知道这个网站一般不允许效率；但是，我已经尝试了很多，我不只是想在这里偷懒。我知道 Haskell 有很多奇怪的怪癖可以迫使它变得看似低效。我已经有一段时间没有在其中编写代码了，我怀疑我遇到了其中一个怪癖。

Answer 1

首先，我们必须找出慢的原因。

处理数字时，应使用提供所需精度和范围的最快表示。您的代码没有顶级类型签名，因此 GHC 将 coprime_subsets 的类型签名推断为

coprime_subsets :: Integral a => a -> a -> [[a]]

这允许 GHC 为您选择一个Integral，它会很乐意选择Integer，这比Int 慢得多。使用Integer，程序会花费大量时间来计算 gcd。强制 GHC 使用 Int 将运行时间从 6 秒缩短到 1 秒，因为 GHC 可以直接对机器整数进行整数数学运算。

注意：始终提供顶级类型签名也是一个好习惯。即使编译器没有从中受益，人类也经常这样做。

现在，进入问题的核心。跑步

main = print . length $ coprime_subsets (100 :: Int) 120
main :: IO ()

启用分析（stack build --profile 用于堆栈）并将+RTS -p -h 传递给可执行文件（-p 用于时间，-h 用于空间）给我们一个细分：

COST CENTRE            MODULE  SRC                            %time %alloc

subsets                Coprime src/Coprime.hs:(4,1)-(5,52)     52.5  100.0
coprime                Coprime src/Coprime.hs:11:1-26          25.5    0.0
coprime_list           Coprime src/Coprime.hs:(19,1)-(21,41)   18.5    0.0
coprime_subsets_helper Coprime src/Coprime.hs:(27,1)-(29,69)    1.8    0.0
mutually_coprime       Coprime src/Coprime.hs:(14,1)-(16,43)    1.7    0.0

当我们使用 Integer 时，绝大多数 (~78%) 时间都花在了 coprime 测试上。现在大部分都花在构建 powerset 上，所以让我们先看看那里。

接下来，我们要明白为什么会慢。

通常有三种策略可以加快速度：

1. 提高其渐近复杂度。
2. 提高其常数因子。
3. 减少调用次数。

哪些可能适用于subsets？ (1) 是一个显而易见的选择。构造幂集是O(2^n)，因此这里的任何渐近改进确实会非常有帮助。我们能找到吗？理论上，我们应该可以。正如丹尼尔所建议的，这个问题等价于最大团问题，它也是指数的。然而，极大团有一个基数较小的解，这意味着我们应该也能找到这个问题的渐近改进。

降低其渐近复杂度（以及我们调用它的次数）的关键在于，我们生成绝大多数子集只是为了稍后在检查它们的共质性时将它们丢弃。如果我们一开始可以避免生成错误的子集，我们将生成更少的子集并执行更少的互质性检查。如果修剪结果的数量与整个计算树的大小成正比，这将产生渐近改进。这种剪枝在函数式算法优化中很常见；您可以在Richard Bird's sudoku solver 中看到一个有启发性的示例。事实上，如果我们可以编写一个只生成非互质子集的函数，我们就已经解决了整个问题！

终于，我们准备好修复它了！

我们可以通过修改原始生成器subsets 来过滤掉非互质项：

coprimeSubsets [] = [[]]
coprimeSubsets (x:xs)
  = coprimeSubsets xs ++ map (x :) (coprimeSubsets (filter (coprime x) xs))

（如果我们认真考虑的话，我们可能可以在这里使用一些巧妙的折叠，但显式递归也很好。）

一旦我们这样做了，我们可以在大约 0.1 秒内找到[100..120] 的互质子集，提高了一个数量级。成本中心报告讲述了这个故事：

COST CENTRE    MODULE          SRC                            %time %alloc

MAIN           MAIN            <built-in>                      42.9    0.5
coprimeSubsets Coprime         src/Coprime.hs:(33,1)-(35,75)   28.6   67.4
CAF            GHC.IO.Encoding <entire-module>                 14.3    0.1
coprime        Coprime         src/Coprime.hs:13:1-26          14.3   31.1

现在我们实际上大部分时间都在做 IO 等。此外，我们coprime 测试的调用次数现在只有 3,848 次，提高了几个数量级。我们还可以在 3 秒内找到 [100..150] 的互质子集，相比之下……好吧，我没有等它完成，而是至少几分钟。

寻找加速的下一个地方可能是记忆 coprime 函数，因为这个问题涉及多次计算相同的参数。

Answer 2

我建议两个主要变化：

1. 只为每对数字计算一次gcd，而不是重复计算。
2. 仅扩展“可接受的”子集，而不是预先构建所有子集。

（这两个建议都与懒惰有关；我认为你的猜测是一个红鲱鱼。）下面我实施这两个建议；结果很简洁，我认为：

coprime_subsets' [] = [[]]
coprime_subsets' (x:xs)
    = coprime_subsets' xs
    ++ map (x:) (coprime_subsets' (filter ((1==) . gcd x) xs))

我们可以检查这是否在 ghci 中计算出相同的答案：

> coprime_subsets 10 20 == coprime_subsets' [10..20]
True
> coprime_subsets 100 120 == coprime_subsets' [100..120]
True

显然我的电脑比你的快得多：coprime_subsets 100 120 在这里即使在 ghci 中也只需不到 16 秒即可完成。当然，我的版本即使在 ghci 中也需要 0.02 秒，所以... =)

如果你只关心最大互质子集，你可以让它更快。这里的主要变化是在递归前面添加了filter：

maximal_coprime_subsets [] = [[]]
maximal_coprime_subsets (x:xs)
    = filter (any ((>1) . gcd x)) (maximal_coprime_subsets xs)
    ++ map (x:) (maximal_coprime_subsets (filter ((1==) . gcd x) xs))

尽管这个重复 gcds 的次数要多得多，但即使在 ghci 中它也可以在大约 0.01 秒内完成 [100..120]，因此它比之前的实现要快很多。

这里的时间差异似乎是输出 coprime_subsets' 只是花了更长的时间来打印，哈！

Answer 3

我不确定这是否是一种足够有效地解决问题的方法，但即使在 GHCI 中，它也会在 0.02 秒内返回所提供的两个整数值中的所有子集。

coprimes :: Int -> Int -> [[Int]]
coprimes x y | y <= x    = []
             | otherwise = (x : filter ((== 1) . gcd x) (tail [x..y])) : (coprimes . head . tail) [x..y] y

*Main> coprimes 100 120
[[100,101,103,107,109,111,113,117,119],[101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120],[102,103,107,109,113,115],[103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120],[104,105,107,109,111,113,115,119],[105,106,107,109,113,116,118],[106,107,109,111,113,115,117,119],[107,108,109,110,111,112,113,114,115,116,117,118,119,120],[108,109,113,115,119],[109,110,111,112,113,114,115,116,117,118,119,120],[110,111,113,117,119],[111,112,113,115,116,118,119],[112,113,115,117],[113,114,115,116,117,118,119,120],[114,115,119],[115,116,117,118,119],[116,117,119],[117,118,119],[118,119],[119,120]]
(0.02 secs, 811,576 bytes)