Project euler 21

pe_21.clj

;; 初めて解いたときのものです。  

;; 約数の和を求めるために、素因数分解してそれを利用してます。
;; 思いつきをそのまま実装しちゃってる感じです。

;; Problem 21  : 2011/5/5
;;"Elapsed time: 1260.253748 msecs"
;; prime list

;; まずは素数列を作ってしまいます。
(def prime-list (atom []))

(defn is-prime? [target]
  (loop [plist @prime-list]
    (let [check-num (first plist)]
      (cond (empty? plist) true
            (> check-num (Math/sqrt target)) true
            (zero? (rem target check-num)) false
            true (recur (rest plist))))))

(defn create-prime-list-under [n]
  (loop [target 2]
    (if (>= target n)
      true
      (if (is-prime? target)
        (do
          (reset! prime-list (conj @prime-list target))
          (recur (inc target)))
        (recur (inc target))))))

(create-prime-list-under 100000)

;; 素因数分解します。
;; 素数のリストの小さい方から順に割っていくのを再帰で書いてます。
;; 割れた場合はその素数は残します。数が1になったら終りで、割った素数のリストを返します。

;; factors 
(defn factors [n]
  ;; (foctors 12) => (2 2 3)
  (loop [factor-list ()
         target n
         prime-list @prime-list]
    (if (or (= target 1) (empty? prime-list))
      (reverse factor-list)
      (let [one-prime (first prime-list)]
        (if (zero? (rem target one-prime))
          (recur (cons one-prime factor-list) (/ target one-prime) prime-list)
          (recur factor-list target (rest prime-list)))))))

;; p^n (実はpがn個のリスト(p p ... p))について、p^n+p^(n-1)+p^(n-2)+...+p^0を計算します。
;; 乗の和っていうネーミングですが、よくないですね。

(defn multi-plus [lst]
  ;; multi-plus : (2 2 2) -> (+ 8 4 2 1)
  (if (empty? lst)
    1
    (+ (apply * lst) (multi-plus (rest lst)))))

;; 約数の和を求める関数です。ちょっとわかりにくい記述なのが気になります。
;; 約数に、自分自身も入ってしまうので、最後に引いています。

(defn sum-of-divisor [n]
  (let [ftr (factors n)]
    (- (reduce *
            (for [tgt (distinct ftr)]
              (multi-plus (filter #(= % tgt) ftr))))
       n)))

;; 回答を作る処理
(reduce +
  (filter (fn [x] (let [y (sum-of-divisor x)
                        z (sum-of-divisor y)]
                     (and (= x z) (not (= x y)))))
		(range 2 10000)))

pe_21_2.clj

;; これまでの成果を取り入れて、作りなおしてみました
;;必要な素数列の作成も込みで以下の時間ですので、倍くらい早くなっています。
;;"Elapsed time: 1552.285574 msecs"

;; 見習ってテスト入れてみました
(require '[clojure.test :as test])

;; 前に作った素数列の作成のためのツールたち
(defn make-is-prime? []
  (let [prm-list (atom [2])]
    (fn [target]
      (if (not-any? #(zero? (rem target %))
                    (take-while #(<= % (Math/sqrt target)) @prm-list))
        (do (swap! prm-list conj target)
            true)
        false))))

(defn create-primes [is-prime?]
  (filter is-prime? (cons 2 (iterate #(+ % 2) 3))))

;; クロージャーにすることで、factors関数が自前の素数列を持てるようにしました。
(let [primes (create-primes (make-is-prime?))]
  ;; factors 
  (defn factors [n]
    (loop [factor-list ()
           target n
           prime-list primes]
      (if (or (= target 1) (empty? prime-list))
        (reverse factor-list)
        (let [one-prime (first prime-list)]
          (if (zero? (rem target one-prime))
            (recur (cons one-prime factor-list) (/ target one-prime) prime-list)
            (recur factor-list target (rest prime-list))))))))


;; reducetionsを使うことで、とてもすっきり。
(defn multi-plus [lst]
  (reduce + 1 (reductions * lst)))

;; 最近は、->> マクロが使えるようになりました。
;; partition-by を使うことで、とてもすっきり。
(defn sum-of-divisor [n]
  (->> (factors n)
       (partition-by identity ,,)
       (map multi-plus ,,)
       (reduce * ,,)
       (+ (- n) ,,)))


;; tests
(test/is (= (factors 24) [2 2 2 3]))

(test/is (= (multi-plus [2 2 2]) (+ 8 4 2 1)))
(test/is (= (multi-plus [5]) (+ 5 1)))

(test/is (= (sum-of-divisor 12) (+ 1 2 3 4 6)))
(test/is (= (sum-of-divisor 7) (+ 1)))


;; project euler 21
(reduce +
  (filter (fn [x] (let [y (sum-of-divisor x)
                        z (sum-of-divisor y)]
                     (and (= x z) (not (= x y)))))
		(range 2 10000)))

pe_21_3.clj

;; 別解です。
;;"Elapsed time: 923.9876 msecs"
;;激速です。

;; 約数の和のリストを生成
(defn sigma1-list [^long size]
  (loop [i 2
         arr (long-array size 1)]
    (if (>= i (/ size 2))
      arr
      (recur (inc i)
             (loop [a arr
                    j (* i 2)]
               (if (>= j size)
                 a
                 (recur (do (aset a j (+ (aget a j) i))
                            a)
                        (+ j i))))))))

(test/is (= (seq (sigma1-list 20))
            (seq (long-array [1 1 1 1 3 1 6 1 7 4 8 1 16 1 10 9 15 1 21 1]))))

;;友愛数を抽出
(defn find-amicable [size arr]
  (for [i (range 1 size)
        :let [j (aget arr i)]
        :when (and (not= i j)
                   (<  j size)
                   (= (aget arr j) i))]
    i))

(test/is (= (find-amicable 300 (sigma1-list 300))
            [220 284]))


;; pe21回答用処理
(let [size 10000]
  (->> (sigma1-list size)
       (find-amicable size ,,)
       (reduce + ,,)))

pe_21_3_解説.md

問題

d(n)をn未満の真の約数の和としたとき、d(a)=b かつ d(b)=a であるようなa とbは友愛数と呼ばれます。 10000以下の全ての友愛数の和を求めなさい。

解説

pe_21_1.cljは以前解いたときのものです。pe_21_2.cljは改めて見てみて気に 入らないところを直したものです。

解きかた

対象となる数それぞれの約数の和を求めるわけですが、約数そのものを求めな くても約数の和を求めることができます。Wikipediaとか こことかに解説があり ますが、その数を素因数分解したものを使います。

素因数分解するために作った関数が、factorsです。 小さい素数から順に割れるだけ割って、次の素数を調べます。

次に、上の方法で必要な数列を作ります。 3^4であれば、(3^0+3^1+3^2+3^3+3^4)=(1+3+9+27+81) が欲しいのですが、こ れを作るのが、multi-plusです。

素因数分解した各成分についてそれぞれmulti-plusを作って掛けると、約数の 和になります。

ここまでできたら回答です。

範囲の数について、約数の和を求め、さらにその数の約数の和を求めて、元の数 と同じになれば友愛数です。

別解を作りました。

Elapsed time: 923.9876 msecs" 上々ですね。

約数の和のデータを作る

自分で思いついたのだか、人から聞いたのだかわすれましたが、ある範囲の数 について、その約数の和をすべて求めるには、こんな方法もあります。

以下は、1〜8までの数の約数を並べたものです。

1 : 1
2 : 1 2
3 : 1   3
4 : 1 2   4
5 : 1      5
6 : 1 2 3   6
7 : 1        7
8 : 1 2   4   8

こうやって並べてしまうとそれだけで分ってしまうのですが、縦に見ていくと、 2を約数に持つものは2の倍数で、3を約数に持つものは3の倍数です。

この性質を使えば、機械的に表を埋めていくことができます。この先は、5,10,15...と5の倍数のところに5を入れていく、6の倍数のところに...と進めていけばいいということです。

そして、これは、ほとんどエラトステネスの篩の処理と同じです。

これを使って作ったのがsigma1-list関数です。

中身は普通に二重のループを回しているだけですが、約数の和がわかればいいので、上記の表と違って該当の数を加算しています。

long-array が返ります。clojure的には、シーケンスにして返すのが筋でしょうが、 この後の処理を考えてこうしてあります。

find-amicable : 受け取ったarrayから友愛数を持つものを抽出する

• iでチェックしたjは、iに含めなければ高速化しそうですが、どうしたらいいかわからない
• 友愛数を求めるという問題であれば、ペアにして出力すべきでしょうが、今回は合計してしまうので無視。

おまけ

約数の和のリストを作るときに、iではなく、1を足すと、約数の個数が求まります。

ということで、パラメータで渡すようししてみました。 約数関数の添字と同じく、σx として、 σ0 は約数の個数、σ1は約数の和がでます。 ただし、一般的な理解と異なり、この処理ではNはNの約数に入っていません。

(defn sigmax-list [x ^long size]
  (loop [i 2
         arr (long-array size 1)]
    (if (>= i (/ size 2))
      arr
      (recur (inc i)
             (loop [a arr
                    j (* i 2)]
               (if (>= j size)
                 a
                 (recur (do (aset a j (+ (aget a j) (long (Math/pow i x))))
                            a)
                        (+ j i))))))))

;; 別解です。
;;"Elapsed time: 923.9876 msecs"
;;激速です。

このアルゴリズムと実装ならもっと速いはず，と思ってみてみると，find-amicable 関数が type hinting されていませんでした．

(defn find-amicable [size ^longs arr]
  (for [i (range 1 size)
        :let [j (aget arr i)]
        :when (and (not= i j)
                   (<  j size)
                   (= (aget arr j) i))]
    i))

これで <10ms になるはずです．あと，for マクロのところを loop-recur で書き直すか，areduce を使えば，さらに速くなりそうですが，これはこれでそこまでしなくてもと思います．↓はそこまでした版の例，

(defn- asum
  [^longs xs]
  (areduce xs i ret 0
           (+ ret (let [j (aget xs i)]
                    (if (and (not= i j)
                             (< j (alength xs))
                             (= i (aget xs j)))
                      i
                      0)))))

(let [size 10000]
  (-> size sigma1-list asum))

@kohyama さん

これはすばらしいアルゴリズム 感動です.

これは、前にも紹介した、Rich Hickeyさんのエラトステネスの篩の借用です。
このような、「必要であれば、mutableなデータでも躊躇無く使う」というパランスの取れた姿勢が、clojureを魅力的にしているんだと思います。
ランダムアクセス+破壊的代入が必要なので、シーケンスでは効率が悪いでしょうね。

しっかり reductions 使われてました.

これまで、「どんな使いかたがあるのかわからん！」と思っていた、reductionsですが、わかってみるといろいろ使えます。「知っている」のと「使える」の差は大きいです。

,,は、半分苦しまぎれですね。 condのソースのインデントはemacs任せなのと、condの条件を文は並べて書きたくないので、どれが条件でどれが文だか分らなくなるし、矢印マクロは、関数の引数が足りないみたいで気持悪いんです。

@tnoda さん

type hint の効果、大きいですねぇ。普段、よほど のことが無いかぎり、速度のことは気にしていないのですが、やってないことはできるようにならないのは、実務でやという程思い知らされていますので、今後はもうちょっと気を付けてみようと思いました。

素因数分解のmulti-plusも、約数の和をエラトステネスの篩的に一気に求めてしまうsigma1-listも、すごいです。すてきです。
私が23を解いているときに約数の和を求めるのに苦しんでいたのですが、とても速くなりそうなので、ぜひ使わせていただきます。

ひとつ、配列を埋めるloop-recurで、引数に毎回配列を渡しているのが省けそうにみえたのが気になってやってみたのですが、「副作用を終わらせたあとにrecur」するためのdoがちょっと気持ち悪いことになってしまいました。

(defn sigma1-list [^long size]
  (let [a (long-array size 1)]
    (loop [i 2]
      (if (>= i (/ size 2))
        a
        (do (loop [j (* i 2)]
              (if (>= j size)
                a
                (do (aset a j (+ (aget a j) i))
                  (recur
                    (+ j i)))))
          (recur (inc i))
          )))))

doの中は一刻も早く抜けたいという印象があるのですが、これだと最外loop抜けるまですべてdoの中で、しかも何重にもネストします。これならたぶん、元のほうがきれいですね……。今後ぜひ参考にします。

@plaster (doto (long-array 10000) do-sigma1) みたいにして使うんですけど，こんなのはいかがでしょう↓

(defn- do-sigma1
  [^longs a]
  (loop [i 1 j 1]
    (if (< i (alength a))
      (if (< j (alength a))
        (do
          (aset a j (+ i (aget a j)))
          (recur i (+ i j)))
        (recur (inc i) (inc i))))))

loop-head を一つにするのがポイントです．実行時間は N=10000 で <1ms です．Clojure 速い．

---

と，ここでみなさんにお聞きしたいのですが，ここで書いた do-sigma1 関数は副作用しかない関数です，みなさんは，

• 副作用しかない関数を書きますか?
• 副作用しかない関数を書くときの決め事（スタイル）は何かありますか?
  • 私の場合は名前が do-something のように do-prefixed なことが多いです．

主に扱っている副作用は

• TCPソケット通信 (サーバ側)
• シリアルポート通信 (センサや制御機器との通信)
• 上記二つの読み書きバッファとなる配列の読み書き
• オンメモリDB代わりの atom の読み書き

などです.
副作用した結果どうだったか返すところまで含めて一つの塊にする感じでやってます.
が, 副作用対象によって状況が変わるので, なかなか, スタイルのようなものが確立できずに居ます.
後, 通信相手があるものは単体テストできず悩んでます.

「配列と破壊を使った効率的な計算」は, なるべく避けてシーケンスでやってきましたが, 今回, 沢山勉強させていだいたので今後はもう少し積極的に使っていこうと思ってます.

副作用しかない関数は、

• 状態保存のため
• 情報保存のため

に必要と判断すれば躊躇せず使います。データパッシングに置き変えることはあまり考慮しません。
PEを解くときには探索空間をatomにして、それを操作するための関数を作ります。
経験は多くありませんが、実務では、セッション管理でいくつか作りました。

命名については、末尾に「!」を付けるのがClojureの作法だと思っていましたので、そうしています。
でも、 asetなど一部の関数や、java由来の関数はこの規則に沿っていないので、あまり統一感は無いですよね。

個人的には、do が付いているものは逐次処理のための関数で、その対象に副作用のある関数を呼ぶことを想定しているものだと思うので、ちょっと違和感があります。