Project-Euler第69題

尤拉函式 · 發表 2019-01-28 23:21:22

摘要：大學的時候挺喜歡解Project Euler 上的題目的，尤其是它不在乎答題者使用哪一門程式語言，甚至還有很多參與者是使用pen&paper來解題的。去年開始重新開始做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多分鐘。

全文完

Project-Euler第69題

您可能也會喜歡…