Common Lisp 入門 の番外編です。今回はデータ圧縮法のひとつである「レンジコーダ」を取り上げます。なお、このドキュメントは拙作のページ Memorandum で説明した「Lisp でレンジコーダ」をまとめたものです。内容は重複していますが、ご了承くださいませ。
「レンジコーダ (RnageCoder) 」は 1998 年に Michael Schindler が発表し、「高性能、高速、特許フリー」の方法として注目を集めるようになりました。Michael Schindler のレンジコーダは計算の途中で「桁上がり」が発生します。ところが、ロシアの Dmitry Subbotin が発表した「桁上げのないレンジコーダ」は、その名のごとく桁上がりが発生しません。現在、レンジコーダは主にこの 2 種類の形式が存在するようです。
レンジコーダは原理的には算術符号と同じ方法です。性能は算術符号に比べるとわずかに劣化しますが、実現方法はとても簡単で実行速度も高速です。もちろん、ハフマン符号よりも高性能です。今回はレンジコーダのプログラムを Common Lisp で作成しましょう。これから作成するプログラムは学習用なので実用性はありませんが、実際にプログラムを作ることでレンジコーダの理解は深まると思います。
最初に、レンジコーダの基本的な考え方について説明しましょう。ここで説明するレンジコーダは「桁上がり」が発生するバージョンです。
算術符号は区間 [0, 1) を分割していきますが、レンジコーダは [0, 1) を分割するのではなく、最初に大きな区間たとえば [0, 1000) を設定して、それを小さな区間に分割していくことで符号化を行います。レンジコーダは整数で演算するので、記号列が長くなると当然ですが区間が狭くなって分割できなくなります。そのときは区間を引き伸ばすことで対応します。
たとえば、[0, 1000) を分割していくと [123, 124) になりました。もうこれ以上分割できませんね。そこで、区間をたとえば 100 倍して [12300, 12400) を分割することにします。このとき、区間全体の大きさは [0, 1000) ではなく、それを 1000 倍した [0, 100000) と考えるわけです。
単純に考えると、区間を表すために多倍長整数が必要になりますが、区間を引き伸ばすタイミングを定めることにより、通常の整数演算でレンジコーダをプログラムすることができます。また、区間全体の大きさも覚えておく必要はありません。レンジコーダは分割した区間の幅 (range) と下限値だけで符号化することができます。復号の処理でも、符号化と同じタイミングで range を引き伸ばしていくことで、符号語を記号列に復号することができます。レンジコーダは区間の下限値を符号語として出力します。
それでは具体的に説明しましょう。最初は区間の幅 range を 0x1000000 に設定し、下限値 low は 0 に初期化します。区間は [low, low + range) と表すことができるので、最初の区間は [0, 0x1000000) となります。また、range の初期値が 0x1000000 なので、low の値は 0 から 0xffffff までの範囲 (24 bit) になります。
記号の出現確率により区間を分割するところは算術符号と同じです。レンジコーダは range が一定の値より小さくなった時点で、range を引き伸ばすところがポイントです。レンジコーダでは、range が初期値の 1/256 (0x10000) より小さくなったら 256 倍します。これは処理をバイト単位で行うための工夫です。次の例を見てください。
[0x123456, 0x123456 + 0xabcd) = 256 倍 => [0x12345600, 0x12345600 + 0xabcd00)
いま low の値が 0x123456 で range の値が 0xabcd だとします。0xabcd < 0x10000 なので range を 256 倍します。このとき、low の値もいっしょに 256 倍すればいいわけです。これで区間を正しく表すことができますが、このままでは low の値がどんどん大きくなる一方ですね。そこで、low の値を一定の範囲内 (24 bit) に収めることを考えます。
range の値は 24 bit の範囲内に収まるので、low の計算は 24 bit の足し算になります。桁上がりの処理を工夫すれば、low を 24 bit で保持することが可能です。たとえば、次のように low の上位 8 bit (0x12) をバッファへ出力します。
[0x12345600, 0x12345600 + 0xabcd00) => [0x345600, 0x345600 + 0xabcd00) low (0x12345600) の上位 8 bit (0x12) をバッファへ出力 => (0x12)
値をバッファに溜めておけば、桁上がりには簡単に対応することができます。また、桁上がりが発生しないように工夫することができれば、上位 8 bit (0x12) をそのまま符号語としてファイルなどへ出力することができます。あとは、記号を読み込んで区間の分割と引き伸ばしを繰り返して、最後に low の値 (24 bit) を出力します。
簡単な例を示しましょう。記号列 "dcbbaaaa" を符号化します。記号の出現確率は次のようになります。
| a | b | c | d | |
|---|---|---|---|---|
| 出現確率 | 1/2 | 1/4 | 1/8 | 1/8 |
| 下限値 | 0 | 4/8 | 6/8 | 7/8 |
| 上限値 | 4/8 | 6/8 | 7/8 | 8/8 |
符号化の過程は次のようになります。
low = low + (range * 記号の下限値)
range = range * 記号の出現確率
[ low, range] [ low, range] (数値は 16 進数)
[ 0,1000000] - d -> [e00000, 200000]
[e00000, 200000] - c -> [f80000, 40000]
[f80000, 40000] - b -> [fa0000, 10000]
[fa0000, 10000] - b -> [fa8000, 4000] 256 倍して fa を出力
[800000, 400000] - a -> [800000, 200000]
[800000, 200000] - a -> [800000, 100000]
[800000, 100000] - a -> [800000, 80000]
[800000, 80000] - a -> [800000, 40000] low [80, 00, 00] を出力
符号語 => [fa, 80, 00, 00]
図：符号化の過程（レンジコーダ）
記号を読み込むたびに、range の値は小さくなり low の値は増えていきます。d, c, b, b まで記号を読み込むと、range は 0x4000 になり 0x10000 より小さくなります。ここで range と low を 256 倍して、low の上位 8 ビット (0xfa) を出力します。次に記号 a を読み込みます。range の値は小さくなりますが、a の下限値が 0 なので low の値は増えません。最後に low の値を出力して終了です。符号語は [0xfa, 0x80, 0, 0] になります。
次は復号について説明します。下限値 low と幅 range は符号化と同様に 0 と 0x1000000 に初期化します。符号語を code とすると、最初 low は 0 なので [0, range) の範囲で code に対応する記号を探すことになります。見つけた記号を c1 とすると、low と range の値は符号化と同様に次式で更新します。
low1 = low (0) + (range * 記号 c1 の下限値) range1 = range * 記号 c1 の出現確率
こんどは [low1, low1 + range1) の範囲で code に対応する記号を探します。ここで code から下限値の増分を引き算した値 code1 を求めてみます。すると、次の図に示すように code1 は区間 [0, range1) の符号語に対応していることがわかります。つまり、次は [0, range1) の範囲で code1 に対応する記号を探せばよいのです。
0 low1 code low1+range1 range
├──────┼────┼───────┼───────┤
┌──────┘ │ │
│ ┌──────┘ │
│ │ ┌──────┘
↓ ↓ ↓
├────┼───────┤
0 code1 range1
図：区間の更新
このように、符号語 code から下限値の増分を引き算することで、区間を [low, low + range) から [0, range) に変換することができるわけです。したがって、復号処理では下限値 low の値を覚えておく必要はありません。
range が 0x10000 より小さくなったら range を 256 倍するのは符号化と同じです。このとき符号語 code も 256 倍して、新しい符号語を 1 バイト読み込んで code に加算します。これで符号語を復号することができます。
それでは、復号の過程を具体的に説明しましょう。次の図を見てください。
code = code - (range * 記号の下限値)
range = range * 記号の出現確率
符号語を 3 バイト [fa, 80, 00] 読み込み code を初期化
[ code, range] [ code, range] (数値は 16 進数)
[fa8000, 1000000] - d -> [1a8000, 200000]
[1a8000, 200000] - c -> [ 28000, 40000]
[ 28000, 40000] - b -> [ 8000, 10000]
[ 8000, 10000] - b -> [ 0, 4000] 256 倍する
符号語を 1 バイト (00) code に加算
[ 0, 400000] - a -> [ 0, 200000]
[ 0, 200000] - a -> [ 0, 100000]
[ 0, 100000] - a -> [ 0, 80000]
[ 0, 80000] - a -> [ 0, 40000]
記号列 => "dcbbaaaa"
図：復号の過程（レンジコーダ）
最初に range と code を初期化します。code の範囲は 24 bit なので、3 バイト読み込んで 0xfa8000 に初期化します。次に、「記号の下限値 <= code / range < 記号の上限値」を満たす記号を探します。この場合、記号は d になります。そして、range を記号 d の出現確率で縮小して、code から (range * d の下限値) を引き算します。今度は 0x200000 の幅の中で 0x1a8000 に相当する記号を探すわけです。
d, c, b, b まで復号すると、range は 0x4000 になり 0x10000 より小さくなります。ここで range と code を 256 倍して、新しい符号語を 1 バイト読み込んで code に足し算します。この場合、符号語は 0 なので code の値は増えません。あとは、同じ処理を繰り返して記号列 "dcbbaaaa" を求めることができます。
それではプログラムを作りましょう。最初に記号の出現確率を求める関数を作ります。記号と記号列はシンボルとリストで表します。レンジコーダは整数で演算するので、出現確率は各記号の個数と記号の総数から求めます。記号と個数は連想リストに格納して変数 *count-table* に、記号の総数は変数 *total* にセットします。たとえば、記号列 (a b c c d d d d) は次のようになります。
*count-table* => ((a 1 7 8) (b 1 6 7) (c 2 4 6) (d 4 0 4)) *total* => 8
先頭要素が記号、2 番目の要素が記号の個数、3, 4 番目の要素が区間（下限値と上限値）を表します。記号の下限値と上限値は分数ではなく累積度数で表しています。*count-table* を作成する関数 make-count-table は次のようになります。
リスト：出現頻度表の作成
(defun make-count-table (buffer)
(setq *count-table* nil
*total* 0)
(let ((sum 0) work data)
; カウント
(dolist (code buffer)
(incf *total*)
(if (setq data (assoc code work))
(incf (cdr data))
(push (cons code 1) work)))
; *count-table* の作成
(dolist (code work)
(push (list (car code) ; 記号
(cdr code) ; 個数
sum ; 下限値
(+ sum (cdr code))) ; 上限値
*count-table*)
(incf sum (cdr code)))))
make-count-table は 算術符号 のプログラムとほとんど同じなので説明は不要でしょう。
次は記号列を符号化する関数 range-encode を作ります。次のリストを見てください。
リスト：レンジコーダの符号化
; 定数の定義
(defconstant *MAX-RANGE* #x1000000)
(defconstant *MIN-RANGE* #x10000)
(defconstant *MASK* #xffffff)
; 新しい range を計算する
(defun calc-range (range data)
(truncate (* range (second data)) *total*))
; 下限値を計算する
(defun calc-low (range data)
(truncate (* range (third data)) *total*))
; 上限値を計算する
(defun calc-high (range data)
(truncate (* range (fourth data)) *total*))
; 桁上がりの処理 (over は 0 or 1)
(defun overflow (code &optional (over 1))
(cond ((zerop over) code)
((null code) (list over))
((= 255 (car code))
(cons 0 (overflow (cdr code) 1)))
(t (cons (+ (car code) over) (cdr code)))))
; 符号化
(defun range-encode (buffer)
(let ((range *MAX-RANGE*) (low 0) data code)
(dolist (c buffer)
(setq data (assoc c *count-table*))
(incf low (calc-low range data))
(setq range (calc-range range data))
; 桁上がりのチェック
(when (<= *MAX-RANGE* low)
(setq code (overflow code))
(decf low *MAX-RANGE*))
; range の拡張
(while (< range *MIN-RANGE*)
; コードを出力
(push (ash low -16) code)
; 8 bit shift
(setq low (logand (ash low 8) *MASK*)
range (ash range 8))))
; low を出力
(dotimes (x 3 (reverse code))
(push (ash low -16) code)
(setq low (logand (ash low 8) *MASK*)))))
range-encode は buffer から記号をひとつずつ取り出して、区間の幅 range を狭めて下限値 low の値を計算します。range の初期値は *MAX-RANGE* (#x1000000) で、*MIN-RANGE* (#x10000) より小さくなったならば range と low の値を 256 倍するとともに、low の上位 8 ビットを符号語として変数 code のリストに格納します。
range と low の計算は関数 calc-range と calc-low で行います。range は記号の出現確率で縮小すればいいので (range * 記号の総数) / *total* を calc-range で計算します。low の増分は区間の下限値なので (range * 記号の下限値) / *total* を calc-low で計算します。Common Lisp の関数 / は割り切れないと結果を分数で返すので、小数点以下を切り捨てるために関数 truncate を使っていることに注意してください。
low の値が *MAX-RANGE* 以上になったならば、桁上がりの処理を関数 overflow で行います。符号語は code に逆順で格納されているので、最初の要素に 1 を足して値が 256 になったならば、その要素を 0 にして次の要素に 1 を足し算します。overflow はこの処理を再帰呼び出しで実現しています。難しい処理ではないので、詳細はプログラムリストお読みください。
幅 range が *MIN-RANGE* より小さくなったならば、range と low を 256 倍します。最初に low の上位 8 bit を code に格納します。low の値は 24 bit の範囲内なので、関数 ash で low を 16 bit 右シフトすれば上位 8 bit の値を取り出すことができます。次に、range と low を ash で 8 bit 左シフトします。これで値は 256 倍されます。このままでは low の値が 24 bit の範囲に収まらないので、論理積を求める logand で符号語として出力した部分を取り除きます。
記号列を最後まで読み込んだら、low の値を code に出力します。low の上位 8 bit から順番に出力していることに注意してください。これでプログラムは完成です。
簡単な実行例を示します。
(make-count-table '(a b c c d d d d)) => nil (range-encode '(a b c c d d d d)) => (250 128 0 0) (make-count-table '(a b b b b b b b b b)) => nil (range-encode '(a b b b b b b b b b)) => (230 102 102)
次は復号のプログラムを作りましょう。次のプログラムを見てください。
リスト：レンジコーダの復号
; 記号を見つける
(defun find-range-code (code range)
(find-if #'(lambda (data)
(and (<= (calc-low range data) code)
(< code (calc-high range data))))
*count-table*))
; 復号
(defun range-decode (n code-list)
(let ((range *MAX-RANGE*) (code 0) data buffer)
; code の初期化
(dotimes (x 3)
(setq code (+ (ash code 8) (pop code-list))))
; 復号
(dotimes (x n (reverse buffer))
(setq data (find-range-code code range))
; コードを出力
(push (first data) buffer)
(decf code (calc-low range data))
(setq range (calc-range range data))
; range の拡張
(while (< range *MIN-RANGE*)
(setq range (ash range 8))
(setq code (+ (ash code 8) (pop code-list)))))))
関数 find-range-code は、記号の出現確率 *count-table* から符号語 code に対応する記号を探します。記号の下限値と上限値（累積度数）を calc-low と calc-high で range の幅に変換し、code がその範囲内にある記号を find-if で探します。
関数 range-decode の引数 n が復号する記号の総数で、code-list が入力データ（符号語のリスト）を表します。range は *MAX-RANGE* に初期化します。符号語 code の範囲は 24 bit なので、code-list から 3 バイト分読み込んで初期化します。あとは、find-range-code で記号を求めて buffer にセットし、range と code の値を更新します。range が *MIN-RANGE* よりも小さくなったならば、range と code の値を 256 倍します。このとき、code-list から符号語を 1 バイト分読み込んで code に加算します。あとはこれを n 回繰り返すだけです。
それでは実行してみましょう。
(make-count-table '(a b c c d d d d)) => nil (range-encode '(a b c c d d d d)) => (250 128 0 0) (range-decode 8 '(250 128 0 0)) => (a b c c d d d d) (make-count-table '(a b b b b b b b b b)) => nil (range-encode '(a b b b b b b b b b)) => (230 102 102) (range-decode 10 '(230 102 102)) => (a b b b b b b b b b)
きちんと復号できましたね。
ところで、適応型レンジコーダも簡単に作成することができます。興味のある方は挑戦してみてください。
Common Lisp 入門 の番外編です。今回は二分木の操作の中でちょっと複雑な「データの削除」を取り上げます。
二分木からデータを削除する場合、データの挿入と違ってちょっと面倒です。削除するデータが「葉」の場合は、それを削除するだけでよいのですが、データが木の途中にある場合は、二分木の条件を満たすように削除しないといけません。まず最初に「葉」を削除する場合を説明しましょう。次の図を見てください。
14 14
／ ＼ ／ ＼
／ ＼ ／ ＼
12 16 => 12 16
／ ＼ ／ ＼ ／ ＼ ／ ＼
11 13 15 17 11 13 NIL 17
↑
15 を削除する 削除
図：データの削除（葉の場合）
上図の 二分木でデータ 15 を削除します。15 は「葉」にあたるので、それを削除するだけでよいですね。親節 16 の左側の子 (left) に終端 (NIL) を代入するだけです。
次に、子がひとつある場合を考えてみましょう。次の図を見てください。
14 14
／ ＼ ／ ＼
／ ＼ ／ ＼
12 16 => 12 15
／ ＼ ／ ／ ＼
11 13 15 11 13
16 を削除する
図：データの削除（子がひとつの場合）
節 16 は子 15 がひとつあります。この節を削除する場合、その子である 15 と置き換えれば 二分木の条件は満たされます。これも簡単ですね。
問題は、子がふたつある節を削除する場合です。この場合、削除するデータの右部分木の中から最小値のデータ [*1] を探して、それと削除するデータと置き換えれば、二分木の条件を満たすことができます。次の図を見てください。
14 15 <- 最小値と置き換え
／ ＼ ／ ＼
／ ＼ ／ ＼
12 16 => 12 16
／ ＼ ／ ＼ ／ ＼ ／ ＼
11 13 15 17 11 13 NIL 17
↑
14 を削除する 削除
図：データの削除（子がふたつの場合）
節 14 はふたつの子を持っています。節 14 を削除する場合、右部分木の中の最小値 15 と 14 を置き換えて、15 を格納していた節を削除します。節が最小値を格納している場合、その節の左側の子は存在しません。つまり、削除する節は上図のように「葉」の場合か、子はあっても右側の子ひとつしかないので、その節を削除することは簡単です。
それでは、二分木からデータを削除するプログラムを Common Lisp で作りましょう。次のプログラムを見てください。
リスト：2 分木の操作
; 節の定義
(defstruct Node
data left right)
; データの挿入
(defun insert-tree (node num)
(cond ((null node)
(setq node (make-Node :data num)))
((< num (Node-data node))
(setf (Node-left node)
(insert-tree (Node-left node) num)))
((< (Node-data node) num)
(setf (Node-right node)
(insert-tree (Node-right node) num))))
; 最後は node を返す
node)
; 表示
(defun print-tree (n node)
(when node
(print-tree (1+ n) (Node-left node))
(dotimes (x n) (princ " "))
(format t "~D~%" (Node-data node))
(print-tree (1+ n) (Node-right node))))
節は構造体で定義します。格納するデータは整数値とします。データの挿入は関数 insert-tree で、2 分木の表示は print-tree で行います。どちらの関数も簡単ですね。2 分木の操作は拙作のページ xyzzy Lisp Programming の Common Lisp 入門 二分探索木 で詳しく説明しています。よろしければ参考にしてください。
insert-tree と print-tree の実行例を示します。
(setq *root* nil)
(dolist (x '(40 20 10 30 50))
(setq *root* (insert-tree *root* x)))
(print-tree 0 *root*)
10
20
30
40
50
ルートの節は変数 *root* に格納します。*root* は nil に初期化しておくことをお忘れなく。print-tree は昇順にデータを出力するので、左部分木が上に右部分木が下に表示されます。
次に、木の中から最小値を探して、それと削除するデータと置き換える関数を作成します。次のリストを見てください。
リスト：最小値を探して置き換える
(defun delete-min-tree (node replace-node)
(cond ((null (Node-left node))
(setf (Node-data replace-node) (Node-data node)
node (Node-right node)))
(t
(setf (Node-left node)
(delete-min-tree (Node-left node) replace-node))))
; node を返す
node)
関数 delete-min-tree の引数 node が探索する部分木、replace-node が削除するデータを格納している節です。二分木の場合、最小値を探すことはとても簡単です。左の子を順番にたどっていき、左の子がない節に行き着いたとき、その節のデータが最小値になります。
最初に、node の左の子がないかチェックします。そうであれば、この節のデータが最小値です。replace-node の data を node のデータに書き換えて、node の右の子を返します。これで、親節の左部分木が書き換えられ、最小値を持つ節が削除されます。node が葉の場合、node の右の子は nil なので、単純に削除されることになります。
node の左の子があれば、delete-min-tree を再帰呼び出しします。その返り値を node の左の子にセットして、node をそのまま返します。なお、この関数は削除関数のサブルーチンとして使うことを前提にしているので、node に空の木 (nil) を与えると動作しません。ご注意くださいませ。
それでは、データを削除する関数 delete-tree を作りましょう。次のリストを見てください。
リスト：データの削除
(defun delete-tree (node num)
(when node
(cond ((< num (Node-data node))
(setf (Node-left node)
(delete-tree (Node-left node) num)))
((< (Node-data node) num)
(setf (Node-right node)
(delete-tree (Node-right node) num)))
(t ; データ発見
(cond ((null (Node-left node))
(setf node (Node-right node)))
((null (Node-right node))
(setf node (Node-left node)))
(t
(setf (Node-right node)
(delete-min-tree (Node-right node) node)))))))
; node を返す
node)
引数 node が節で、num が削除するデータです。node が nil ならば空の木なので、何もしないで node を返します。削除するデータが見つからない場合や、空の木を与えた場合がこれに相当します。次に、二分木の中から num を探索します。num が node の data より小さければ左部分木を、大きければ右部分木をたどります。
データが見つかった場合、node の子の有無をチェックします。まず、左の子がない場合は右の子を返します。ここで、node が「葉」の場合は、右の子が nil なので単純に削除されることになります。次に、左の子があって右の子がない場合は左の子を返します。
最後の t 節が、子が 2 つある場合の削除です。delete-min-tree を呼び出して右部分木の中から最小値を探して、データと最小値を置き換えます。delete-min-tree の返り値は右の子にセットします。最後に node を返します。
それでは、簡単な実行例を示しましょう。
(setq *root* nil)
(dolist (x '(50 30 20 10 40 60 70))
(setq *root* (insert-tree *root* x)))
(print-tree 0 *root*)
10
20
30
40
50
60
70
; 葉 (70) の削除
(setq *root* (delete-tree *root* 70))
(print-tree 0 *root*)
10
20
30
40
50
60
; 子をひとつ持つ節 (20) の削除
(setq *root* (delete-tree *root* 20))
(print-tree 0 *root*)
10
30
40
50
60
; 子をふたつ持つ節 (30) の削除
(setq *root* (delete-tree *root* 30))
(print-tree 0 *root*)
10
40
50
60
正常に動作していますね。