pe-23-2.clj

;;; Project Euler Problem 23 solution
;;; http://projecteuler.net/problem=23

(use 'clojure.test)

(def limit 28123)

(defn gen-sigma1-list
  "original implementation by @ypsilon-takai: https://gist.github.com/4284814"
  [^long size]
  (loop [i 2
         a (long-array size)]
    (if (< (* 2 i) size)
      (recur (inc i)
             (loop [j (* 2 i)
                    a a]
               (if (< j size)
                 (recur (+ j i)
                        (do (aset a j (+ (aget a j) i))
                          a))
                 a)))
      a)))

(defn gen-divisors [n] (filter #(zero? (mod n %)) (range 1 n)))
(defn abundant? [n] (< n (apply + (gen-divisors n))))

(defn gen-abundants [size]
  (let [sigma1-list (gen-sigma1-list size)]
    (filter #(< % (aget sigma1-list %))
            (range 1 size))))

(is (not-any? abundant? (range 1 12)))
(is (abundant? 12))

(defn gen-sums
  "candidates-ascending 中の2要素のすべての組み合わせについて和を計算する(thanks to @kohyama)"
  [candidates-ascending]
  (for [x candidates-ascending
        y candidates-ascending
        :while (<= y x)
        :let [sum (+ x y)]
        :while (<= sum limit)
        ]
    sum))

;; 補集合X-Aの要素の総和は、全体集合Xの要素の総和から部分集合Aの要素の総和を引いたものに等しい
;; thanks to @tnoda
(defn solve []
  (let [abundants (filter abundant? (range 12 (inc limit)))
        abundant-sums (set (gen-sums abundants))]
    (- (apply + (range (inc limit)))
       (apply + abundant-sums))))

(defn solve2 []
  (let [abundants (gen-abundants (inc limit))
        abundant-sums (set (gen-sums abundants))]
    (- (apply + (range (inc limit)))
       (apply + abundant-sums))))

;; (is (= (solve) 4179871))

;pe-23.core=> (time (solve))
;"Elapsed time: 25569.034676 msecs"
;4179871
;pe-23.core=> (time (solve2))
;"Elapsed time: 1789.932635 msecs"
;4179871

pe-23.clj

;;; Project Euler Problem 23 solution
;;; http://projecteuler.net/problem=23

(use 'clojure.test)

(def limit 28123)

(defn gen-divisors [n] (filter #(zero? (mod n %)) (range 1 n)))
(defn abundant? [n] (< n (apply + (gen-divisors n))))

(is (not-any? abundant? (range 1 12)))
(is (abundant? 12))

(defn gen-sums
  "candidates-ascending 中の2要素のすべての組み合わせについて和を計算する"
  [candidates-ascending]
  (if (empty? candidates-ascending)
    nil
    (let [ [lhs & larger-candidates] candidates-ascending ]
      (concat
        (for [rhs candidates-ascending
              :let [sum (+ lhs rhs)]
              :while (<= sum limit)
              ]
          sum)
        (lazy-seq (gen-sums larger-candidates))))))

(defn solve []
  (let [abundants (filter abundant? (range 12 (inc limit)))
        abundant-sums (set (gen-sums abundants))]
    (apply + (remove #(contains? abundant-sums %)
                   (range 1 (inc limit))))))

;; (is (= (solve) 4179871))

解説

「過剰数の和で表せない数すべてを求める」ために、「過剰数の和で表せる数すべてを求める」ことをしています。
28123以下の正の整数だけを考えればよいので、それを全体集合と考えて補集合を求めます。

• 28123以下の過剰数すべてを昇順に求めます: (filter abundant? (range 12 (inc limit)))
  • 私の手元の環境での実行時間全体(約28秒)の8割(約23.5秒)がここです。
    • (time (count (filter abundant? (range 12 (inc limit))))) => "Elapsed time: 23426.675034 msecs" 6965
    • abundant? が n以下すべてで試し割りするnaive実装なので、もっといい方法がありそうな気がしています。
• そのうちの2要素の和すべてを求めます: gen-sums
  • 以下のようなささやかな高速化をしています:
    • 28123を超えるものは求めても使わないので、計算を打ち切っています: :while (<= sum limit)
    • lhs と rhs を入れ替えただけの組み合わせをそもそも探さないようにしています。
      • これをかんたんに実現するため、二重ループの代わりに再帰をしています。
        • 7000段近くの再帰なので、Stack overflowしないようにlazy-seq します。lazy-seqの入れどころがだんだん分かってきたかもしれません。
        • vectorなりにして添字アクセスならたぶん再帰要りませんが、添字つかうと間違えそうで、複雑さにもよりますが、こういうほうが私は好きです。
        • もしかして、全部をひとつのforで簡単に書けたりしないかなあ、というのもちょっと気になっています。
• あとは、要素の有無を高速判定できるように set に入れて、補集合を(remove #(contains? abundant-sums %) ...) で求めています。

「真の約数の和」については,
@ypsilon-takai さん担当の Problem 21 にも出てきますので, 参考にされると良いかと思います.

昇順の有限数列 (foo とします) 内の二要素の組み合わせ列挙は

(for [a foo
      b (take-while #(<= % a) foo)
      ...]
  ...)

では如何でしょうか.

@kohyamaさん
組み合わせ列挙、すばらしいです。この順に取り出せば自然に書けるんですね。
まったく思いつかなかったのですが、ほしかったのはまさにこれです。反映した版を追加しました。
ありがとうございます！

Problem 21 も見てみます。

@ypsilon-takaiさんの Problem 21 からほぼそのままもらってきたsigma1-listを使ったsolve2を書きました。
全体で10倍以上速くなってます。
過剰数生成部分だけに限ると、100倍オーダーの速度差です。さすがです……

pe-23.core=> (time (count (filter abundant? (range 12 (inc limit)))))
"Elapsed time: 23248.896077 msecs"
6965

pe-23.core=> (time (count (gen-abundants (inc limit))))
"Elapsed time: 272.861643 msecs"
6965

👍 分かりやすかったです．

あとは、要素の有無を高速判定できるように set に入れて、補集合を(remove #(contains? abundant-sums %) ...) で求めています。

これ私も binary search 代わりによくやります．実際 O(log n) でオーダー同じですしね．ただ，今回に限っては，

(- (apply + (range (inc limit))
   (apply + abundant-sums))

のほうが速かったかもしれません．

そのうちの2要素の和すべてを求めます

2 要素の和を全て求めるよりも，一つ一つ過剰数の和で書けるか確認していったほうが速いかもと思って書いたのが↓ です：

(ns tnoda.projecteuler.problem-23
    (:use clojure.test))

(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))))))

(def ^:private abundant-numbers
  (delay (->> (doto (long-array 20162)
                do-sigma1)
              (keep-indexed #(if (< (* 2 %) %2) %))
              (apply sorted-set))))

(with-test
  (defn- sum?
    [x]
    "Returns true if x can be written as the sum of two abundant-numbers."
    (->> (subseq @abundant-numbers <= (quot x 2))
         (some #(@abundant-numbers (- x %)))
         boolean))
  (is (= (sum? 12) false))
  (is (= (sum? 24) true)))

(defn solver
  []
  (->> (range 1 20162)
       (remove sum?)
       (apply +)))

で，これも 2 秒ほどとあまり速くありません．Clojure のコレクションではなく Java の配列を使うと 100ms 切れると思いますが，同じ配列を使うなら，@plaster の「2 要素の和を全て求める」方法のほうが速くなるはずです．

@kohyamaさんの提示のfor文の変数束縛の方法ですが、僕もこれが使えるのを知るまでは、whileで書いてました。
変数は順番に束縛されるのでこういつ使いかたができるんですね。

以前の自分の回答を見てみると、結果を出すのに2時間もかかっているらしい。 リベンジしないといけませんね。

@tnodaさん

ただ，今回に限っては，

(- (apply + (range (inc limit))
   (apply + abundant-sums))

のほうが速かったかもしれません．

なるほど！補集合を求めてから総和を計算するのと、総和をそれぞれ計算して差を取るの、まったく等価ですね。
この性質は見落としてました。コードもシンプルになりますし、素晴らしいとおもいます。反映します。