欧拉50计划:算法非常慢,无法理解原因



我正在使用Project Euler学习Haskell。我是Haskell的新手,在想出一个不需要花费大量时间的算法时遇到了很多麻烦。我估计这个项目需要140亿年才能找到解决方案。

问题:

小于一百万的哪个素数可以写成最数的和连续质数?

这是我的来源。我漏掉了isPrime。我把它贴出来是因为解决这个问题的效率太低了。我认为问题在于切片链和primeChains调用,但我不确定它是什么。我以前用c++解决了这个问题。但无论出于何种原因,在Haskell中高效的解决方案似乎超出了我的能力范围。

编辑:我已经包括了isPrime。

import System.Environment (getArgs)
import Data.List (nub,maximumBy)
import Data.Ord (comparing)
isPrime :: Integer -> Bool
isPrime 1 = False
isPrime 2 = True
isPrime x
        | any (== 0) (fmap (x `mod`) [2..x-1]) = False
        | otherwise = True
primeChain :: Integer -> [Integer]
primeChain x = [ n | n <- 1 : 2 : [3,5..x-1], isPrime n ]
slice :: [a] -> [Int] -> [a]
slice xs args = take (to - from + 1) (drop from xs)
    where from = head args
          to = last args
subsequencesOfSize :: Int -> [a] -> [[a]]
subsequencesOfSize n xs = let l = length xs
                          in if n>l then [] else subsequencesBySize xs !! (l-n)
    where
        subsequencesBySize [] = [[[]]]
        subsequencesBySize (x:xs) = let next = subsequencesBySize xs
                                    in zipWith (++) ([]:next) (map (map (x:)) next ++ [[]])
slicedChains :: Int -> [Integer] -> [[Integer]]
slicedChains len xs = nub [x | x <- fmap (xs `slice`) subseqs, length x > 1]
    where subseqs = [x | x <- (subsequencesOfSize 2 [1..len]), (last x) > (head x)]
primeSums :: Integer -> [[Integer]]
primeSums x = filter (ns -> sum ns == x) chain
    where xs = primeChain x
          len = length xs
          chain = slicedChains len xs
compLength :: [[a]] -> [a]
compLength xs = maximumBy (comparing length) xs
cleanSums :: [Integer] -> [[Integer]]
cleanSums xs = fmap (compLength) filtered
    where filtered = filter (not . null) (fmap primeSums xs)
main :: IO()
main = do
    args <- getArgs
    let arg = read (head args) :: Integer
    let xs = primeChain arg
    print $ maximumBy (comparing length) $ cleanSums xs

你的基本问题是你没有根据目前为止找到的最佳解决方案来修剪你的搜索空间。

我可以从您使用maximumBy查找最长序列的事实中看出这一点。

例如,如果在搜索过程中找到一个连续的4个素数序列,它们的和是素数<10^6,你不需要检查任何以大于250000的素数开始的序列。

要做这种修剪,你必须跟踪到目前为止找到的解决方案,并将候选序列的测试与它们的生成交叉进行,以便到目前为止找到的最佳解决方案可以提前停止搜索。

slicedChains中存在几个低效率。Haskell列表是通过链表实现的。这个视频很好地概述了链表以及它们与数组的区别:(link)

代码中的下列表达式将会影响w.r.t.效率:

* nub的运行时间是二次的

* length x > 1 - length的复杂度为O(n),其中n为列表的长度。更好的写法是:

lengthGreaterThan1 :: [a] -> Bool
lengthGreaterThan1 (_:_:_) = True
lengthGreaterThan1 _       = False

* subsequencesOfSize 2 [1..len]可以写得更简洁:

[ [a,b] | a <- [1..len], b <- [a+1..len] ]

,这也将确保

所以你需要的第一件事是一个好的数据结构:一旦你找到一个长度为n的序列,你就不关心更短长度的序列了,所以你的主要需求是:(1)跟踪总和,(2)跟踪集合中的素数,(3)删除最小元素,(4)添加一个新的最大元素。关键是摊销,其中支付的大成本很少,以至于您可以假装它是每个过程的小成本。数据结构如下所示:

data Queue x = Q [x] [x]
q_empty (Q [] []) = True
q_empty _ = False
q_headtails (Q (x:xs) rest) = (x, Q xs rest)
q_headtails (Q [] xs) = case reverse xs of y:ys -> (y, Q ys [])
                                           []   -> error "End of queue."
q_append el (Q beg end) = Q beg (el:end)

所以解构列表是可能的,但有时会触发一个O(n)的操作,但这没关系,因为当它这样做时,我们不必为另一个n步骤做这个操作,所以它平均为每一步一个操作。(你也可以用一个严格的列表来做。)

为了节省长度操作和列表项的求和,您可能也需要缓存这些操作:

type Length = Int
type Sum = Int
type Prime = Int
data PrimeSeq = PS Length Sum (Queue Prime)
headTails (PS len sum q) = (x, PS (len - 1) (sum - x) xs)
  where (x, xs) = q_headtails q
append x (PS len sum xs) = PS (len + 1) (sum + x) (q_append x xs)

算法如下:

  1. 缓存以
  2. 开头的PrimeSeq的副本
  3. 继续添加素数并测试素数,直到得到10^6。
  4. 如果发现新的素数序列较长,则替换缓存
  5. 当你遇到10^6时,返回到缓存,从队列的前面取出一个素数,然后根据需要重复。

你的素数代是二次的(isPrime 101测试rem 101 100 == 0,即使10101需要测试的最大数字——实际上7就足够了)。

然而,即使有了它,一个足够简单的基于列表的代码在2秒内找到答案(在Intel Core i7 2.5 GHz上,在GHCi中解释)。通过对代码进行修正以利用上面提到的优化(另外,仅通过素数进行测试),它需要0.1s

同样,f x | t = False | otherwise = Truef x = not t是一样的。

我们被PE网站要求不要给你任何提示。

但是一般来说,在Haskell中效率的关键,由于它的惰性,是在尽可能少的重复工作的情况下生成。举个例子,我们可以在一个进程中生成一堆片段,而不是从头开始单独计算列表的每个片段,

 slices :: Int -> [a] -> [[a]]
 slices n = map (take n) . iterate tail -- sequence of list's slices of length n each

另一个原则是,尝试解决更一般问题,您的问题就是一个实例。

写了这样一个函数后,我们可以通过为它的参数尝试不同的值来玩它,从小到大,以一种探索性的方式解决问题。已知21个连续素数。22呢? 27日?1127 ?…我已经说得够多了。

如果它开始花费太多时间,我们可以通过增长分析的经验阶数来评估整个解决方案所需的运行时间。

尽管使用未优化的isPrime代码可以很快找到解决方案,但探索过程可能会非常慢,但使用优化的代码则足够快:

 primes :: [Int]
 primes = 2 : filter isPrime [3,5..]
 isPrime n = and [rem n p > 0 | p <- takeWhile ((<= n).(^2)) primes]

最新更新