Project Euler的解题笔记

关于问题1的描述：

Multiples of 3 and 5

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23. Find the sum of all the multiples of 3 or 5 below 1000.

也就是求小于1000的3或5的倍数的数的和。鉴于这个题目挺简单，我就直接在Emacs里面用Elisp给解决了（我这个解题思路可能有点怪:-D）：

(funcall
 (lambda (max)
   (let ((m3 3) (m5 5) (sum 0))
     (while (< m3 max)
       (setq sum (+ sum
                    (if (= (% m3 5) 0) 0 m3)
                    (if (< m5 max) m5 0)))
       (setq m3 (+ m3 3))
       (setq m5 (+ m5 5)))
     sum))
 1000)

个人觉得这个解法比PE论坛上大部分从3开始遍历到1000的解法要好，因为这样迭代直接去掉了不是3或者5的倍数的验证过程（自从看了SICP之后，对迭代越来越感兴趣了:-p）。程序中之所以存在判断m3是不是5的倍数的那部分，是因为m3有可能是5的倍数，这样会导致其被加两遍，所以，碰到这种情况要将其除掉。

但这个并不是最优解。每个问题解决之后，PE会给出一个“标准”答案，这个问题的“标准”答案如下：

sum_multi(1000) = sum_multi_3(1000) + sum_multi_5(1000) - sum_multi_15(1000)
sum_multi_3(1000) = 3 + 6 + 9 + ... + 999
                  = 3 × (1 + 2 + ... + 333)
                  = 3 × (1 + 333) × 333 ÷ 2
sum_multi_5(1000) 和 sum_multi_15(1000) 类推

依据上述解法，写出Elisp的程序如下（这个解法的性能要好上XX倍 :-D）：

(defun sum(factor max)
  (let ((count (/ max factor)))
    (* factor
       (/ (* count
             (+ 1 count))
          2))))

(funcall
 (lambda (max)
   (- (+ (sum 3 max)
         (sum 5 max))
      (sum 15 max)))
 999)

关于问题2的描述：

Even Fibonacci numbers

Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, …

By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.

也就是求斐波那契数列所有项值不超过四百万并且为偶数的和。同样，直接在Emacs里面用Elisp解决（解题思路同样比较怪:-D）：

(funcall
 (lambda (max)
   (let ((3i 2) (3i+1 3) (3i+2 5) (sum 0))  ;; 这里的3i+1, 3i+2是变量，不是运算
     (while (<= 3i max)
       (setq sum (+ sum 3i))
       (setq 3i (+ 3i+1 3i+2))
       (setq 3i+1 (+ 3i 3i+2))
       (setq 3i+2 (+ 3i 3i+1)))
     sum))
 4000000)

因为注意到序号为3的倍数（我习惯以1,1,2来开始斐波那契数列，所以序号为3的倍数，如果以问题描述中为例，这个规律应该是2+3n）的项是偶数，所以只需要处理这些项就可以了。同问题1，采用迭代的方式，计算基本都是加法，连PE论坛中大多数普通循环解法中的求模运算都没用到。

PS：因为这个解法就是“标准”答案的思路，所以这里就不分析“标准”答案了。 :-p

关于问题3的描述：

The prime factors of 13195 are 5, 7, 13 and 29.

What is the largest prime factor of the number 600851475143 ?

也就是求数字600851475143的最大质数因子。

这个题目就有点讲究了，数字比较大，而分解大数也没什么好的办法，基本是需要一轮一轮的循环，而判断一个数是不是质数，基本也是需要一轮一轮循环去除，直到它的平方根为止。

但是，这里面还是有可以优化的地方，就是用来测试一个数是否为质数的算法，本来应该是O(n)的复杂度，我们可以将其优化至O(log n)，原理，就是著名的费马小定理¹：

费马小定理：如果n是一个素数，a是小于n的任意正整数，那么a的n次方与a模n同余。

根据费马小定理，可以引申出费马检查算法：选取一个小于n的a，如果n和a满足费马小定理，那么n是素数的可能性就很大，如果不满足，那么n一定不是素数。在满足定理的情况下，我们可以再选择其它的a来测试，这样，在选取一定次数的a之后，如果定理一直满足，我们就可以认为n是素数。

费马检查不是一个精确的检查，而只是表明一个数是素数的概率，但是，由于对于非素数n，大多数的a < n都不会满足定理，所以，如果有一个a能通过检查，n是素数的概率就大于50%，如果有两个a能通过检查，n是素数的概率就大于75%，因此，只需要少量的a，就可以让概率达到很高的值。

下面的代码是费马检查的Elisp实现：

(defun even? (n)
  (= (% n 2) 0))

(defun square (n)
  (* n n))

(defun expmod (base exp m)
  (cond ((= exp 0) 1)
        ((even? exp) (% (square (expmod base (/ exp 2) m)) m))
        (t (% (* base (expmod base (- exp 1) m)) m))))

(defun fermat-test (n)
  (defun try-it (a)
    (= (expmod a n n) a))
  (try-it (+ 1 (random (- n 1)))))

(defun fast-prime? (n times)
  (cond ((= times 0) t)
        ((fermat-test n) (fast-prime? n (- times 1)))
        (t nil)))

fast-prime? 函数接收两个参数，一个要用来测试的数n，一个是执行费马检查的次数。整个检查的耗时部分在 expmod 函数，但是它的复杂度也只有O(log n)。

基于上面的基础，我们就可以得出问题3的解法（Elisp实现）：

(funcall
 (lambda (n test-times)
   (let ((small-factor 1) (large-factor n)
         (found nil) result)
     (while (and (not found) (<= small-factor large-factor))
       (when (= (% n small-factor) 0)
         (when (fast-prime? small-factor test-times)
           (setq result small-factor))
         (setq large-factor (/ n small-factor))
         (when (fast-prime? large-factor test-times)
           (setq result large-factor)
           (setq found t)))
       (setq small-factor (+ 1 small-factor)))
     result))
 600851475143 2)

这个实现是比较丑陋的，递增地进行因数测试相当耗时间，幸好是两个因数同时测试，这样可以节省不少时间。另外，Lisp在这里的优势就体现出来了：对于600851475143这样的大数可以直接处理，而在C或者其它一些语言中，就不得不考虑溢出问题了。另外需要说明的是，上面的费马检查的次数设定为2，已经是足够了。

下面来看看“标准答案”，直接先贴上伪代码：

n = "the evil big number"

if n mod 2 = 0 then
    lastFactor = 2
    n = n div 2
    while n mod 2 = 0
        n = n div 2
else
    lastFactor = 1

factor = 3
maxFactor = sqrt(n)
while n > 1 and factor <= maxFactor
    if n mod factor = 0 then
        n = n div factor
        lastFactor = factor
        while n mod factor = 0
            n = n div factor
        maxFactor = sqrt(n)
   factor=factor+2
if n = 1 then
    output lastFactor
else
    output n

“标准答案”基于两个事实：1. 任何正整数都可以被分解为多个素数因子的的积（如果把1也当作素数的话）；2. 一个正整数，最多只能有一个大于其平方根的素数因子。

因此，“标准答案”就是一个不停解质因数的过程，小学的时候就学过，可是我看这个算法还是看来好久才明白原理，可怜了我那无下限的智商。。。

关于问题4的描述：

A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 × 99.

Find the largest palindrome made from the product of two 3-digit numbers.

即求两个3位数相乘所能得到的最大回文数。所谓回文数，即从最高位到最低位的数字顺序与最低位到最高位的数字顺序完全一样，例如12321。

这个算是比较简单的一个题目，依旧使用Elisp：

(defun check-palindromic (n)
  (let (numbers)
    (setq numbers (list (% n 10)))
    (setq n (/ n 10))
    (while (> n 0)
      (nconc numbers `(,(% n 10)))
      (setq n (/ n 10)))
    (equal numbers (reverse numbers))))

(funcall
 (lambda (min max)
   (let ((s max) (l max) (result 0) temp)
     (while (>= s min)
       (setq l s)
       (while (>= l min)
         (setq temp (* s l))
         (when (and (< result temp)
                    (check-palindromic temp))
           (setq result temp))
         (setq l (1- l)))
       (setq s (1- s)))
     result))
 100 999)

首先定义一个 check-palindromic 函数用于判断一个数是不是回文数，做法就是取出所有数字组成一个列表，然后检查反转后的列表和原列表是否相同。程序主体部分定义了一个lambda函数，以1为步长进行双重循环逐个判断即可。这里需要注意的就是类似10=2×5=5×2两个因数位置互换的情况，需要在双重循环中过滤掉。

至于“标准答案”中的最优解，采用分解因式的方法，确定满足要求的数必定可以分解为11和另外一个数的乘积（请自行证明），所以两个因子中有一个必定是11倍数，从而能让普通解法中的循环的递增/递减步长从1变为11，这大大减少了运行时间。

但不得不说的是，这个最优解是有前提的，即已经限定满足要求的数是六位数，所以不具有普遍性，因此我也就不贴最优解的解法了，因为普通解对这个题目已经是足够了，反而感觉最优解有点画蛇添足。

关于问题5的描述：

2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.

What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?

即求从1到20的所有数的最小公倍数。

难道要我告诉你们，我是用手算出这道题的么！！我写了段程序试了试，结果不对！！于是觉得1到20范围也不大，果断手算之，然后答案就对了！！

其实是我最初理解错了题意，写成了求所有数的质因子的并集的最小公倍数，后来发现，所有数的最小公倍数和所有数的质因子的并集的最小公倍数并不是一回事，比方说2和8，它们的最小公倍数是8，但质因子并集的最小公倍数却是2，就是因为在对8分解质因数的过程中，相同的质因子2被忽略了，但需要注意的是，相同的质因子也不能全部保存，例如2、4和8，最小公倍数是8，但如果保存所有分解的质因子，会发现乘积会是64，这显然是不正确的。因此，需要对这种分解出来只有一个质因子的数，即某一个质因子的次幂，保证有且只有一个质因子被乘到乘积中。

虽然题目已经用手解决了（PS：我真邪恶 :-p），但还是贴出解题程序，可以在解一些大数的时候使用：

(defun prime-factors (num)
  (let ((i 2) (factors '(1)))
    (while (>= num i)
      (when (= (% num i) 0)
        (setq num (/ num i))
        (add-to-list 'factors i t)
        (while (= (% num i) 0)
          (setq num (/ num i))))
      (setq i (+ i (if (= i 2) 1 2))))
    factors))

(funcall
 (lambda (max)
   (let ((i 2) (primes '(1)) (result 1) factors)
     (while (<= i max)
       (setq factors (prime-factors i))
       (dolist (f (cdr factors))
         (when (or (not (memq f primes))
                   (= (length factors) 2)) ;; 长度为2表明这是一个幂类型的数
           (setq result (* result f))
           (add-to-list 'primes f t)))
       (setq i (1+ i)))
     result))
 20)

然后看了一下“标准答案”，和我上面的原理是一样，也是首先列出所有质因子的并集，但它在处理每个质因子应该乘的次数的时候，思路稍有不同，以70为例，我在遇到2，4，8，16，32，64的时候，都会乘以2，而它则是一次确定质因子2应该的相乘次数，即 floor(log(70) / log(2)) ，因为程序是不可能一眼就看出来64是小于70的最大的2的次幂的，所以，通过对70求2为底的对数，再向下取整，就是2应该相乘的次数。其它质因子的相乘次数可以用相同的方法求出来，但对于平方已经大于70的质因子以及比它大的质因子，就没必要求了，因为它们需要相乘的次数肯定为1。

关于问题6的描述：

… Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

问题的介绍有点长，所以我省略了前面的例子，意思就是求前100个自然数的和平方与平方和的差。

这个问题其实不算是编程问题，稍作推导就知道：

${\displaystyle S(n)=(1+2+3+...+n{)}^2-<!--\; -\; -->(1^2+2^2+3^2+...+n^2)=[\frac{n(1+n)}{2}}{]}^2-<!--\; -\; -->{\displaystyle \frac{n(n+1)(2n+1)}{6}=...=\frac{n^4}{4}+\frac{n^3}{6}-<!--\; -\; -->\frac{n^2}{4}-<!--\; -\; -->\frac{n}{6}}$

所以，这个题目手算也是可以的…如果非要用程序求的话…

(funcall
 (lambda (n)
   (-
    (+ (/ (expt n 4) 4)
       (/ (expt n 3) 6))
    (+ (/ (expt n 2) 4)
       (/ n 6))))
 100)

“标准答案”就不介绍了，没有这个简单。:-p

关于问题7的描述：

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.

What is the 10 001st prime number?

即求第10001个素数。

还记得问题3里面的费马素数测试吗？我心里暗喜，刚好用前面的代码就可以了，于是赶紧写出以下代码：

(defun solve-problem7 (n)
  (let ((i 2)
        (prime 3))
    (while (< i n)
      (setq prime (+ prime 2))
      (when (fast-prime? prime 2) ;; 这里的fast-prime?的定义在问题3里面可以看到
        (setq i (+ i 1))))
    prime))

(solve-problem7 10001)

告诉你们，我在求第10001个素数之前，我还做过测试的！！我试着求过前面15个素数，结果都是正确的！！

于是，我信心满满地把上面求出来的结果——104327输入，结果错误。。。我再运行一遍，结果是104347，奇怪，再运行一遍，结果是104183。。。怎么结果还是随机的？一看代码，发现是费马测试的次数太少——每个用来测试的数只做了两次费马测试，虽然每次测试正确的概率很大，但总数太多，总会有漏网之鱼，所以每次的结果都不一样，于是把费马测试的次数改成10，这样，漏网之鱼的概率就很小了，试着运行了三次，结果都是104579。

于是，再次信心满满地把104579输入，结果还是错误。。。我再运行了三次，结果都是104579，那这肯定不是因为费马测试次数的原因了。这时，我才突然想起来影响结果的罪魁祸首——Carmichael数。

所谓Carmichael数，就是能完全通过费马测试（注意是完全通过，不是概率性通过）但不是素数的数。所以，就算是把费马测试的次数增加到一千次一万次，这些数也会安然无恙地通过费马测试。

所以，因为Carmichael数的存在影响了我们的计数，因而真正的结果肯定要比104579大。比104579小的Carmichael数有（这个数列的下一个是115921）：

561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657, 52633, 62745, 63973, 75361, 101101

可以看到，有17个Carmichael数骗了我们，于是，从这里找到素数表，从104579再向后数17个偏移量，结果就出来了，我就不再去写程序求了。:-p

“标准答案”的素数判定原理和普通的循环求模判定差不多，只是稍做了一些优化，这里就不再介绍了。

关于问题8的描述：

Find the greatest product of five consecutive digits in the 1000-digit number.

73167176531330624919225119674426574742355349194934 96983520312774506326239578318016984801869478851843 85861560789112949495459501737958331952853208805511 12540698747158523863050715693290963295227443043557 66896648950445244523161731856403098711121722383113 62229893423380308135336276614282806444486645238749 30358907296290491560440772390713810515859307960866 70172427121883998797908792274921901699720888093776 65727333001053367881220235421809751254540594752243 52584907711670556013604839586446706324415722155397 53697817977846174064955149290862569321978468622482 83972241375657056057490261407972968652414535100474 82166370484403199890008895243450658541227588666881 16427171479924442928230863465674813919123162824586 17866458359124566529476545682848912883142607690042 24219022671055626321111109370544217506941658960408 07198403850962455444362981230987879927244284909188 84580156166097919133875499200524063689912560717606 05886116467109405077541002256983155200055935729725 71636269561882670428252483600823257530420752963450

即从这个1000位数中，找出乘积最大的5个连续数字。

这个题目也是不需要代码的，既然要求5个连续数字乘积最大，那么当然是要这5个数字中9越多越好，剩下的尽量大，而且不可以出现0。所以，直接在这个数字中搜索9，然后以它为中点，加左边两位，右边两位，看这5个数字包含9的数量以及剩下的数字的大小。这样，两分钟之内，应该就可以轻松找到结果了。（PS：我就是这样找到的 :-p）

去这个题目的论坛看了一下，很多人都是像我这样猜出来的，看来英雄所见略同嘛。:-D

另外，这个题目没有提供“标准答案”，看来出题者也是觉得写代码循环来求没什么意思，还不如猜呢。:-D

关于问题9的描述：

A Pythagorean triplet is a set of three natural numbers, a < b < c, for which, a² + b² = c²

For example, 3² + 4² = 9 + 16 = 25 = 5².

There exists exactly one Pythagorean triplet for which a + b + c = 1000.
Find the product abc.

描述中所谓的a, b, c就是满足勾股定理的数，也称为毕达哥拉斯三元组。题目要求的就是满足和等于1000的这样的三元组的乘积。

这其实就是一个解方程的过程，三个未知数，两个限定方程，所以最终只能化简得到一个函数，下面就是我最终化简得到的函数（为了一般性，我把限定条件中的1000换成了常数s）：

$a={\displaystyle \frac{s}{2}}\times <!--\; \times \; -->{\displaystyle \frac{s-<!--\; -\; -->2b}{s-<!--\; -\; -->b}}$

另外还有一个隐含条件，就是所有的数字必须是整数，所以上面的等式右面的分子必须可以被分母整除，根据这个约束，有以下代码：

(funcall
 (lambda (sum)
   (let ((b 1) (found nil) (half (/ sum 2))
         a1 a2)
     (while (and (< b half) (not found))
       (setq a1 (* half (- sum (* 2 b))))
       (setq a2 (- sum b))
       (if (= (% a1 a2) 0)
           (setq found t)
         (setq b (1+ b))))
     (* (/ a1 a2) b (- sum (/ a1 a2) b))))
 1000)

“标准答案”咱就不看了，没看太懂。。

关于问题10的描述：

The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.

Find the sum of all the primes below two million.

描述挺简单，就是求小于两百万的所有素数的和。

又是素数。。关于素数的测试，自然不能用费马测试了，不可能把小于两百万的所有Carmichael数都列出来排除掉。。

不过，在SICP的练习1.28中，提到一种改进的素数测试方法，叫Miller-Rabin检查，可以避开Carmichael数，刚好我也在做SICP的习题，于是把这个题做了一下，然后有下面的代码（因为SICP的代码都是用Scheme写的，所以图方便，就不把相关的函数一一porting到Elisp了，直接用Scheme写了）：

(define (find-prime-sum max)
  (cond ((> max 2)
         (+ (if (miller-rabin-prime? max 5) max 0)
            (find-prime-sum (- max (if (even? max) 1 2)))))
        (else 2)))

(find-prime-sum 2000000)

上面用到的 miller-rabin-prime? 过程，在我的SICP习题解可以找到，这里就不再抄一遍了。

至于“标准答案”，我搞出一个能正确判断素数的函数就很不容易了，结果它宣称能在一秒内算出答案，太伤心了，不看鸟。。

附：上面的话当然是气话，我后来查了一下，“标准答案”用的解法叫“埃拉托斯特尼筛法”，在wikipedia上有详细的介绍。介绍中还详细解释了这种解法和黎曼猜想的联系，看来数学各个分支都是相通的。。

问题11的描述比较长，不好贴，贴出传送门：

http://projecteuler.net/problem=11

即求一个二维矩阵上，横线、竖线或对角线上连续四个数字最大的积。这个问题和问题8有些类似，不同的是，这次是二维的了。

同样，这个问题也不需要代码，只需要到矩阵中搜索9即可，并且要是十位数字是9的，然后就靠感觉去看横竖斜的数字的乘积大小了，试了五次，正确结果就出来了。

这个题目也没有“标准答案”，没什么意思的题目，也跟数学没什么关系，所以不需要“标准答案”。。