Y箱

最大のパンデジタルな素数を求める問題。
ここでのパンデジタルな数というのは、1からnまでの全ての数を使った数ということなので、最大でも9桁(1から9)ということになります。

また、1から9まで全部足すと45になるので、9桁のパンデジタルな数は必ず3の倍数になっちゃいます。計算すると、3の倍数にならないのは、7桁と4桁(と1桁だけど、1は素数じゃない)だけ。


あとは、7桁以下の素数の大きなほうから順にパンデジタルかどうか確認していけばいいのでしょうけど、あえて、パンデジタルな数を作って素数かどうか判定する方式にした。
最初のコードは、あえてした割にはエレガントじゃないコードになっちゃった。

select-numsとnum-listx関数で、1からxまでの数の順列を作ろうとしてるんだけど、一般化できなかった。 マクロしかないような気がするけど、どうなんだろうと思って調べたら、contrib.combinatorics にあるみたい。permutations。
これを使ったら簡単だね。
短かいけど、そんなに速さは変らない。

;;
;; Problem 41  : 2011/6/9
;; "Elapsed time: 5.017956 msecs"

(use 'clojure.contrib.combinatorics)

(take 1
      (flatten 
       (for [n [7 4]]
         (filter is-prime? (map list-to-num 
                                (permutations (range n 0 -1)))))))
;;

初めにつくったのがこれ。

select-numsは、[col digit-list]を受けとって、colで指定された順にdigit-listの数字を取っていく関数。

(select-nums [0 0 0 0] [1 2 3 4]) -> [1 2 3 4]
(select-nums [3 2 1 0] [1 2 3 4]) -> [4 3 2 1]
(select-nums [1 1 1 0] [1 2 3 4]) -> [2 3 4 1]

;;
;; Problem 41  : 2011/6/9
;;"Elapsed time: 5.947124 msecs"

(use 'clojure.contrib.math)

(defn select-nums 
  ([col] (select-nums col [1 2 3 4 5 6 7 8 9]))
  ([col digit-list]
     (loop [num-list digit-list
            coll col
            res-list []]
       (if (empty? coll)
         res-list
         (let [tgt-num (nth num-list (first coll))]
           (recur (vec (remove #(= % tgt-num) num-list))
                  (rest coll)
                  (conj res-list tgt-num )))))))

(defn num-list7 []
  (for [a3 (range 7) a4 (range 6)
        a5 (range 5) a6 (range 4) a7 (range 3) a8 (range 2)]
    (select-nums [a3 a4 a5 a6 a7 a8 0]  [7 6 5 4 3 2 1])))

(defn list-to-num [digit-list]
  (apply + (map #(* %1 (expt 10 %2)) 
                (reverse digit-list) (iterate inc 0))))

(time (doall (take 1 (filter is-prime? 
                             (map list-to-num (num-list7))))))

;;

 
これまでもPEを解いていて、memoize関数によるメモ化が効いていないなぁと思うことが何度かあったのだけれど、そのままにしてあった。
でも、Problem 57 を解いているときにこれを解決しないと問題が解けないことが判明。
ちょっと調べてみたら、memoizeが思ったようなmemoizeをしていないことがわかって、僕の思うようなmemoizeをしてくれる関数を作る新しいマクロを作って解決。

なので、そのことについて。


たとえば、フィボナッチ数を求める関数。(実験のためにプリント文を仕込んである。)

;;
(defn fib [n]
  (do
     (println "-->" n)
     (if (<= n 1)
        n
        (+ (fib (dec n)) (fib (- n 2))))))
;;

メモ化する前に実行するとこんな感じ

user> (fib 4)
--> 4
--> 3
--> 2
--> 1
--> 0
--> 1
--> 2
--> 1
--> 0
3

行き着くところまでネストしていく様子が見えます。 さて、これを、memoize関数でメモ化します。

(def fib (memoize fib))

で、実行してみます。

user> (fib 4)
--> 4
--> 3
--> 2  ※
--> 1
--> 0
--> 1
--> 2  ★
--> 1
--> 0
3

続けて (fib 4) を実行すると、

user> (fib 4)
3

ちゃんとメモ化してるじゃん....じゃない!

ってことで、長くなるので
続きました。


最初のときの出力を見てほしい。★のところは、(fib 2) が呼ばれているんだけど、これは※のところですでに一度呼ばれているからメモ化が効くべきで、この先には行って欲しくない。ほかの (fib 1) (fib 0)も同様。

さらに(fib 5)をやってみるともっとがっかり。

user> (fib 5)
--> 5
--> 4
--> 3
--> 2
--> 1
…略…

おいおい。(fib 4) も (fib 3)もわかってるんだから、ぱっと出せよ…。

計算を始めちゃうと、それまでにメモされている情報にアクセスできていないってことですね。
これじゃあメモ化っていってもねぇ。だめじゃん。

ということで、この実装では、フィボナッチ数列のような漸化式型のテータの生成には使えないってことがわかりました。

素数判定みたいに、前の結果が関係ないものについてなら有効でしょうけど。

これじゃあ問題が解けないんで、ソースを見て考えた。 原因は、メモ化のためのメモ領域が関数の外にあることだと推測した。

memoizeのソースは下のもので、メモ化したい関数を引数「f」で与えるとメモ化された関数が返る。

(defn memoize [f]
  (let [mem (atom {})]
    (fn [& args]                           ;; あ
      (if-let [e (find @mem args)]
        (val e)
        (let [ret (apply f args)]          ;; い
          (swap! mem assoc args ret)  
          ret)))))

memoizeは「あ」のところで新しい関数をつくって、メモ領域の処理と「f」をラップして返している。
さて、再帰は「い」の「f」の中で行なわれているわけだけど、その中で呼びだされるのはあくまで「f」であって、memoizeで作成された関数ではないわけ。「f」からこいつを呼べればいいのだけれど、そんなことは不可能。

解決するには、メモ化したい関数とメモ領域を同じクロージャに入れればよい。

ということで、どっかを探せば誰かが作ってるとは思うけど、メモ化した関数を定義するマクロを作ってみた。実験のために、メモを消去する関数とメモの中身を表示する関数も一緒につくるようにした。

;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; create memoised func

(defmacro defn-memo 
  "Creates memoised function, and accessor of the memo which named -data,
   memo deletion function named -clear"
  [name arg-list & body]
  `(let [mem# (atom {})]

     (defn ~(symbol (str (str name) "-data")) []
       @mem#)

     (defn ~(symbol (str (str name) "-clear")) []
       (reset! mem# {}))

     (defn ~name ~arg-list
       (if-let [e# (find @mem# ~arg-list)]
          (val e#)
          (let [ret# ~@body]
              (swap! mem# assoc ~arg-list ret#)
              ret#))))))
;;

最後のやつが関数の実体。
このように、関数が呼び出されたあとにメモの中を見るようにしておけば、再帰的に呼ばれたときにも、過去の情報にアクセスできる。

使いかたは普通のdefnと同じ。単にdefnをdefn-memoに書き替える。

;;
(defn-memo mem-fib [n]
  (do
     (println "-->" n)
     (if (<= n 1)
        n
        (+ (mem-fib (dec n)) (mem-fib (- n 2))))))

使ってみる

user> (mem-fib 4)
--> 4
--> 3
--> 2
--> 1
--> 0
3

変っていないようだけど、メモ化してないバージョンと比べると、
 

user> (fib 4)
--> 4
--> 3
--> 2
--> 1
--> 0
--> 1 ★
--> 2
--> 1
--> 0
3

mem-fibは★以降の計算で前の値にヒットしているからその先を計算していないことがわかる。 さらに(fib 5) とかやってみると、 差は歴然。

user> (mem-fib 5)
--> 5
5

user> (fib 5)
--> 5
--> 4
--> 3
--> 2
--> 1
--> 0
--> 1
--> 2
--> 1
--> 0
--> 3
--> 2
--> 1
--> 0
--> 1
5

めでたしめでたし。


メモの中身は、mem-fib-data でアクセスできます。 atomなので、その気になれば変更できます。

user> (mem-fib-data)
{[5] 5, [4] 3, [3] 2, [2] 1, [0] 0, [1] 1}

消すこともできます。

user> (mem-fib-clear)
{}
user> (mem-fib-data)
{}

消しちゃうと、

user> (mem-fib 5)
--> 5
--> 4
--> 3
--> 2
--> 1
--> 0
5

計算しなおす。

ただ、メモ化したとはいえ、未計算のところは計算しなくちゃならなくて、定義が末尾再帰になってないから、このまま大きな数を求めようとするとスタックを食いつぶしてエラーになったりします。

(mem-fib 10000) ;; やっちゃだめ

こういう場合は、

(last (map mem-fib (range 10001)))

みたいにやる。(printlnを取っておかないと...)

user> (time (last (map mem-fib (range 10001))))
"Elapsed time: 60.550712 msecs"
33644764876431783266 … 7366875

たぶん、memoizeしたfibで同じことをやっても思ったようにはうごかない。

(last (map fib (range 10001)))

データはatomだから、pmapでも大丈夫。

user> (time (last (pmap mem-fib (range 10001))))
"Elapsed time: 227.005108 msecs"
33644764876431783266 …

でも、早くない。前に書いたように、1つずつ別プロセスでやってしまうとオーバーヘッドが増えちゃってだめ。

ここのところ忙しくてすすんでない。

1ケタが9こで 9ケタ分
2けたが90こで180ケタ分
3けたが900こで2700ケタ分
っていう具合になっていることを利用して求めようかと思ったけど、めんどくさくなって、普通に解いちゃった。

;;
;; Problem 40  : 2011/6/5
;; "Elapsed time: 6926.687051 msecs"
(use 'clojure.contrib.math)

(defn fractional-part-of-irrational []
  (mapcat #(seq (str %)) (iterate inc 1)))

(loop [seq (map list (fractional-part-of-irrational) (iterate inc 1))
       tgt (map #(expt 10 %) (range 7))
       res 1]
  (if (empty? tgt)
    res
    (let [[digit num] (first seq)]
      (if (= num (first tgt))
        (recur (rest seq) (rest tgt) (* (Character/digit digit 10) res))
        (recur (rest seq) tgt res)))))

;;

Project Euler を解いているときに、素数生成とかの細かい関数を作っているのだけれど、10個以上になってきたので、まとめることにした。
あらためて見てみると、変だったり、もう少し汎用的にしたほうがよかったりするものが目について、直し始めちゃって泥沼化寸前。

そのなかで、同じような計算をたくさんしているところで気になっていたのが、CPUの使用率が50%を超えないこと。mapなんかは自動的に並列処理してくれるのだと思っていたからちょっと意外だった。黒猫本(と呼ばれているかどうかは知らない)にも、STMとかの機能の説明はあるけど、並列化するための方法については書かれていなかった。「自分でJavaのThreadを呼ばなければならんのかい!」と思っていたら、pmapというのを見つけた。

で、実際どうなのかを調べてみた。


なんか適当に時間のかかる処理ということで、適当に2のX乗をY個足す処理にした。

mapでやるとこんな感じ。

(let [nums 10000]
  (time 
   (reduce + (map (fn [x] (rem (expt 2 x) 100)) (repeat 10000 nums)))))

"Elapsed time: 5561.606014 msecs"

5.6秒くらい。 このときのCPUの使用率は50%。


mapをpmapに変えると、

(let [nums 10000]
  (time 
   (reduce + (pmap (fn [x] (rem (expt 2 x) 100)) (repeat 10000 nums)))))

"Elapsed time: 3480.631756 msecs"

お、速い。 3.5秒。 36%高速化。

タスクマネージャで見ると、ちゃんと100%使用率になる。
あと、そこでjavaを片方のCPUにだけ割りあてると、もとと同じく、5秒半ぐらいかかる。


さて、これで安心してはいけない。
ここの記事を見たりすると、ただ単に1つ1つ並列化すればいいというものではないらしい。
そりゃあ当然だ。細切れにすれば、その分Over headが増えるわけだからねぇ。
どんな風に分けたらいいのかなぁ。

僕のPCはデュアルコアだから、2つに分けてみると、

(let [nums 10000]
  (time 
   (reduce + (pmap #(reduce + (map (fn [x] (rem (expt 2 x) 100)) %))
                   (partition 5000 (repeat 10000 nums))))))

"Elapsed time: 3092.340227 msecs"

3.1秒。 44%高速化。 8ポイント向上。


ほかも適当にやってみると、
10分割 "Elapsed time: 3064.371031 msecs"
100分割 "Elapsed time: 3099.338324 msecs"
1000分割 "Elapsed time: 3255.111753 msecs"

今回は10分割が一番速いという結果。


並列化のための関数は、他にもfutureとかいろいろあるみたいなので、これからちょっとずつつついていこうと思っとります。

循環させても素数のままの素数の問題。

循環させた数のリストを作る関数を作った。0の扱いでひっかかったりして、結局文字列化して操作する方式になった。

あとは全数検査
途中に偶数とか5とかが入っていたら対象からはずすとかしたらちょっと速くなるかも。

;;
;; Problem 35  : 2011/5/19
;; "Elapsed time: 18328.781172 msecs"

(def *prime-list* (atom []))
(create-prime-list-under 1000000)

;;get all ratate num
(require '[clojure.contrib.string :as ccstr])

(defn rotate-num-str [num-str]
    (str (ccstr/drop (dec (.length num-str)) num-str)
         (ccstr/butlast 1 num-str) ))

(defn get-all-rotate-num [n]
  (loop [new-num (rotate-num-str (str n))
         res-list [n]]
      (if (= new-num (str n))
        res-list
        (recur  (rotate-num-str new-num) 
                (conj res-list (Integer/valueOf new-num))))))

(count 
 (map first
      (filter #(every? true? (map is-prime? %))
              (map get-all-rotate-num @*prime-list*))))
;;