Clojure での素因数分解プログラム

これ読んでたら Clojure で実装したくなったのでやってみた*1。

最終的なコードはこれ。

(ns prime-factorization.core
  (:gen-class))

(def prime-numbers
  ((fn f [x]
     (cons x
           (lazy-seq
            (f (first
                (drop-while
                 (fn [n]
                   (some #(zero? (mod n %))
                         (take-while #(<= (* % %) n) prime-numbers)))
                 (iterate inc (inc x))))))))
   2))

(defn prime-factorization [n]
  (if (< n 2)
    [n]
    (loop [dividend n
           prime-factors []
           primes (take-while #(<= % (Math/sqrt n))
                              prime-numbers)]
      (if-let [divisor (first primes)]
        (if (zero? (mod dividend divisor))
          (recur (/ dividend divisor)
                 (conj prime-factors divisor)
                 primes)
          (recur dividend
                 prime-factors
                 (rest primes)))
        (if (= dividend 1)
          prime-factors
          (conj prime-factors dividend))))))

(defn -main [num]
  (println (map str
                (prime-factorization (BigInteger. num)))))

遅延シーケンス使って素数リスト作ってるあたりがかなり大人気ない。この素数リストのコードはココから拝借。
Clojure - 職業プログラマは素数の夢を見るか? - Qiita

元の記事で書かれているような Timer は作らなかった(めんどくさい)。ので、適当にベンチマークだけかけた。
結果。

prime-factorization.core> (use 'criterium.core)
nil
prime-factorization.core> (bench (prime-factorization 8624871152))
WARNING: JVM argument TieredStopAtLevel=1 is active, and may lead to unexpected results as JIT C2 compiler may not be active. See http://www.slideshare.net/CharlesNutter/javaone-2012-jvm-jit-for-dummies.
WARNING: Final GC required 6.637604849087532 % of runtime
Evaluation count : 21840 in 60 samples of 364 calls.
             Execution time mean : 2.980871 ms
    Execution time std-deviation : 194.092642 µs
   Execution time lower quantile : 2.677770 ms ( 2.5%)
   Execution time upper quantile : 3.404867 ms (97.5%)
                   Overhead used : 10.291900 ns

Found 2 outliers in 60 samples (3.3333 %)
	low-severe	 1 (1.6667 %)
	low-mild	 1 (1.6667 %)
 Variance from outliers : 48.4644 % Variance is moderately inflated by outliers

平均的な実行時間約 3ms らしい。速い。

ちなみにマシンスペックこんな感じで、

この上に浮いてる VM の中で実行した結果です。

適当に書いたコードはココに置いておきますね。実行可能 jar を置いてあるので java 環境しかなくても動かせます*2。

% java -jar target/prime-factorization-0.1.0-SNAPSHOT-standalone.jar 1000
(2 2 2 5 5 5)

余談

最初はエラストテネスの篩をめちゃくちゃ素直に実装したんだけど、元の記事でテストされているような巨大な数を扱おうとすると range とかでリストを作ったタイミングで OutOfMemory とか StackOverflow とか出て死ぬ。まぁデカイからそのタイミングで死ぬことに何の疑問もないけど、遅延シーケンスでどうやって実装したらいいんだろうって悩んだ*3。

;;ボツコード
(defn sieve-of-eratosthenes [n]
  (assert (>= n 2))
  (loop [primes []
         numbers (range 2 (inc n))]
    (if-let [m (first numbers)]
      (if (<= m (inc (Math/sqrt n)))
        (recur (conj primes m)
               (filter #((complement zero?) (mod % m))
                       numbers))
        (apply conj primes numbers)))))

なんか僕が Clojure 書くと Clojure ぽくないコードになってしまう…。

理由は再帰を書いてしまっているからだと思う。同じエラストテネスの篩を実装した例と比べるとなんか僕のほうが圧倒的にひどい、つらい。
例えばこれとかすごく読みやすい。

(defn sieve [[xs ps]]
  (let [[p & more] xs]
    [(remove #(zero? (rem % p)) xs) (cons p ps)]))
 
(defn primes [n]
  (if (< n 2)
    []
    (->> [(range 2 (inc n)) nil]
         (iterate sieve)
         (drop-while #(< (ffirst %) (Math/sqrt n)))
         first
         (apply concat))))

無限遅延シーケンスによる解法は優雅で，しかも速い．しかし，遅い．とはいえ結局何とかなる． - tnoda-clojure

がんばろうって思った。

追記

少しだけ書き直したら読みやすくなった。

(defn- iter [[dividend prime-factors [divisor & rest-primes :as primes]]]
  (let [divisor (or divisor dividend)]
    (if (zero? (mod dividend divisor))
      [(/ dividend divisor) (conj prime-factors divisor) primes]
      [dividend prime-factors rest-primes])))

(defn prime-factorization' [n]
  (if (< n 2)
    [n]
    (let [primes (take-while #(<= % (Math/sqrt n)) prime-numbers)]
      (->> [n [] primes]
           (iterate iter)
           (drop-while #(not= (first %) 1))
           first
           second))))

ただ同じ条件でベンチマークしたら、コードが綺麗な方が遅くなった。読みやすさと速さを両立するのは難しい…。

*1:多分昔、 Scheme かなんかで同じようなの実装した気がするけど、、、

*2:誰が動かすんだ感

*3:結局、他人のコード借りてるわけだけど