2011年6月25日土曜日

重複組み合わせ(個数に制限がある場合) in Elisp

重複組合せ」より引用:

n 種のものから r 個取り出す組合せ(並べ方)の総数を求めよ。
ただし、同種のもの同士に区別はなく、それぞれ a1、a2、a3、・・・、an 個ずつ存在し、 n=a1+a2+a3+・・・+an≧r とする。

要するに、卑近な例を挙げれば「靴下どれでも5足で1,000円。ただしどれも在庫わずか」みたいな状況における、組み合わせの数でしょうか。こういう数のことを一般的に重複組合せ(ちょうふくくみあわせ、repeated combination)と呼び、個数の制限が無い場合は nHmm+n-1Cm = m+n-1Cn-1 で計算することができるとのこと(数学で習ったかもしれないけど、nHmなんて記号は覚えていない)。

以下では、個数制限ありの場合の問題を解くプログラムを書いてみた。n種類の物が並んでいる中から、(a)先頭の物を1個取り出す場合、(b)取り出さない場合、のそれぞれについて関数を再帰的に適用して、和を取っていけばいいだけ。

;;; lat: 個数のリスト(a1, a2, ..., an)
;;; r: 取り出す数
(defun solve (lat r)
  (cond
   ((null lat) 0)                       ;リストが空なので 0
   ((zerop r) 1)                        ;最後まで取り出したので 1
   (t (+ (if (<= 1 (car lat))           ;先頭の物を取り出す場合
             (solve
              (cons (1- (car lat)) (cdr lat)) ;先頭の物を1個消費
              (1- r))                   ;あと r - 1 個取り出す必要がある
           0)                           ;在庫切れで取り出せないから 0
         (solve (cdr lat) r)))))        ;先頭の物を取り出さない場合

;;; テスト
(solve '(1) 1)
=> 1

(solve '(1 1 1) 1)
=> 3
;; 以下の3通り
;; 1 0 0 
;; 0 1 0
;; 0 0 1

(solve '(2 2 2) 3)
=> 7
;; 以下の7通り
;; 2 1 0
;; 2 0 1
;; 1 2 0
;; 1 1 1
;; 1 0 2
;; 0 2 1
;; 0 1 2

(solve '(10 10 10) 10)
=> 66

最後の 66 の信憑性に関しては、こちらのページ「重複組み合わせ」に、3種類のものが10個ずつある中から10個とる場合の数は 66 だと書いてあったので間違いないと思われる。

ここでは自然に再帰で考えたけれども、漸化式を導出して動的計画法を使うと O(n×r) で解ける。具体的な方法は「Amazon.co.jp: プログラミングコンテストチャレンジブック: 秋葉 拓哉, 岩田 陽一, 北川 宜稔: 本」に書いてある。式の展開が多少面倒かもしれないが、紙と鉛筆を用意して臨めば問題ないレベルだった(数学嫌いでなければ)。

分割数 in Elisp

いわゆる「組み合わせ数学」の問題。整数 n を m 以下の整数の和で表す方法が何通りあるかを求めるというもの。考え方は「動的計画法」が参考になる。

単なる再帰で解く

(defun partition-num (n m)
  (cond
   ((or (= n 0) (= m 1)) 1)
   (t (+ (if (>= n m) (partition-num (- n m) m) 0)
         (partition-num n (1- m))))))
=> partition-num
(partition-num 3 3)
=> 3

(partition-num 4 3)
=> 4

(partition-num 9 9)
=> 30

動的計画法的に配列を使って解く

my-make-arrayという関数で多次元配列とそのアクセサみたいな関数を作っている。要するに、Emacs Lisp で多次元配列を使うのは面倒くさいということ。

;;; 多次元配列と、その要素を読み書きする関数(ref, set)を作る関数
;;; 内部的にはデータを1次元配列でもつ。クロージャと連想配列でできている。
(defun my-make-array (&rest dim)
  (lexical-let ((dim dim) (v (make-vector (apply '* dim) 0)))
    (flet ((index (args) (cond
                          ((null dim) 0)
                          (t (+ (* (apply '* (cdr dim)) (car idx))
                                (index (cdr dim) (cdr idx)))))))
      `((ref ,(lambda (&rest args) (elt v (index dim args))))
        (set ,(lambda (val &rest args) (aset v (index dim args) val)))))))
=> my-make-array

;;; 関数を適用する関数
(defun my-call (obj method &rest args)
  (apply (second (assoc method obj)) args))
=> my-call

;;; 配列を使って分割数を解く関数
(defun partition-num-dp (n m)
  (let ((dp (my-make-array (1+ n) (1+ m))))
    (loop for i from 0 to n do
          (loop for j from 1 to m do
                (my-call dp 'set i j
                         (if (or (<= j 1) (zerop i)) 1
                           (+ (my-call dp 'ref i (1- j))
                              (if (>= i j)
                                  (my-call dp 'ref (- i j) j)
                                0))))))
    (my-call dp 'ref n m)))

(partition-num-dp 4 3)
=> 4

(partition-num-dp 9 9)
=> 30

2011年6月22日水曜日

じゃんけんプログラム in Elisp

深夜になぜかEmacs Lispでジャンケンしようと思いたった…

いろいろな戦略を、たとえば相手の真似をする「反射戦略」や、グー・チョキ・パーを順番に出す「繰り返し戦略」などを、コールバック関数で渡す形で作る。いわゆるStrategyパターンと言ってよいだろうか。

メインの関数

;;; メインの関数
(defun janken (f1 f2)
  (switch-to-buffer "*janken-mode*")
  (let (first-val                       ;プレイヤー1の手
        second-val                      ;プレイヤー2の手
        result    ;勝敗を表す数(0: あいこ、1: P2の勝ち、2: P1の勝ち)
        history   ;手の履歴
        (vals '(g c p)))             ;手の選択肢(グー、チョキ、パー)
    (flet (                          ;ローカルな関数を定義
           (janken-prn-result ()     ;表示する関数
             (insert (format "%S : %S, %S\n"
                             first-val
                             second-val
                             (case result (0 'even) (1 'P2-win) (2 'P1-win)))))

           (janken-judge (v1 v2 f)                ;勝敗を判定する関数
             (if f (funcall f v1 v2)) ;事前処理。手の履歴を保存したり
             (destructuring-bind (x y)
                 (mapcar (lambda (v) (position v vals)) (list v1 v2))
               (% (+ (- x y) 3) 3))) ;三すくみの関係を表現する式

           (janken-loop ()                       ;ループする関数
              (while
                (zerop
                 (setf result
                       (janken-judge (setf first-val (funcall f1))
                                     (setf second-val (funcall f2))
                                     (lambda (x y) (setf history (cons (list x y) history))))))
                (janken-prn-result))
              (janken-prn-result)
              (message "%S" history) ;勝負がついた時点で*Messages*バッファに履歴を出力
              (when (y-or-n-p "Another game? ") ;継続するか問い合わせ
                (setf result nil)
                (janken-loop))))
      (janken-loop))))

なお、じゃんけんの「三すくみの関係」の表現方法はここに書いてあった:あるふぁ~版電算機部: 三すくみ。コードの(% (+ (- x y) 3) 3)に対応する。

戦略の関数

とりあえず5種類。

;手動入力
(defun janken-do-manually ()
  (intern (char-to-string (read-char "Press g or c or p: "))))

;ランダム戦略
(defun janken-do-randomly ()
  (elt vals (random (length vals))))

;繰り返し戦略(グー、チョキ、パーの順)
(defun janken-do-sequentially ()
  (elt vals (% (length history) 3)))

;pつまりパーを多めに出す戦略(確率 2/3)
(defun janken-do-weighted ()
  (let ((r (random 6)))
    (if (> r 2) (setf r 2))
    (elt vals r)))

;相手の最後の手を真似る戦略(反射戦略)
(defun janken-do-imitatively (firstp)
  (if history (funcall
               (lambda (x) (if firstp (second x) (first x)))
               (first history)) (janken-do-randomly)))

自由変数を含む形で書いているので(Elispなので動的スコープ)、メインの関数の中身と照らし合わせながら読まないと意味不明なはず。

実行例その1(手動 vs. ランダム戦略)

;次の式を評価する
(janken 'janken-do-manually 'janken-do-randomly)

;*janken-mode*バッファへの出力
g : c, P1-win
p : c, P2-win
p : g, P1-win
g : g, even
c : c, even
g : p, P2-win

実行例その2(パー多め戦略 vs. 反射戦略)

;次の式を評価する。反射戦略の関数は引数が要るので、ラムダ式を作る必要あり。
(janken 'janken-do-weighted (lambda () (janken-do-imitatively nil)))

;*janken-mode*バッファへの出力
p : g, P1-win
p : p, even
p : p, even
c : p, P1-win
p : c, P2-win
c : p, P1-win
p : c, P2-win
g : p, P2-win
p : g, P1-win

P2のほうはちゃんと真似している。これでよし。

…とまあこんな感じで、自動対戦させるプログラムくらいEmacs上ですぐに作れてしまうのであった。

2011年6月13日月曜日

最長増加部分列の長さ in Elisp

最長増加部分列(Longest Increasing Subsequence, LIS)の長さを求める問題を、Lisp風に解いただけ。

とりあえず手を動かして例題を解いてみる

これも動的計画法の有名な問題なので、漸化式で考える方法が普通かもしれない。が、自分の場合は例題を手作業で解いているときに頭の中で行っている計算をプログラムに書き直した感じなので、再帰的な解法になっている(空リストのときはゼロ、要素1個のリストの場合は1を返すという前提で開始すれば自然に再帰的な思考になる、はず)。

例題としては Wikipediaに載っていたもの (0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15) を考える。この数列のLISはユニークではなく (0, 2, 6, 9, 13, 15) と (0, 4, 6, 9, 11, 15) の二種類あり、LISの長さは 6。

;; 最後の要素から順番に、各要素を先頭とするLIS長を求めていく。
;; 計算を具体的に手でやってみると、次のような流れになる(リストに作用させる関数をfとおく)。
(f 15) => 1
(f 7 15) => 2              ;7 < 15 なので (f 15) + 1 => 1 + 1
(f 11 7 15) => 2           ;11 < 15 なので (f 15) + 1 => 1 + 1
(f 3 11 7 15) => 3         ;3より大きい11,7,15で始まるLISはそれぞれ長さが2,2,1。maxの2に1加えて3
(f 13 3 11 7 15) => 2      ;13 < 15 なので (f 15) + 1 => 1 + 1
(f 5 13 3 11 7 15) => 3    ;5より大きい13か11か7で始まるLISが最長で2。それに1加えて3
(f 9 5 13 3 11 7 15) => 3  ;9より大きい13か11で始まるLISが最長で2。それに1加えて3
(f 1 9 5 13 3 11 7 15) => 4    ;1より大きい9や5で始まるLISが最長で3。それに1加えて4
(f 14 1 9 5 13 3 11 7 15) => 2    ;14 < 15 なので (f 15) + 1 => 1 + 1
(f 6 14 1 9 5 13 3 11 7 15) => 4  ;6より大きい数のうち、9で始まるLISが最長で3。それに1加えて4
(f 10 6 14 ...略) => 3   ;10より大きい11,13のLISが最長で2。それに1加えて3
(f 2 10 6 ...略) => 5    ;2より大きい6のLISが最長で4。それに1加えて5
(f 12 2 10 ...略) => 3   ;12より大きい13, 14のLISが最長で4。それに1加えて3
(f 4 12 2 ...略) => 5    ;5より大きい数のうち、6で始まるLISが最長で4。それに1加えて5
(f 8 4 12 ...略) => 4    ;8より大きい12,10,9のLISが最長で3。それに1加えて4
(f 0 8 4 ...略) => 6     ;0より大きい数のうち、2で始まるLISが最長で5。それに1加えて6

Elispで

上で書いた解法をEmacs Lisp化する。使っている組み込み関数はcond, null, or, max, apply, mapcar, remove-if-notremove-if-notで全てremoveされてしまう場合に備え、orを使って空リストのリスト'(())を返すところが今ひとつだがしょうがない。また、補助関数のsublistsというのはたとえば(1 2 3)((1 2 3) (2 3) (3))のような形、つまり各要素を先頭とするリストのリストに変換するためのもの。これも別の形できれいに書けそうだが今は思いつかない。

(defun lis-length (lst)
  (cond
   ((null lst) 0)                       ;空リストなら0
   (t (1+                               ;サブ問題を解いて、その中の最大値に1を足す
       (apply 'max
              (mapcar
               'lis-length
               (or
                (remove-if-not          ;先頭の数より大きい数から始まるリストだけ抽出
                 (lambda (sublst) (> (car sublst) (car lst)))
                 (sublists (cdr lst)))  ;先頭の数を除いたリスト全て
                '(()))))))))
=> lis-length

(defun sublists (lst)
  (cond
   ((null lst) nil)
   (t (cons lst (sublists (cdr lst))))))
=> sublists

;;; テスト
(sublists '(1 2 3))
=> ((1 2 3) (2 3) (3))

(lis-length nil)
=> 0

(lis-length '(15))
=> 1

(lis-length '(7 15))
=> 2

(lis-length '(11 7 15))
=> 2

(lis-length '(0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15))
=> 6

あとはこの関数にメモ化の機能を追加してやれば、性能が上がって動的計画法で解いたのと同じになる(はず)。

2011年6月7日火曜日

個数制限つき部分和問題 in Elisp

プログラミングコンテストチャレンジブック演習「個数制限付き部分和問題」 - 日々精進に問題文が載っている。これも動的計画法の問題。

とりあえず性能を考えずに書いてみた場合

再帰で自然に書ける。

;;; 個数制限つき部分和問題
;;; - 数と個数のリスト ((数 個数) (数 個数) (数 個数) ...) 
;;; - 総和 K
;;; を与える。総和Kを満たす組合せが存在する場合は t を、その他の場合は nil を返す。
(defun solve (lst K)
  (cond
   ((zerop K) t)                        ;ちょうどゼロ => t
   ((< K 0) nil)                        ;マイナス => 解なし
   ((null lst) nil)                     ;もう追加できない => 解なし
   (t (or
       (if (< 0 (cadar lst))            ;数を追加できる(0個以上ある)
           (solve (decrease-first-num lst) (- K (caar lst))))
       (solve (cdr lst) K)))))          ;次の数を試す
=> solve
;;; 補助関数
(defun decrease-first-num (lst)
  (cons (list (caar lst) (1- ( cadar lst))) (cdr lst)))
=> decrease-first-num
;;; 補助関数の動作確認
(decrease-first-num '((1 4) (2 2) (3 3)))
=> ((1 3) (2 2) (3 3))
;;; 動作確認
(solve '((3 3) (5 2) (8 2)) 17)
=> t
(solve '((3 3) (5 2) (8 2)) 15)
=> nil

性能改善版(動的計画法あるいはメモ化)による

まず、関数solveの引数が取り得る範囲を減らすため、数の「個数」の部分を1ずつ減らしたりせず固定する。すると第一引数のlstを最大n種類になるので、第二引数のKと合わせて最大 n * K の組み合わせで全ての引数が表現でき、メモ化などで計算量を抑えることができる。その代わり関数の戻り値を t, nil の2値ではなく整数にして、リストの先頭の数が何個余るかを表すようにする(Kになる組み合わせがないときは -1 で表す)。

とりあえず数が1種類の場合と空リストの場合だけ考えると、

(defun solve* (lst K)
  (cond
   ((null lst) -1)                      ;空リストなら解なし
   ((zerop K) (second (car lst)))       ;K=0なら数を1個も消費しないので、個数をそのまま返す
   ((< K (first (car lst))) -1)         ;数が大き過ぎるなら解なし
   ((and (<= (first (car lst)) K)       ;Kから数を引くことができる?
         (< 0 (solve* lst (- K (first (car lst)))))) ;Kから数を1個引いた場合の問題を解いたら、まだ数が余る?
    (1- (solve* lst (- K (first (car lst)))))) ;余った数から1を引くと、解になる
   (t -1)))                             ;その他の場合は解けない
=> solve*

;;; テスト
(solve* '() 2)
=> -1

(solve* '((2 1)) 1)
=> -1

(solve* '((2 1)) 2)
=> 0 ;解あり

(solve* '((2 1)) 3)
=> -1

この発想になれたら、数が何種類もある場合も考慮する。先頭の要素をスキップしても解けたという場合を追加するだけなのであまり変わらないが、次のようになる。

(defun solve* (lst K)
  (cond
   ((null lst) -1)                      ;空リストなら解なし
   ((or (zerop K)                       ;K=0?
        (<= 0 (solve* (cdr lst) K)))    ;現在の数をスキップして解いたら、解けた?
    (second (car lst)))                 ;現在の数を1個も消費しないので、個数をそのまま返す
   ((< K (first (car lst))) -1)         ;数が大き過ぎるなら解なし
   ((and (<= (first (car lst)) K)       ;Kから数を引くことができる?
         (< 0 (solve* lst (- K (first (car lst)))))) ;Kから数を1個引いた場合の問題を解いたら、まだ数が余る?
    (1- (solve* lst (- K (first (car lst)))))) ;余った数から1を引くと、解になる
   (t -1)))                             ;その他の場合は解けない
=> solve*

;;; テスト
(solve* '((2 2) (4 2) (6 2) (8 2)) 40)
=> 0 ;解あり

(solve* '((2 2) (4 2) (6 2) (8 2)) 39)
=> -1

あとは、例によってこの再帰的関数をメモ化するとよい。

2011年6月6日月曜日

0-1ナップサック問題その2(重さWが巨大な場合) in Elisp

プログラミングコンテストチャレンジブック演習「ナップサック問題その2」 - 日々精進にも書いてある、動的計画法の問題。重さの総和Wが非常に大きい場合、重さを変数とする漸化式で n * W のテーブルを作ると計算量の制限を越えてしまうことがある。そういう場合は発想を変えて、価値を変数とする漸化式を考えよう、という話。すると、計算量は

n sum of v i equals 1 to n

と表されるようになり、巨大な W の値に依存しなくなる。

コード(Elisp)

配列は使わず、メモ化再帰で作る(メモ化の機能は後で付加する)。

solve-with-knapsack*がお膳立てをしてから、再帰的関数knapsack*を実行する。knapsack*が全て実行し終わると、再びsolve-with-knapsack*のほうが結果を整理して、最終的な解を出す。

;; 価値 v を満たす最小の重さを求める関数。
;; 解が無い場合は異常値(999999)を返す。
(defun knapsack* (lst v)
  (cond
   ;; v がゼロでないということは、解なし
   ((null lst) (if (zerop v) 0 999999))
   ;; 品物の価値が大き過ぎる場合
   ((< v (cadar lst)) (knapsack* (cdr lst) v))
   ;; 品物を入れない場合と入れる場合を試し、軽いほうを採用
   (t (min
       (knapsack* (cdr lst) v)
       (+ (caar lst) (knapsack* (cdr lst) (- v (cadar lst))))))))
knapsack*

;;; ちょっと確認
(knapsack* '((10 1) (15 2) (30 3)) 0)
=> 0
(knapsack* '((10 1) (15 2) (30 3)) 3)
=> 25
(knapsack* '((10 1) (15 2) (30 3)) 7)
=> 999999

;;; knapsack* を利用して問題を解く関数
(defun solve-with-knapsack* (lst w)
  (apply 'max                           ;価値が最大のものを選択
         (mapcar 'car                   ;価値だけを抽出
                 (remove-if             ;重過ぎる結果は削除
                  (lambda (x) (> (cdr x) w))
                  (mapcar
                   ;; (0 1 2 ... Sum vi)を満たす重さを求めて、(価値 . 重さ) のドット対にする
                   (lambda (v) (cons v (knapsack* lst v)))
                   ;; 0 から 価値の総和までの数列(0 1 2 ... Sum vi)
                   (loop
                    for i from 0 to (apply '+ (mapcar 'second lst))
                    collect i))))))
=> solve-with-knapsack*

;;; 確認
(solve-with-knapsack* '((50 1) (40 2) (30 3) (20 4) (10 5)) 60)
=> 12

実行時間をざっと調べる

;;; 計測するための関数
(defun profile (func &rest args)
  (let ((start (current-time)))
    (apply func args)
    (time-to-seconds (time-subtract (current-time) start))))
=> profile

;;; 重さと価値のリストを乱数で作り *lst* という変数に入れる
(defun gen-lst (n)
  (setf *lst*
        (loop for i from 1 to n collect (list (random 1000) (random 10)))))
=> gen-lst

;;; リストの要素数 n と実行時間の関係をざっと見る。Wは1000で固定。
;;; まだ再帰的関数をメモ化していないので、 2のn乗で増えていく。
(gen-lst 12)
=> ((801 6) (262 7) (923 5) (131 6) (221 0) (520 5) (221 9) (41 5) (209 6) (743 2) (791 2) (583 0))
(profile 'solve-with-knapsack* *lst* 1000)
=> 0.247282

(gen-lst 13)
=> ((977 2) (756 1) (256 3) (628 9) (743 2) (497 3) (551 0) (101 1) (51 4) (882 6) (181 4) (492 6) (649 4))
(profile 'solve-with-knapsack* *lst* 1000)
=> 0.451455

(gen-lst 14)
=> ((387 5) (154 5) (316 1) (803 0) (330 0) (450 9) (41 1) (257 5) (550 8) (575 6) (988 6) (69 8) (294 2) (330 7))
(profile 'solve-with-knapsack* *lst* 1000)
=> 1.229625

;;; 以下はメモ化した場合。nに対して鈍感になる。
(defun memo (fn name key test)
  "Return a memo-function of fn."
  (lexical-let ((_fn fn) (_key key) (table (make-hash-table :test test)))
    ;(setf (get name 'memo) table)
    #'(lambda (&rest args)
        (let ((k (funcall _key args)))
          (if (gethash k table)
              (gethash k table)
            (setf (gethash k table) (apply _fn args)))))))
=> memo

(defun* memoize (fn-name &key (key #'first) (test #'eql))
  "Replace fn-name's global definition with a memoized version."
  (setf (symbol-function fn-name)
        (memo (symbol-function fn-name) fn-name key test)))
=> memoize

(memoize 'knapsack* :key 'identity :test 'equal)
=> (lambda (&rest --cl-rest--) (apply (lambda (G6503 G6504 G6505 &rest args) (let ... ...)) (quote --table--) (quote --_key--) (quote --_fn--) --cl-rest--))

(gen-lst 14)
=> ((909 9) (234 2) (563 3) (890 8) (520 6) (951 8) (754 4) (810 7) (596 0) (334 1) (876 9) (799 8) (537 8) (735 8))
(profile 'solve-with-knapsack* *lst* 1000)
=> 0.011307

(gen-lst 15)
=> ((557 5) (365 2) (66 2) (794 7) (616 2) (9 1) (237 7) (364 2) (405 6) (976 5) (356 2) (558 7) (224 7) (876 0) (348 6))
(profile 'solve-with-knapsack* *lst* 1000)
=> 0.010048