2012年11月10日土曜日

オイラーの贈物(1.2)とLispと順列

二項展開(binomial expansion)の展開後の各項の係数(二項係数(binomial coefficient))は 階乗(factorial)を利用して以下の用に書けます。

${}_n C _r \equiv \frac{n!}{r!(n - r)!}$
latexの数式
${}_n C _r \equiv \frac{n!}{r!(n - r)!}$
;; Common Lisp

(defun recursive-factorial (n)
  (check-type n (integer 0 *))
  (if (zerop n)
      1
      (* n (recursive-factorial (1- n)))))

(defun factorial (n)
  (check-type n (integer 0 *))
  (if (zerop n)
      1
      (loop
         :named loop
         :for i from 1 to n
         :for result = i then (* i result)
         :finally (return-from loop result))))

(defvar *memoized-factorial* (make-hash-table))
(setf (gethash 0 *memoized-factorial*) 1)
(setf (gethash 1 *memoized-factorial*) 1)
(defun memoized-factorial (n)
  (check-type n (integer 0 *))
  (let ((x (gethash n *memoized-factorial*)))
    (if x
        x
        (let ((y (* n (memoized-factorial (1- n)))))
          (setf (gethash n *memoized-factorial*) y)
          y))))

(defun binomial-coefficient (n r)
  (check-type n (integer 0 *))
  (check-type r (integer 0 *))
  (assert (<= r n))
  ;; n! / r!(n - r)!
  (/ (memoized-factorial n)
     (* (memoized-factorial r)
        (memoized-factorial (- n r)))))

この式はn個の中からr個を取る組み合わせ(combination)の数を求める式でもあります。

要素の順番も考慮する順列(permutation)の数の場合は、以下のようになります。

${}_n P _r \equiv \frac{n!}{(n - r)!}$
latexの数式
$ {}_n P _r \equiv \frac{n!}{(n - r)!}$
(defun permutation (n r)
  (check-type n (integer 0 *))
  (check-type r (integer 0 *))
  (assert (<= r n))
  ;; n! / (n - r)!
  (/ (memoized-factorial n)
     (memoized-factorial (- n r))))

多くのプログラミング言語では順列を生成するための機能が用意されているようです。

RubyのArrayクラス(array.cで定義)とC++のstd::next_permutationを使ってみます。

> p [1, 2, 3].permutation(2).to_a
[[1, 2], [1, 3], [2, 1], [2, 3], [3, 1], [3, 2]]

## nPr = 3P2 = 3! / (3-2)! = 3! = 6
> p [1, 2, 3].permutation(2).to_a.size
6

> p [1, 2, 3].combination(2).to_a
[[1, 2], [1, 3], [2, 3]]

## nCr = 3C2 = 3! / 2!(3-2)! = 3! / 2 = 3
> p [1, 2, 3].combination(2).to_a.size
3
#include <algorithm>
#include <cstdio>

void print_array(int arr[3]){
        printf("%d, %d, %d\n", arr[0], arr[1], arr[2]);
}

int main(void){
        int arr[3] = {1, 2, 3};
        std::sort(arr, arr+3);
        int count = 0;
        do {
                print_array(arr);
                count++;
        } while(std::next_permutation(arr, arr+3));

        printf("count = %d\n", count);
        return 0;
}

// 1, 2, 3
// 1, 3, 2
// 2, 1, 3
// 2, 3, 1
// 3, 1, 2
// 3, 2, 1
// count = 6

これらのアルゴリズムをCommon Lispで書いてみます。

なお、Common Lispには順列を操作するためのライブラリとして cl-permutation があります。 プログラム中で順列を生成したりしたい場合はわざわざ自作せずありがたく利用させて頂きましょう。

;; from Ruby (array.c)
(defun permute0 (n r p index used vals result)
  (dotimes (i n)
    (when (zerop (bit used i))
      ;; 順列のindex番目の要素として元の配列のi番の要素を選択
      (setf (svref p index) i)
      (if (< index (1- r))
          (progn
            ;; 選択済みの要素に対応するフラグをONにする
            (setf (bit used i) 1)
            (permute0 n r p (1+ index) used vals result)
            ;; 選択済みフラグをOFFにする
            (setf (bit used i) 0))
          (progn
            ;; 添字の配列から値の配列を作成して結果配列に追加
            (vector-push 
             (loop :for j :across p :collect (elt vals j))
             result))))))

(defun permute (seq r)
  (let* ((n (length seq))
         (result (make-array (permutation n r) :fill-pointer 0 :adjustable t))
         (used (make-array n :element-type 'bit :initial-element 0))
         (p (make-array r :element-type '(integer 0 *))))
    (permute0 n r p 0 used seq result)
    result))
> (permute '(1 2 3) 2)
#((1 2) (1 3) (2 1) (2 3) (3 1) (3 2))
;; from C++ (std::next_permutation)
(defun next-permutation (arr)
  (let ((len (length arr)))
    (unless (or (= len 0) (= len 1))
      (loop
         :for pos :from (1- len) :downto 1
         :when (< (aref arr (1- pos)) (aref arr pos))
         :do (progn
               (rotatef (aref arr (1- pos))
                        (aref arr
                              (position-if
                               #'(lambda (x) (< (aref arr (1- pos)) x))
                               arr
                               :from-end t)))
               (setf (subseq arr pos)
                     (nreverse
                      (make-array (- len pos)
                                  :displaced-to arr
                                  :displaced-index-offset pos)))
               (return-from next-permutation t))))))
> (defmacro do-while (test &body body)
          `(loop
              :do (progn ,@body)
              :while ,test))
> (let ((tmp (vector 0 1 2)))
          (do-while (next-permutation tmp)
            (print tmp)))
#(0 1 2) 
#(0 2 1) 
#(1 0 2) 
#(1 2 0) 
#(2 0 1) 
#(2 1 0)

他にも色々なアルゴリズムが存在するようです。 その中のひとつとして、階乗進数 を利用したプログラムを書いてみます。

(defun factoradic-permutation (n width)
  (let ((fact (make-array width :initial-element 1))
        (mantissa (make-array width :initial-element 0)))
    ;; 階乗を計算して配列に設定
    (loop
       :for i :from 1 :below width
       :for acc := i :then (* acc i)
       :do (setf (aref fact i) acc))
    ;; 階乗進数の仮数を計算して配列に設定
    (loop
       :for i :from (1- width) :downto 1
       :for acc := n :then (mod acc (aref fact (1+ i)))
       :do (setf (aref mantissa i)
                 (floor acc (aref fact i))))
    ;; 階乗進数を利用して順列を作成 (配列mantissaを使いまわす)
    (let ((tmp (loop :for i :from 0 :below width :collect i)))
      (loop
       :for i :from (1- width) :downto 0
       :for x := (aref mantissa i)
       :do (setf (aref mantissa i) (nth x tmp)
                 tmp (delete (nth x tmp) tmp))))
    (nreverse mantissa)))
> (dotimes (n 6)
>  (print (factoradic-permutation n 3)))
#(0 1 2) 
#(0 2 1) 
#(1 0 2) 
#(1 2 0) 
#(2 0 1) 
#(2 1 0)

2012年10月25日木曜日

オイラーの贈物とLisp(1.1.1)

オイラーの贈物 を読みながら、登場する数式などをCommon Lispで書いてみます。

1.1.1 自然数と素数(P6)より、有名なエラトステネスの篩と、 Wikipediaの 素数の項目 に載っていたウラムの螺旋を出力するコードです。 画像の作成には Vecto を利用しました。

(defun eratosthenes-sieve (n)
  "Return a sequence which indicates whether index is the prime number or not."
  ;; `n' is greater than or equal to 2
  (check-type n (integer 2 *))
  (let ((seq (make-array (1+ n) :initial-element t)))
    (setf (svref seq 0) nil) ; 0 is not a prime number.
    (setf (svref seq 1) nil) ; 1 is not a prime number.
    ;; The outer loop can stop at the square root of `n'.
    (loop :for i :from 2 :to (floor (sqrt n))
       :when (svref seq i)
       :do
       ;; The inner loop can start at the square of `i'.
       ;; (Multiples of `i' which are less than the square of `i' are already set to `nil'.)
       (loop :for j :from (* i 2) :to n :by i
          :do (setf (svref seq j) nil)))
    seq))


(ql:quickload "vecto")

(defun draw-uram-spiral (edge-length output-file-name &key (pixel 2))
  (assert (<= 1 edge-length))
  (let* ((limit (expt edge-length 2))
         (primes (eratosthenes-sieve (expt edge-length 2)))
         (picture-edge-length (* pixel edge-length)))
    (vecto:with-canvas
        (:width picture-edge-length :height picture-edge-length)
      (vecto:set-rgb-fill 1.0 1.0 1.0)
      (vecto:rectangle 0.0 0.0 picture-edge-length picture-edge-length)
      (vecto:fill-path)
      (vecto:set-rgb-fill 0.0 0.0 1.0)
      (let ((idx 2)
            (step 0)
            (x (floor edge-length 2))
            (y (floor edge-length 2)))
        (loop
           :while (<= idx limit)
           :do 
           (dotimes (_ (1+ (floor step 2)))
             (when (<= idx limit)
               (case (mod step 4)
                 (0 (incf x))
                 (1 (decf y))
                 (2 (decf x))
                 (3 (incf y)))
               (when (svref primes idx)
                 (vecto:rectangle (* x pixel) (* y pixel) pixel pixel)
                 (vecto:fill-path))
               (incf idx)))
           (incf step)))
      (vecto:save-png output-file-name))))
> (draw-uram-spiral 100 "uram.png")

2012年10月16日火曜日

coutへの出力をstringstreamへの出力に切り替える

標準出力(std::cout)への出力を別のストリーム(std::stringstream)への出力に切り替えてみます。

#include <iostream>
#include <sstream>

int main(void){
  std::stringstream ss;
  std::streambuf *backup = std::cout.rdbuf();

  // 標準出力への出力をstringstreamへの出力に切り替える
  std::cout.rdbuf(ss.rdbuf());

  std::cout << "A";

  std::cout.rdbuf(backup);
  std::cout << "stringstream = " << ss.str() << std::endl;
  return 0;
}

C++ならrdbuf、C言語ならfreopenなどを使えば良さそうです。

scratchバッファをorg-modeにする

org-babelがとても便利なので、scratchバッファをデフォルトでorg-modeにしてみます。

Emacsの終了時に自動的にファイルに保存されるようにしておけば、 ものぐさでorg-rememberを使いこなせない私でも日々org-modeを活用できそうです。

(defvar *scratch-file* "~/.scratch.org")

;; 初期化時の処理
(defun init-scratch-buffer ()
  (let ((buf (get-buffer "*scratch*")))
    (when buf
      (save-excursion
        (with-current-buffer buf
          (erase-buffer)
          (org-mode)
          (insert
           (format "* [%s]"
                   (format-time-string "%Y/%m/%d %H:%M:%S"))))))))
;; 終了時、バッファ削除時にバッファの内容を保存する処理
(defun save-scratch-buffer ()
  (let ((buf (get-buffer "*scratch*")))
    (when buf
      (save-excursion
        (with-current-buffer buf
          (append-to-file (point-min) (point-max) *scratch-file*))))))

(defun save-scratch-kill-emacs-hook ()
  (save-scratch-buffer))

(defun save-scratch-kill-buffer-hook ()
  (when (equal (current-buffer) (get-buffer "*scratch*"))
    (save-scratch-buffer)))

;; hook登録
(add-hook 'after-init-hook 'init-scratch-buffer)
(add-hook 'kill-emacs-hook 'save-scratch-kill-emacs-hook)
(add-hook 'kill-buffer-hook 'save-scratch-kill-buffer-hook)

2012年9月22日土曜日

Haskell入門書的クイックソート

Haskellの入門書に乗っていそうなクイックソートをClojureとCommon Lispで書いてみます。

- Clojure
;; defnやletで分配束縛ができます
;; group-byで関数を適用した結果の値によってグループ分けができます
;; ハッシュテーブル(map)は指定されたキーに対応する値を取得する関数にもなります
(defn qsort [cmp [piv & rst :as coll]]
  (if (empty? coll) []
    (#(concat (qsort cmp (%1 true)) [piv] (qsort cmp (%1 false)))
     (group-by #(boolean (cmp %1 piv)) rst))))

- Common Lisp
;; remove-if-not (filter)
(defun qsort-1 (cmp lst)
  (when lst
    (destructuring-bind (piv &rest rest) lst
      (flet ((f (x) (funcall cmp x piv)))
 (append (qsort-1 cmp (remove-if-not #'f rest))
  (list piv)
  ;; (remove-if-not (complement #'f) rest)
  (qsort-1 cmp (remove-if #'f rest)))))))

;; loopマクロ
(defun qsort-2 (cmp lst)
  (when lst
    (loop
       :with piv = (first lst)
       :for x in (rest lst)
       :if (funcall cmp x piv)
       :collect x into lesser
       :else
       :collect x into greater
       :finally (return
    (append (qsort-2 cmp lesser)
     (list piv)
     (qsort-2 cmp greater))))))

2012年9月10日月曜日

Clojure+leiningen+Apache POI

いつか業務でこっそり使うことを夢見てExcelをいじってみます。

1. leiningenでプロジェクト作成
 lein new poitest

2. 作成されたプロジェクトのディレクトリのproject.cljを編集
(defproject poitest "0.1.0-SNAPSHOT"
  :description "Apache POI Test"
  :dependencies [[org.clojure/clojure "1.4.0"]
                 [org.apache.poi/poi "3.8"]
                 [org.apache.poi/poi-ooxml "3.8"]]
  :main poitest.core)

3. 依存解決
 lein deps

4. コードを書く
(ns poitest.core
  (:gen-class))

(import '(org.apache.poi.xssf.usermodel
          XSSFSheet
          XSSFWorkbook
          XSSFRow
          XSSFCell)
        '(org.apache.poi.ss.usermodel
          WorkbookFactory))

(import '(java.io
          FileInputStream
          FileOutputStream))

(defn load-workbook [path]
  (-> path FileInputStream. WorkbookFactory/create))

(defn rows [^XSSFSheet  sheet]
  (let [nrows (.getPhysicalNumberOfRows sheet)]
    (letfn [(f [i]
              (if (<= nrows i)
                nil
                (cons (.getRow sheet i)
                      (lazy-seq (f (inc i))))))]
      (f 0))))

(defn create-9x9 [path]
  (let [wb (XSSFWorkbook.)
        sh (.createSheet wb)]
    (dorun
     (for [y (range 9)]
       (.createRow sh y)))
    (dorun
     (for [x (range 9) y (range 9)]
       (-> (.createCell (.getRow sh y) x)
           (.setCellValue (str (* (inc x) (inc y)))))))
    (.write wb (FileOutputStream. path))))


(defn -main
  "9x9を書き込んだExcelファイルを作成->ファイルを読み込み3列目の要素を表示"
  [& args]
  (create-9x9 "test.xlsx")
  (let [wb (load-workbook "test.xlsx")
        sh (.getSheetAt wb 0)]
    (dorun 
     (for [r (rows sh)]
       (println (-> (.getCell r 2) .getStringCellValue))))))

5. 実効したりコンパイルしたり
 lein run
 lein uberjar
※ 9/22 追記
こちらのほうがシンプルそう
(defn rows [sheet]
  (keep #(.getRow sheet %)
        (range
         (.getFirstRowNum sheet)
         (inc (.getLastRowNum sheet)))))

2012年9月9日日曜日

Clojure + Emacs環境を作る

いつのまにかすごく簡単にできるようになってました。

leiningen のインストール。Clojure本体もダウンロードしてくれるらしいです。
  1. leiningenのスクリプト(lein or lein.bat)をダウンロードしてくる
  2. 実効パスの通った場所に置いて、実効権限付与(chmod u+x)
  3. leinスクリプトを実行する。(lein self-install)(要 curl or wget)
 Emacs の設定。swank-clojureのGithubのページに、新しいユーザはnreplかRitzを使うと良い、と書いてあるので、nrepl.elを入れてみます。
  1. package.elのレポジトリとしてmarmaladeを登録する
  2. package-list-packagesでパッケージ一覧を表示する
  3. clojure-modeとnreplをインストールする
  4. M-x nrepl-jack-in  でemacsの中でleiningenのreplが起動

※9/9追記

clojure-modeのバッファでeldocを有効化

(add-hook 'clojure-mode-hook
   (lambda ()
     (nrepl-eldoc-enable-in-current-buffer)))