Project-Euler第69题

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

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

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

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

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

1
2
3
4
5
(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,就知道后续的2a3a等等,都不再需要计算其与n的最大公约数了。改进后的算法代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
(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都进行质因数分解效率不高,先生成一个一百万以内的素数表吧

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
(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,遍历所有小于它的素数并对相应的倍数所在的比特置一就可以了,示例代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
(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多分钟。

全文完