Project Euler Problem 24

問題24: 0〜9の10個の数字で辞書順に順列を作る。100万番目のパターンは何か[=>問題文]
　順列生成ロジックはライブラリに入れたいですね。さてどう作ろう。
　順列生成というと、私の場合再帰関数で作るのが一番最初に思い浮かびます。でも再帰で順列を作ると再帰の一番深いところでパターン生成が完了するんですよね。schemeなら継続とか大域脱出とか使うかも知れませんが、Emacas Lisp ではあまりスマートにパターンを取り出す方法が思いつきません。そこで発想を転換。再帰の深いところからパターンを取り出すのではなく、パターンができたら関数を呼んでもらうことにしましょう。要するにコールバックです。

問題24の解答（再帰＋コールバック）

;; Emacs Lisp
(require 'cl)

(defun permutation-walk (array walker)
  (let ((continue t)
        (dest (subseq array 0)))
    (labels
        ((repeat (src dest-index)
                 (if (zerop (length src))
                     (unless (funcall walker dest)
                       (setq continue nil))
                   (do ((i 0 (1+ i))
                        (next-index (1+ dest-index)))
                       ((or (>= i (length src))
                            (not continue)))
                     (aset dest dest-index (aref src i))
                     (repeat (vconcat (subseq src 0 i)
                                      (subseq src (1+ i)))
                             next-index)))))
      (repeat array 0))))

(defun problem024-walk (n)
  (let* ((save)
         (walker (lambda (str)
                   (if (zerop (decf n))
                       (progn (setq save str) nil)
                     t))))
    (permutation-walk "0123456789" walker)
    save))

(problem024-walk 1000000)

　permutation-walk は第一引数に順列パターンの元となる配列、第二引数にコールバック関数をとります。再帰で順列パターンを作成し、一つできるたびにコールバック walker を評価します。walker が「真」を返す限りは順列を作り続け、nilを返すか、順列を最後まで作り終わると処理は終了します。
　problem024-walk では、n をカウントダウンし、ゼロになった時点で渡された配列を save して nil を返して順列生成を中断させています。
　walker に渡されるのは順列生成で使っているワーク配列なので、walker側でバッファの内容を変更したり保存したりする場合、コピーをとる必要があります。

STLのやりかた

　呼び出すたびに次の順列パターンを返すイテレータを作るにはどうすればいいんでしょうか。「次の順列を作る」順列生成ライブラリというと、真っ先に思い浮かぶのは STL(C++)の next_permutation です。というわけで gcc の テンプレートライブラリを見てみると、こんな感じのロジックでした。

// C++/STL
bool next_permutation(std::vector<int> &vec){

  std::vector<int>::iterator first = vec.begin();
  std::vector<int>::iterator last = vec.end();

  if (first+1 == last) return false;

  std::vector<int>::iterator i = last -1;

  for(;;){
    std::vector<int>::iterator ii = i;
    --i;
    if (*i < *ii){
      std::vector<int>::iterator j = last;
      while (!(*i < *--j)){}
      std::iter_swap(i, j);
      std::reverse(ii, last);
      return true;
    }
    if (i == first){
      std::reverse(first, last);
      return false;
    }
  }
}

　実際はもっと読みにくく汎用的ですが、わかりやすく改変しています。なぜこれでいいのか完全には理解できていませんが、i, ii, j の3つのインデックス(イテレータ)で配列内の要素を比較し、大小関係によって要素を swap したり reverseしたりして次のパターンを生成しているようです。これを Emacs Lispでアレンジしてみました。

別解（クロージャでイテレータを作る）

;; Emacs Lisp
(require 'cl)

(defun make-permutation-iterator (vect)
  ;; クロージャの環境
  (lexical-let*
      ((vect (subseq vect 0))         ; 順列の元となる配列
       (last (length vect))           ; 配列の長さ
       (handler)                      ; ハンドラクロージャ保存用
       ;; 配列内の2要素を入れ替える関数
       (swap-vect (lambda (i j)
                    (unless (= (aref vect i) (aref vect j))
                      (let ((tmp (aref vect i)))
                        (setf (aref vect i) (aref vect j))
                        (setf (aref vect j) tmp)))))
       ;; 配列内の ahead 以降の要素を反転する関数
       (reverse-vect (lambda (ahead)
                       (do ((a ahead (1+ a)) (b (1- last) (1- b)))
                           ((>= a b))
                         (funcall swap-vect a b))))
       ;; 順列の次のパターンを生成する関数
       (rest-handler (lambda ()
                       (let ((i (1- last)))
                         (do ((result t))
                             ((not (eq result t)) result)
                           (let ((ii i))
                             (decf i)
                             (cond
                              ((< (aref vect i) (aref vect ii))
                               (let ((j last))
                                 (while (>= (aref vect i) (aref vect (decf j))))
                                 (funcall swap-vect i j)
                                 (funcall reverse-vect ii)
                                 (setq result vect)))
                              ((zerop i)
                               (funcall reverse-vect 0)
                               (setq result nil)
                               (setq handler (lambda () nil)))))))))
       ;; 最初に一回だけ評価する関数（vect をそのまま返す）
       (first-handler (lambda ()
                        (if (< 1 last)
                            (setq handler rest-handler)
                          (setq handler (lambda () nil)))
                        vect)))         ; lexical-let* ここまで
    ;; クロージャ（イテレータ）を返す
    (if (> 1 last)
        (lambda () nil)
      (progn
        (setq handler first-handler)
        (lambda () (funcall handler))))))

(defun problem024-iter (n)
  (let ((iter (make-permutation-iterator "0123456789")))
    (while (< 0 (decf n))
      (funcall iter))
    (funcall iter)))

(problem024-iter 1000000)

　make-permutation-iterator はクロージャを返します。このクロージャは評価するたびに次の順列パターンを返します。最後までパターンを作り終えると、以降は nil を返し続けます。
　動的スコープの Emacs Lisp では本来クロージャに環境が紐付けされませんが、cl モジュールで提供される lexical-let や lexical-let* を使うことで環境を持ったクロージャを作ることができます。
　こちらも iter が返すのは ワーク配列なので変更や保存をするときにはコピーをとる必要があります。

どっちが速い？

　見たところどっこいどっこいかと思いましたが・・・

Function         Avg (sec)
===============  ===========
problem024-walk   13.359000
problem024-iter    6.979333

　ずいぶん差がつきました。約2倍です。walk は vconcat や subseq が頻繁に評価されるのでこれがボトルネックになっているのかもしれません。

Project Euler Problem 23

問題23: 2つの過剰数の和で書き表せない正の整数の総和を求めよ.(28123より大きい整数は2つの過剰数の和で表せることがわかっている）。[=>問題文]
　だんだん手ごたえのある問題になってきましたね。

問題23の解答

math-lib はこちら。

;; Emacs Lisp
(require 'math-lib) ; primes-sieve, get-divisors
(require 'cl)

(defun problem023 (n)
  (let* ((abundant-numbers  (make-vector n 0))
         (last-aset-index 1)
         (primes (primes-sieve (1+ (truncate (sqrt n)))))
         ;; 真の約数の和を計算する関数
         (sum-of-divisors
          (lambda (n)
            (- (apply #'+ (get-divisors n primes)) n)))
         ;; 2つの過剰数の和で表せるかチェックする関数
         (sum-of-abundants-number-p
          (lambda (n)
            (let* ((result nil)
                   (i 2)
                   (j last-aset-index)
                   (ni (aref abundant-numbers i))
                   (nj (aref abundant-numbers j)))
              (while (and (null result) (<= i j))
                (let ((sum (+ ni nj)))
                  (cond
                   ((< sum n) (setq ni (aref abundant-numbers (incf i))))
                   ((> sum n) (setq nj (aref abundant-numbers (decf j))))
                   (t (setq result t)))))
              result))))
    ;; メインループ
    (do ((k 1 (1+ k)) (ans 0))
        ((> k n) ans)
      (unless (funcall sum-of-abundants-number-p k)
        (incf ans k))
      ;; 過剰数ならば vector へ追加
      (if (< k (funcall sum-of-divisors k))
          (aset abundant-numbers (incf last-aset-index) k)))))

(problem023 28123)

　配列 abundant-numbers に順に過剰数を保存してゆきます。これを使って「2つの過剰数の和」で表せるかどうかの判定します。判定関数 sum-of-abundants-number-p では、インデックス i を abundant-numbers の先頭からインクリメント、j を(nより小さい過剰数の)末尾からデクリメントしてゆき、それらが指し示す2数の和を n に近づけていきます。(> i j) となったら該当する過剰数のペアは存在しません。なお判定処理のインデックス移動の都合から、abundant-numbers の最初の2要素は 0 固定で、インデックス2から過剰数を埋めています。
　今回 Core2 Quad 6600 と Pentium4(PrescottNorthwood)コアのCeleronで試してみました。結果は順に 3.8秒と 75.8秒。ハードウェアの差が大きく反映されました。 うそでした、多分Celeronマシンは byte-comile してなかったのかな。改めて計測しなおしました。

Function     Avg(sec)    CPU
==========   ========   =======================
problem023   3.766000   Core2 Quad 6600 2.4GHz
problem023   9.170000   Celeron (Pentium4 Northwood コア) 2.4GHz

　それでも同じクロック数であるにもかかわらず 約 2.4 倍の差が出ました。キャッシュ差かマルチコアの能力か。

Project Euler Problem 22

問題22: ファイルに保存されている全ての名前のスコアを計算しその合計を求めよ。[=>問題文]

問題22の解答

;; Emacs Lisp
(require 'cl)

(defun problem022 (file-path)
  (labels
      ((read-names () ;; ファイルから名前リスト読み込み
                   (let ((tmp)
                         (buffer (find-file file-path)))
                     (replace-string "\"" "")
                     (goto-char (point-min))
                     (replace-string "," " ")
                     (setq tmp (buffer-string))
                     (set-buffer-modified-p nil) ; 変更フラグクリア
                     (kill-buffer buffer)
                     (split-string tmp))))
    (let ((names (sort (read-names) #'string<)) ; 読み込んでソート
          (i 0) (code-base (1- ?A)))
      ;; スコア計算
      (apply
       #'+
       (cons 0.0
             (mapcar
              #'(lambda (name)
                  (* (incf i)
                     (reduce #'(lambda (acc c) (+ acc (- c code-base)))
                             name
                             :initial-value 0)))
              names))))))

(problem022 "names.txt")

　ファイルの内容を名前のリストにパースするためには。「カンマで分割」「ダブルクォートの除去」が必要です。Emacs Lisp でのやり方を調べたり試したりした結果、バッファに全部読み込んで replace-string でカレントバッファを編集し、lispで読み込みやすい形に整形してから文字列処理をする方法が一番簡単なようです。
　compile-defun でコンパイルすると、「replace-string は対話的操作で使う関数だぞ」と警告が出ますが、とりあえず問題なく機能しているようです。

リストにしてしまえばもっと楽（？）

　「lispで読みやすい形に整形」をもう一歩進めて、バッファ内文字列をリストの形に編集してみました。

(defun problem022 (file-path)
  (labels
      ((read-names ()
                   (let ((lis)
                         (buffer (find-file file-path)))
                     (insert "(")
                     (replace-string "," " ")
                     (goto-char (point-max))
                     (insert ")")
                     (goto-char (point-min))
                     (setq lis (read buffer))
                     (set-buffer-modified-p nil)
                     (kill-buffer buffer)
                     lis)))
    ;; 以下略

　バッファの最初に "(" を挿入して、","をスペースに置換。そしてバッファ最後に ")"を挿入。あとは read 一発でリストとして全データを読み込みます。

Project Euler 用 汎用計算ライブラリ

;;
;; math-lib.el
;; 利用するときは、このファイルをload-pathの通ったところに置いて
;; バイトコンパイルして、require してください。

(provide 'math-lib)
(require 'cl)

(defun primes-sieve (n)
  "エラトステネスのふるい (n未満の素数を列挙)"
  (let* ((vec (make-vector n nil))
         (vec-len (length vec)))
    (do ((i 3 (+ i 2)) (lis (list 2)))
        ((<= vec-len i) (reverse lis))
      (if (null (aref vec i))
          (push i lis))
      (do ((j i (+ j i)))
          ((<= vec-len j))
        (setf (aref vec j) t)))))

(defun factorize (n primes)
  "素因数分解"
  (let ((dest) (limit (1+ (truncate (sqrt n)))))
    (do ((primes primes (cdr primes))
         (continue t))
        ((or (null continue) (null primes)))
      (let* ((pn (car primes))
             (count (do ((count 0 (1+ count)))
                        ((< 0 (mod n pn)) count)
                      (setq n (/ n pn)))))
        (cond
         ((< 0 count)
          (push (cons pn count) dest))
          ((or (>= pn n) (> pn limit))
           (setq continue nil)))))
    (if (< 1 n)
        (push (cons n 1) dest))
    (reverse dest)))

(defun de-factorize (lis)
  "素因数リストから元の合成数復元"
  (reduce #'(lambda (acc cns)
              (* acc (expt (car cns) (cdr cns))))
          lis
          :initial-value 1))

(defun get-divisors (n primes)
  "約数リスト作成"
  (let ((factors (factorize n primes))
        (dest))
    (labels
        ((repeat (lis acc)
                 (if (null lis)
                     (push acc dest)
                   (destructuring-bind (n . count) (car lis)
                     (do ((i 0 (1+ i)))
                         ((< count i))
                       (repeat (cdr lis) acc)
                       (setq acc (* acc n)))))))
      (repeat factors 1))
    dest))

Project Euler Problem 21

問題21: 10000未満の友愛数の合計を求めよ。[=>問題文]
　素数リスト、素因数分解、約数、メモ化と、これまでに出てきた要素がいろいろ活用できる良い問題ですね。

問題21の解答

;; Emacs Lisp
(require 'cl)
(require 'math-lib)

(defun problem021 (n)
  (let ((primes (primes-sieve (1+ (truncate (sqrt n)))))
        (memo (make-vector n -1))
        (amicable-numbers))
    (do ((m 1 (1+ m)))
        ((= m n))
      (let ((sum (- (apply #'+ (get-divisors m primes)) m)))
        (aset memo m sum)
        (when (and (< sum n)
                   (not (= sum m))
                   (= (aref memo sum) m))
          (push m amicable-numbers)
          (push sum amicable-numbers))))
    (apply #'+ amicable-numbers)))

(problem021 10000)

　素数列挙等のライブラリを改めてまとめた math-lib ファイルを作りました。素数リストを必要とする関数は引数に素数リスト (primes)をとるようにしています。必要な素数の上限がわかる問題なので素数リストは、エラトステネスのふるい (primes-sieve) で作っています。
　整数 m の真の約数の和を memo のインデックス m 番目に保存して、友愛数チェックを行っています。

約数リスト作成関数

　今回作った約数リスト作成関数です。

(defun get-divisors (n primes)
  "約数リスト作成"
  (let ((factors (factorize n primes))
        (dest))
    (labels
        ((repeat (lis acc)
                 (if (null lis)
                     (push acc dest)
                   (destructuring-bind (n . count) (car lis)
                     (do ((i 0 (1+ i)))
                         ((< count i))
                       (repeat (cdr lis) acc)
                       (setq acc (* acc n)))))))
      (repeat factors 1))
    dest))

　素因数分解した後、順列を作る要領で再帰関数 repeat で約数を順次作成しています。作成される約数リストはソートされていません。また n 自身もリストに含まれます（真の約数リストではない）。
　前掲の problem021 では、この get-divisors の作成したリストの要素の和をとった後、対象の数自身を引くことで「真の約数の和」を得ています。

Project Euler 多倍長演算ライブラリ

　今後も使うかもしれないのでここにまとめておきます。
　加算、乗算、累乗、階乗が実装済みです。

多倍長演算ライブラリ

;;
;; bigint.el
;; 利用するときは、このファイルをload-pathの通ったところに置いて
;; バイトコンパイルして、require してください。
;;

(provide 'bigint)
(require 'cl)

(defun integer-to-bigint (n)
  "整数を多倍長整数へ"
  (if (< n 10000)
      (list n)
    (let ((stack) (m))
      (while (> n 0)
        (push (% n 10000) stack)
        (setq n (/ n 10000)))
      (reverse stack))))

(defun string-to-bigint (str)
  "文字列を多倍長整数へ"
  (let ((stack nil)
        (d 1) (c 0) (n 0))
    (dolist (c (reverse (string-to-list str)))
      (incf n (* (- c ?0) d))
      (when (< 1000 (setq d (* d 10)))
        (push n stack)
        (setq n 0 d 1)))
    (reverse (if (< 1 d)
                 (cons n stack)
               stack))))

(defun bigint-to-string (n)
  "多倍長整数を文字列へ"
  (if (= 1 (length n))
      (number-to-string (car n))
    (let ((n (reverse n)))
      (apply #'concat
             (cons
              (number-to-string (car n))
              (mapcar #'(lambda (x)
                          (format "%04d" x))
                      (cdr n)))))))

(defun bigint-length (n)
  "桁数を取得"
  (let ((m (reverse n)))
    (+ (length (number-to-string (car m)))
       (* 4 (length (cdr m))))))

(defun bigint+ (n1 n2)
  "加算"
  (let* ((long-n (if (< (length n1) (length n2)) n2 n1))
         (short-n (if (eq long-n n1) n2 n1))
         (carry 0)
         (lis (mapcar
               #'(lambda (x)
                   (let ((n (+ carry x (pop long-n))))
                     (setq carry (/ n 10000))
                     (% n 10000)))
               short-n)))
    (if (zerop carry)
        (append lis long-n)
      (let ((stack nil)
            (len-remain (length long-n)))
        (while (and (< 0 carry)
                    (<= 0 (decf len-remain)))
          (let ((n (+ carry (pop long-n))))
            (push (% n 10000) stack)
            (setq carry (/ n 10000))))
        (if (< 0 carry)
            (append lis (reverse stack) (list carry))
          (append lis (reverse stack) long-n))))))

(defun bigint* (n1 n2)
  "乗算"
  (let* ((long-n (if (< (length n1) (length n2)) n2 n1))
         (short-n (if (eq long-n n1) n2 n1))
         (result '(0))
         (a 0) (d nil))
    (dolist (a short-n result)
      (when (< 0 a)
        (let* ((carry 0)
               (lis (mapcar
                     #'(lambda (b)
                         (let ((n (+ (* a b) carry)))
                           (setq carry (/ n 10000))
                           (% n 10000)))
                     long-n)))
          (setq result
                (bigint+
                 result
                 (if (zerop carry)
                     (append d lis)
                   (append d lis (list carry)))))))
      (push 0 d))))

(defun bigint-expt (n m)
  "累乗計算。n:bigint, m:0以上の整数"
  (let ((result '(1)) (i 0) (cnt 0))
    (while (zerop (/ m 2))
      (incf cnt)
      (setq m (/ m 2)))
    (while (<= 0 (decf m))
      (setq result (bigint* result n)))
    (while (<= 0 (decf cnt))
      (setq result (bigint* result result)))
    result))

(defun bigint-fact (n)
  "階乗計算。n:0以上の整数"
  (do ((m 2 (1+ m))
       (result (string-to-bigint "1")))
      ((< n m) result)
    (setq result
          (bigint* result
                   (string-to-bigint
                    (number-to-string m))))))

2009/11/28 bigint-length 追加
2009/12/04 bigint* バグ修正、integer-to-bigint 追加

Project Euler Problem 20

問題20: 100! の各桁の数字の合計を求めよ。[=>問題文]

問題20の解答

;; Emacs Lisp
(require 'cl)
(require 'bigint)

(defun problem020 ()
  (let ((str (bigint-to-string
              (bigint-fact 100))))
    (apply #'+
           (mapcar #'(lambda (c) (- c ?0)) str))))

(problem020)

　bigint ライブラリに階乗関数を追加して終了です。