凤凰木的笔记 - 欧拉计划第 50 题解及其优化

题目

Consecutive prime sum

The prime 41, can be written as the sum of six consecutive primes:

41 = 2 + 3 + 5 + 7 + 11 + 13

This is the longest sum of consecutive primes that adds to a prime below one-hundred.

The longest sum of consecutive primes below one-thousand that adds to a prime, contains 21 terms, and is equal to 953.

Which prime, below one-million, can be written as the sum of the most consecutive primes?

简而言之，就是要求出小于 1 百万的素数中，能表示成最长的连续素数之和的素数。

这个问题非常有意思。

初步的想法

先从最自然的想法开始。首先考虑求出 1 百万(以下简称 1M) 内的所有素数，然后对每一个素数，求出由此素数开始的最长连续素数序列，使得连续素数序列的和为 1M 内的素数，最后再根据长度对这些素数序列排序，取出最长的一个出来了事。

多么简洁明了，直奔主题的思路啊…虽然隐隐觉得事实上不会这么顺利，但实在无法抵抗这条思路的诱惑，嗯，来试下吧。

来列个 TODO List:

求素数，
对每个素数 p, 求出满足条件的最长素数序列,
对素数序列排序，取最长一个，求其和，即为本题的解。

为了研究和调试，我们可以把这个问题稍微一般化一点，把 1M 一般化为任意自然数 n (n>1).

求素数的算法我刚写过(看这里)，可以直接抄过来用:

primes :: Integer -> [Integer]
primes 0 = []
primes 1 = []
primes 2 = [2]
primes n = qs ++ [x | x <- [sqrtn..n], and [rem x y /= 0 | y <- qs]] where
    qs =  primes sqrtn
    sqrtn = floor $ sqrt $ fromInteger n + 1

乐观地猜测，这个算法的性能是能满足当下需求的。

接下来，第 2 步，对每个素数 p, 求出满足条件的最长素数序列。

max' :: Integer -> Integer -> [Integer] -> (Integer, Integer, Integer)
max' p n primesN = f 0 0 0 0 (dropWhile (<p) primesN) where
    f m m' s s' [] = (p, m', s')
    f m m' s s' (x:xs)
        | s > n = (p, m', s')
        | elem (s+x) primesN = f (m+1) (m+1) (s+x) (s+x) xs
        | otherwise = f (m+1) m' (s+x) s' xs

这里采用了尾递归的实现。

max' 接受三个参数, 其中，=p= 是目标素数，=n= 是上界, primesN 是不超过 n 的递增素数列表。

max' 返回一个三元组, 三元组内依次为 n, 最大连续素数列表长度, 素数列表之和(一定也是一个素数).

子函数 f 接受 5 个参数：=m= 是当前连续素数列表长度，=m'= 是其“有效值, ”=s= 是当前连续素数列表的和， s' 是其“有效值, 最后一个参数是待遍历的以 =p=开头的素数列表。

f 的返回类型与 max' 一致。

当以 p 开头的素数列表被遍历完了，或者累积和已经超过上界 n 时，=f= 返回累积值；否则，=f= 一边遍历一边检查累积起来的和是不是素数，同时存进累积器里。

接下来实现第 3 步：

p50 n = foldl c (0, 0, 0) $ map (\x -> max' x n ps) $ ps where
    c u@(_, y1, _) v@(_, y2, _) = if y1 < y2 then v else u
    ps = primes n

p50 比较简单，用一个比较函数 c 作左折叠而已。

最后写个 main 调一下:

import System.Environment (getArgs)

main = do
  (x0:_) <- getArgs
  let n = read x0 :: Integer
  putStrLn . show $ p50 n

这样就搞掂了！

编译好来跑一下:

weiwen@w ~/euler/50 % ghc p50-0.hs
[1 of 1] Compiling Main             ( p50-0.hs, p50-0.o )
Linking p50-0 ...
weiwen@w ~/euler/50 % time ./p50-0 10
(2,2,5)
./p50-0 10  0.00s user 0.00s system 81% cpu 0.005 total
weiwen@w ~/euler/50 % time ./p50-0 100
(2,6,41)
./p50-0 100  0.00s user 0.00s system 90% cpu 0.007 total
weiwen@w ~/euler/50 % time ./p50-0 1000
(7,21,953)
./p50-0 1000  0.06s user 0.00s system 98% cpu 0.068 total
weiwen@w ~/euler/50 % time ./p50-0 10000
(3,65,9521)
./p50-0 10000  5.92s user 0.03s system 99% cpu 5.954 total
weiwen@w ~/euler/50 % time ./p50-0 100000
(3,183,92951)
./p50-0 100000  767.78s user 4.78s system 99% cpu 12:52.56 total

果然不妥。我甚至没有勇气跑 n=1M…

不过稍有安慰的是, 可以看到 n=10, n=100, n=1000 时都是正确的。

第一步改进

考察输出的结果，发现当 $ n=10¹ $ 到 $ n=10⁵ $ 时，p50 n 的结果都以小素数开头。

由此猜测，是不是对大素数的枚举都是不必要的呢？

现在还无法证实这一臆想，但仔细想想会发现，至少在对素数 p 作 max' 运算之前，应该检查下有无必要。

以 $ n=1000 $ 为例。遍历至 $ p=7 $ 时，已有结果 $ (7, 21, 953) $, 这表明我们至少能找到一个长达 21 位的连续素数列满足要求。对于大于 7 的素数任意素数 p'，如果其值超过 $ [{1000 }] = 143 $, 那么以 p' 开头的素数序列最长显然不会超过 21, 因此 p' 不必考虑了。

由此，我们来改进 p50, 写成带有尾递归的形式:

p50' :: Integer -> (Integer, Integer, Integer)
p50' n = f (0,0,0) ps where
    ps = primes n
    f r [] = r
    f u@(x1, y1, z1) (x:xs) | y1 > 0 && x > quot n y1 = u
                          | y1 < y2 = f v xs
                          | otherwise = f u xs
                          where v@(x2, y2, z2) = (max' x n ps)

编译测试一下：

weiwen@w ~/euler/50 % ghc p50-1.hs
[1 of 1] Compiling Main             ( p50-1.hs, p50-1.o )
Linking p50-1 ...
weiwen@w ~/euler/50 % time ./p50-1 10
(2,2,5)
./p50-1 10  0.00s user 0.00s system 75% cpu 0.005 total
weiwen@w ~/euler/50 % time ./p50-1 100
(2,6,41)
./p50-1 100  0.00s user 0.00s system 90% cpu 0.006 total
weiwen@w ~/euler/50 % time ./p50-1 1000
(7,21,953)
./p50-1 1000  0.00s user 0.00s system 84% cpu 0.008 total
weiwen@w ~/euler/50 % time ./p50-1 10000
(3,65,9521)
./p50-1 10000  0.05s user 0.00s system 97% cpu 0.052 total
weiwen@w ~/euler/50 % time ./p50-1 100000
(3,183,92951)
./p50-1 100000  2.65s user 0.01s system 99% cpu 2.652 total
weiwen@w ~/euler/50 % time ./p50-1 1000000
(7,543,997651)
./p50-1 1000000  192.96s user 0.02s system 99% cpu 3:12.99 total
weiwen@w ~/euler/50 %

看，用了 3 分 13 秒，我们得到了本题的解：$ p50 1000000 = (7,543,997651) $!

进一步优化

回过头看 max' 的实现，和 s 的值是奇偶相错的，为偶数时，可以不必做素数性判断。因为基于 elem 函数的素数判断代价比较大，我们可以期望这个改进可以提升一点性能：

max'' :: Integer -> Integer -> [Integer] -> (Integer, Integer, Integer)
max'' p n primesN = f 0 0 0 0 (dropWhile (<p) primesN) where
    f m m' s s' [] = (p, m', s')
    f m m' s s' (x:xs)
    | s > n = (p, m', s')
    | even (s+x) = f (m+1) m' (s+x) s' xs -- if s+x is an even number, then ignore it
    | elem (s+x) primesN = f (m+1) (m+1) (s+x) (s+x) xs
    | otherwise = f (m+1) m' (s+x) s' xs

再编译跑一下:

weiwen@w ~/euler/50 % for i in {10,100,1000,10000,1000000}; do time ./p50-2 ${i}; done
(2,2,5)
./p50-2 ${i}  0.00s user 0.00s system 82% cpu 0.005 total
(2,6,41)
./p50-2 ${i}  0.00s user 0.00s system 72% cpu 0.004 total
(7,21,953)
./p50-2 ${i}  0.00s user 0.00s system 72% cpu 0.004 total
(3,65,9521)
./p50-2 ${i}  0.02s user 0.00s system 94% cpu 0.021 total
(7,543,997651)
./p50-2 ${i}  83.16s user 0.01s system 99% cpu 1:23.17 total

可见运行时间缩短 50% 以上，这个手段是显著有效的。

继续优化

重新考虑 max'' 的实现。我们是从最小的 p 开始，从左向右遍历素数列表的。这样做了很多无用功。猜测，如果反过来，从右向左遍历，应该会更加高效。

现在来尝试这个思路。

首先作一个素数列表的累积和列表。以 p=2 为例考察。受到 Haskell 官方 wiki 中求斐波那契数列(Fibonacci sequence)(原文)的启发，我们可以这样构建我们的素数和列表 ss：

   ps  = 2 3  5  7 11 13 17 19 ...
+  ss' =   2  5 10 17 28 41 58 ... 
=  ss  = 2 5 10 17 28 41 58 77 ...

Haskell 的懒惰求值(lazy evaluation)特性, 使得我们可以让一个有限或无限长的列表递归地调用其自身的部分。注意到:

ss' = tail ss

于是

ss = p : x : (zipWith (+) (drop 2 ps) (tail ss)) where x = (head $ tail ps) + p

为了对 ss 的项编号，我们使用 zip [1..] ss 这种伎俩，生成带编号的二元组。

接下来，滤掉超出上界 n 的项，再筛出值为奇数数的项，因为只有奇数才有可能是素数(2 已排除)：

ss' = filter (\(x, y) -> odd y) $ takeWhile (\(x, y) -> y<=n) $ zip [1..] ss

按照刚才拟定的计划，从尾部开始筛出和为素数的项。反序以后再筛，只要筛出一个，就可以把整个列表丢掉不管了。这里再次利用了 Haskell 的惰性求值特性:

ss'' = take 1 $ filter (\(x, y) -> elem y primesN) $ reverse ss'
(long, s) = head ss''

现在我们写出新的 max'' 实现了:

max''' :: Integer -> Integer -> [Integer] -> (Integer, Integer, Integer)
max''' p n primesN = (p, fromIntegral long, fromIntegral s) where
    ps = (dropWhile (<p) primesN)
    ss = p:x:(zipWith (+) (drop 2 ps) (tail ss)) where x = (head $ tail ps) + p
    ss' = filter (\(x, y) -> odd y) $ takeWhile (\(x, y) -> y<=n) $ zip [1..] ss
    ss'' = take 1 $ filter (\(x, y) -> elem y primesN) $ reverse ss'
    (long, s) = head ss''

编译跑一下:

weiwen@w ~/euler/50 % ghc p50-3.hs
[1 of 1] Compiling Main             ( p50-3.hs, p50-3.o )
Linking p50-3 ...
weiwen@w ~/euler/50 % for i in {10,100,1000,10000,1000000}; do time ./p50-3 ${i}; done
(2,2,5)
./p50-3 ${i}  0.00s user 0.00s system 76% cpu 0.004 total
(2,6,41)
./p50-3 ${i}  0.00s user 0.00s system 86% cpu 0.003 total
(7,21,953)
./p50-3 ${i}  0.00s user 0.00s system 77% cpu 0.004 total
(3,65,9521)
./p50-3 ${i}  0.01s user 0.00s system 97% cpu 0.012 total
(7,543,997651)
./p50-3 ${i}  4.88s user 0.02s system 99% cpu 4.900 total

非常显著。n=10^6 时的执行时间从 1 分 20 秒缩短至 5 秒。

连续跑几次，时间大概在 4.8 到 5.5 之间。

到这里，虽然还有几个地方可以进一步优化，但可以暂时先告一段落了。

最终的实现

一共 31 行:

import System.Environment (getArgs)
import Data.List

primes :: Integer -> [Integer]
primes 0 = []
primes 1 = []
primes 2 = [2]
primes n = qs ++ [x | x <- [sqrtn..n], and [mod x y /= 0 | y <- qs]] where
    qs =  primes sqrtn
    sqrtn = floor $ sqrt $ fromInteger n + 1

max' :: Integer -> Integer -> [Integer] -> (Integer, Integer, Integer)
max' p n primesN = (p, fromIntegral long, fromIntegral s) where
    ps = (dropWhile (<p) primesN)
    ss = p : x : (zipWith (+) (drop 2 ps) (tail ss)) where x = (head $ tail ps) + p
    ss1 = filter (\(x, y) -> odd y) $ takeWhile (\(x, y) -> y<=n) $ zip [1..] ss
    ss2 = take 1 $ filter (\(x, y) -> elem y primesN) $ reverse ss1
    (long, s) = head ss2

p50 :: Integer -> (Integer, Integer, Integer)
p50 n = f (0, 0, 0) ps where
    ps = primes n
    f r [] = r
    f u@(x1, y1, z1) (x:xs) | y1 > 0 && x > quot n y1 = u
                            | y1 < y2 = f v xs
                            | otherwise = f u xs
                            where v@(x2, y2, z2) = (max' x n ps)

main = do
  (x0:_) <- getArgs
  putStrLn . show $ p50 $ read x0

继续探索

继续增大 n 的值，求 p50:

eiwen@w ~/euler/50 % time ./p50 10000000
(5,1587,9951191)
./p50 10000000  129.84s user 0.22s system 99% cpu 2:10.07 total
weiwen@w ~/euler/50 % time ./p50 100000000
(7,4685,99819619)
./p50 100000000  3372.38s user 4.53s system 99% cpu 56:17.07 total

把 n 从 10^1 到 10^8 一起算出，我们得到:

weiwen@w ~/euler/50 % for i in {1..8}; do echo `./p50-3 $(( 10 ** i ))`"\t\t\\t\t${i}\t\t\t\t$(( 10 ** i ))" ; done
(2,2,5)               1    10        
(2,6,41)              2    100
(7,21,953)            3    1000
(3,65,9521)           4    10000
(3,183,92951)         5    100000
(7,543,997651)        6    1000000
(5,1587,9951191)      7    10000000
(7,4685,99819619)     8    100000000

如果我们用 P(n) 来表示不超过 n 的素数中，能表示成的最大连续素数和的长度，那么我们有:

\begin{align} P(10^1) & = 2 \\\\ P(10^2) & = 6 \\\\ P(10^3) & = 21 \\\\ p(10^4) & = 65 \\\\ p(10^5) & = 183 \\\\ p(10^6) & = 543 \\\\ p(10^7) & = 1587 \\\\ p(10^8) & = 4685 \\\\ \end{align}