Project-Euler第69题

大学的时候挺喜欢解Project Euler上的题目的，尤其是它不在乎答题者使用哪一门编程语言，甚至还有很多参与者是使用pen&paper来解题的。去年开始重新开始做Project Euler上的题目，而第69题则是最近刚刚解决的一题。惭愧的是，因为不晓得欧拉函数的计算公式（甚至都没有想过欧拉函数有没有可以用来计算的公式），所以这一题我是用暴力计算的方法来解决的。尽管花了40分钟左右才找出了问题的答案，但欧拉函数的计算方法本身还是让我觉得挺有意思的，下面我就来讲讲我在计算欧拉函数方面做的一些尝试。

69题本身很容易读懂，就是要找到一个不大于一百万的正整数n，这个n与以它为参的欧拉函数值的比值能够达到最大。欧拉函数的介绍可以看维基百科，简而言之，就是在不大于n的正整数中与n互质的数的个数，具体的例子可以参见69题描述中给出的表格。

网上可以找到别人的解法，基本的思路是：按从小到大的顺序，对于一百万以内的每个素数，都计算出它们的倍数的欧拉函数值的一部分——即对于素数p计算出1-1/p并与这个位置上原来的值相乘。当遍历完一百万以内的所有素数后，也就计算出每一个位置上的欧拉函数值，再遍历一次就可以计算出比值最大的数字了。

但我今天要讲的是笨方法。

笨方法的关键就是乖巧地计算出每一个数字的欧拉函数值。而其中最笨的，当属挑出每一个不大于n的因数，计算它们与n的最大公约数，并根据这个最大公约数是否为1，来决定是否给某一个计数器加一。示例代码如下

(defun phi (n)
  (let ((count 0))
    (dotimes (i n count)
      (when (= (gcd (1+ i) n) 1)
        (incf count)))))

这个phi可以稍微改进一下。例如，如果一个数a与n不是互质的，那么a的倍数（小于n的那一些）也一定不会与n互质。因此，每当遇到这么一个因数a，就知道后续的2a、3a等等，都不再需要计算其与n的最大公约数了。改进后的算法代码如下

(defun phi (n)
  "通过将不互质的比特设置为1并计算为0的比特的个数来计算phi函数"
  (let ((bits (make-array n :element-type 'bit :initial-element 0))
        (count 0))
    (dotimes (i n)
      (cond ((= (bit bits i) 1)
             ;; 该比特已经为1，说明已经在比它小的倍数被处理时一并被标记了
             )
            ((= i (1- n))
             ;; 只处理比上界要小的数字
             )
            ((/= (gcd (1+ i) n) 1)
             ;; 除了当前这个不互质的数字之外，还需要将这个数字的倍数也一并处理
             (dotimes (j (floor n (1+ i)))
               (let* ((j (1+ j))
                      (m (* (1+ i) j)))
                 (setf (bit bits (1- m)) 1))))
            (t (incf count))))
    count))

为了节省内存空间，这里用了一个bitmap来标记小于n的每一个数字是否与n互质——1表示不互质，0表示互质。

其实并不需要遍历所有比n小的数字，只要遍历所有n的素因子即可。比如，将8分解为素因子，就是3个2，那么对于小于8的所有2的倍数（4和6）都是与8不互质的。基于这个方法，将所有的素因子的倍数所对应的位置为1，再数一下总共有多少个0比特即可。

对每个n都进行质因数分解效率不高，先生成一个一百万以内的素数表吧

(defun primep (n)
  (cond ((< n 2) nil)
    ((= 2 n) t)
    (t
     (let ((bnd (truncate (sqrt n))))
       (labels ((rec (test)
              (cond ((> test bnd) t)
                ((= 0 (rem n test)) nil)
                (t (rec (1+ test))))))
         (rec 2))))))

(defvar *primes-in-1000000* nil)
(defun generate-primes-in-1000000 ()
  (dotimes (i 1000000)
    (when (primep (1+ i))
      (push (1+ i) *primes-in-1000000*)))
  (setf *primes-in-1000000* (nreverse *primes-in-1000000*)))

然后对于一个给定的n，遍历所有小于它的素数并对相应的倍数所在的比特置一就可以了，示例代码如下

(defun phi3 (n)
  "直接用素数表来做筛法"
  (prog
      ((bits (make-array n :element-type 'bit :initial-element 0)))
     (dolist (num *primes-in-1000000*)
       (cond ((> num n)
              (go :return))
             ((zerop (mod n num))
              (labels ((aux (i)
                         (when (< i n)
                           (setf (bit bits i) 1)
                           (aux (+ i num)))))
                (aux (1- num))))))
     :return
     (return-from phi3
       (count-if (lambda (bit) (zerop bit)) bits))))

PS：写这个phi3的时候发现Common Lisp提供了一个prog宏，这个宏倒是真的挺好用。

改进了两轮，其实这仍然是笨方法。即便是用phi3，用来计算题目的答案也花了40多分钟。

全文完