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多分鐘。
全文完