読者です 読者をやめる 読者になる 読者になる

BWT : bzip2 : 修正版

common lisp sbcl algorithm

前回書いたBWT関数はかなり非効率だったので、その修正版を作成した。
修正点は主に以下の三つ。

  1. 対象を文字列からバイト列に変更
  2. 列のローテートは実際には行わない。代わりに、各列の先頭位置を保持する変数を用意し、その値をずらすことで循環列を表現する。
  3. 列末尾への終端要素の追加をなくした。添字が列の長さに達したら、そこが終端位置(暗黙の終端要素)。
;;;;;;;;;;;;;;;
;;;; パッケージ
(defpackage :bwt
  (:use :common-lisp)
  (:export bwt
           bwt-inverse
           to-octets))
(in-package :bwt)

;;;;;;;;;
;;;; 宣言
(declaim (optimize (speed 3) (debug 0) (safety 0) (compilation-speed 0)))

;;;;;;;;;;;
;;;; 型定義
(deftype octet ()         '(unsigned-byte 8))
(deftype simple-octets () '(simple-array octet))
(deftype array-index ()   '(integer 0 #.array-total-size-limit)) ; 配列の添字型

;;;;;;;;;;;
;;;; 構造体: 変換対象となるバイト列と、その列の要素の先頭位置を保持する
(defstruct (bwt-octets (:constructor make-bwt-octets (octets &optional (start 0))))
  (octets #() :type simple-octets)
  (start   0  :type array-index))

(defun to-octets (bwt-octets)
  (bwt-octets-octets bwt-octets))

;; #<1 0 10 ...>形式で、bwt-octetsを表示する
;; 先頭位置の場所は、#<1^0 10...>などのように、'^'で表す。(例は、startが1の場合)
(defmethod print-object ((o bwt-octets) stream)
  (print-unreadable-object (o stream)
    (with-slots (start octets) o
      (loop FOR i FROM 0 BELOW (min (length octets) 200)
        DO
          (if (= i start)
          (princ #\^ stream)
        (unless (zerop i)
          (princ #\Space stream)))
      (princ (aref octets i) stream)
    FINALLY
      (when (= i start)
        (princ #\^ stream))
      (when (< i (length octets))
        (princ "..." stream))))))
> (make-bwt-octets (sb-ext:string-to-octets "bwt用の構造体"))
--> #<^98 119 116 231 148 168 227 129 174 230 167 139 233 128 160 228 189 147>

;; 循環列を作成
(defun gen-all-rotation (bwt-octets)
  (let* ((octets  (bwt-octets-octets bwt-octets))
         (len     (length octets))
         (rotated (make-array (1+ len) :element-type 'bwt-octets
                                       :initial-element bwt-octets)))
    ;; startが異なるそれぞれのbwt-octetsを作成する
    (dotimes (i len rotated)
      (setf (aref rotated (1+ i)) (make-bwt-octets octets (1+ i))))))

> (gen-all-rotation (make-bwt-octets (sb-ext:string-to-octets "abcde")))
--> #(#<^97 98 99 100 101> 
      #<97^98 99 100 101> 
      #<97 98^99 100 101>
      #<97 98 99^100 101> 
      #<97 98 99 100^101> 
      #<97 98 99 100 101^>)

;; bwt-octets同士の比較
(defun bwt-octets< (lft rgt)
  (let* ((octets (bwt-octets-octets lft))
         (len    (length octets))
         (start1 (bwt-octets-start  lft))
         (start2 (bwt-octets-start  rgt)))
    (loop FOR i = start1 THEN (1+ i)
          FOR j = start2 THEN (1+ j) DO
      (cond ((= i len)  #|終端位置|#             (return t))
            ((= j len)  #|終端位置|#             (return nil))
            ((< (aref octets i) (aref octets j)) (return t))
            ((> (aref octets i) (aref octets j)) (return nil))))))

;; bwt関数
(defun bwt (octets &aux (len (length octets)))
  (declare (simple-octets octets))
  (let ((bwt-octets (make-bwt-octets octets)))
    (flet ((last-octet (bos) (aref (bwt-octets-octets bos) (1- (bwt-octets-start bos))))
           (end-octet? (bos) (zerop (bwt-octets-start bos))))
      (loop WITH new-start  = 0
            WITH new-octets = (make-array len :element-type 'octet)
                                                 ;; XXX: ここでのソートが最大のボトルネック
            FOR bos ACROSS (the (simple-array t) (sort (gen-all-rotation bwt-octets) #'bwt-octets<))
            FOR i OF-TYPE fixnum = 0 THEN (1+ i)
        DO 
          (if (end-octet? bos)
              ;; 末尾の要素==(暗黙の)終端要素の場合
              ;; この位置が新しい先頭位置となる
              (setf new-start i
                    i         (1- i))
            ;; それ以外は、普通に列に追加する
            (setf (aref new-octets i) (last-octet bos)))
        FINALLY
          (return (make-bwt-octets new-octets new-start))))))

;; bwt-inverse関数用のソート関数
;; (funcall key 要素)が0〜255を返す配列の要素を、わりかし高速かつ安定にソートする
(defun bucket-sort (key idx-ary)
  (declare (function key)
           ((simple-array array-index) idx-ary))
  (let ((bucket (make-array #x100 :initial-element nil))
        (sorted (make-array (length idx-ary) :element-type 'array-index :initial-element 0)))
    (loop FOR i ACROSS idx-ary DO
      (push i (aref bucket (funcall key i))))  ; (funcall key i)の結果は、0〜255の範囲内
    
    (let ((i -1))
      (declare (fixnum i))
      (loop FOR list ACROSS bucket DO
        (dolist (e (nreverse list))
          (setf (aref sorted (incf i)) e))))
    sorted))

;; bwt逆変換関数
(defun bwt-inverse (bwt-octets)
  (let* ((octets   (bwt-octets-octets bwt-octets))
         (len      (length octets))
         (start    (bwt-octets-start bwt-octets))
         (heads    (make-array len :element-type 'array-index :initial-element 0))
         (original (make-array len :element-type 'octet  :initial-element 0)))
    (declare ((simple-array array-index) heads))

    (dotimes (i len)
      (setf (aref heads i) i))
    (setf heads (bucket-sort (lambda (i) (aref octets i)) heads))

    (loop FOR i FROM 0 BELOW len
                                               ; ↓ 明示的な終端要素をなくしたため、終端要素位置(start)の前方か後方かで、添字を調整する必要がある
          FOR p = (aref heads (1- start)) THEN (aref heads (if (< p start) (1- p) p))
          DO 
            (setf (aref original i) (aref octets p)))
    original))

実行例。
参照: read-binary-file

;; 夏目漱石の『こころ』(UTF-8)に対する結果
;; 見事に13が並んでいる
> (bwt:bwt (read-binary-file "/path/to/kokoro"))
--> #<10 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 
      13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 
      13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 
      13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 
      13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 
      13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13...>

;; 逆変換も問題無し
> (mismatch (bwt:bwt-inverse *) (read-binary-file "/path/to/kokoro"))
--> NIL

まだ、bwt関数内でのソート処理に結構時間が掛かっており改善の余地があるが、取り合えずbzip2を実装する上で差し支えがない程度のものはできた。
ソートに関しては、そのうち適切なアルゴリズムのものを実装するかもしれないが、当面はこのままで進めることにする。