今回はデータ圧縮法のひとつである「算術符号」を取り上げます。「ハフマン符号」は各記号 (文字) の出現確率を調べ、頻繁に現れる記号は短いビットで表し、あまり現れない記号は長いビットで表すことで、データを圧縮する古典的な方法です。古典的とはいっても、ほかのアルゴリズムと簡単に組み合わせることができるため、ハフマン符号は今でも現役のアルゴリズムです。
記号の出現確率だけを利用してデータを圧縮する方法は、ハフマン符号のほかには「算術符号」という方法があります。一般に、算術符号はハフマン符号よりも性能が良いのですが、実現方法が難しくて実行速度がハフマン符号化よりも遅く、なおかつ特許の問題もあって、実際に使われることはあまりありませんでした。
近年になると「レンジコーダ (RangeCoder)」という方法が考案され、注目を集めるようになりました。レンジコーダは原理的には算術符号と同じ方法ですが、参考文献 3 『データ圧縮の基礎から応用まで』によると 『(おそらく) 特許フリー』 とのことで、性能は算術符号に比べるとわずかに劣化しますが、実現方法はとても簡単で実行速度も高速です。もちろん、ハフマン符号よりも高性能です。
レンジコーダのプログラムは、算術符号を理解していると簡単に作ることができます。そして、算術符号は Common Lisp の分数を使うと、算術符号の原理そのままにプログラムを作ることができます。そこで、今回は算術符号のプログラムを Common Lisp で作ってみましょう。これから作成するプログラムは学習用なので実用性はまったくありませんが、実際にプログラムを作ることで算術符号の理解は深まると思います。
最初に「算術符号」の基本的な考え方について説明します。算術符号は記号 (文字) 列全体をひとつの符号語にする方法で、1960 年代に P. Elias によって提案されました。
算術符号は、記号列を実数 0 と 1 の間の区間を用いて表します。たとえば、記号は {a, b, c} の 3 種類があり、出現確率はそれぞれ 0.2, 0.6, 0.2 とします。算術符号は、区間を記号の出現確率に比例した小区間に分割していくことで符号化を行います。それでは、記号列 abbbc を符号化してみましょう。次の図を見てください。
1.0 ┬ ┌→ 0.2 ┬ ┌→ 0.16 ┬ ┌→ 0.136 ┬ ┌→ 0.1216 ┬ ─→ 0.1216 ┬
c │ │ │ │ │ │ │ │ │ │
0.8 ┼ │ 0.16 ┼ ┘ 0.136 ┼ ┘ 0.1216 ┼ ┘ 0.11296 ┼ ┐ 0.1198728 ┼
│ │ │ │ │ │ │ │
b │ │ │ │ │ │ │ │
│ │ │ │ │ │ │ │
0.2 ┼ ┘ 0.04 ┼ ┐ 0.064 ┼ ┐ 0.0784 ┼ ┐ 0.08704 ┼ │ 0.114688 ┼
a │ │ │ │ │ │ │ │ │ │
0.0 ┴ ─→ 0.0 ┴ └→ 0.04 ┴ └→ 0.64 ┴ └→ 0.0784 ┴ └→ 0.11296 ┴
記号 a ab abb abbb abbbc
図 : 符号化の過程 (算術符号)
x 以上 y 未満の区間を [x, y) と表すことにします。区間の初期値は [0, 1) です。記号を読み込んだら区間 [0, 1) を分割します。記号が a ならば区間の 0 から 0.2 までの部分、b ならば 0.2 から 0.8 までの部分、c ならば 0.8 から 1.0 までの部分に分割します。
最初の記号は a なので区間は [0, 0.2) となります。次の記号は b なので、区間 [0, 0.2) の 0.2 から 0.8 までの部分 [0.04, 0.16) が新しい区間となります。このように、記号を読み込むたびに区間を分割していくと、記号列 abbbc を表す区間は [0.11296, 0.1216) となります。そして、実際の算術符号の符号語は、この区間に含まれるひとつの実数を指定します。
ここで符号語を 2 進数で表して、区間内で小数点以下のビット数の少ない値を選ぶことにしましょう。たとえば、0.1171875 を 2 進数で表すと次のようになります。
0.1171875 = 1/16 + 1/32 + 1/64 + 1/128 = (0.0001111)2
0001111 の 7 bit を符号語として出力すれば、記号列 abbbc の 5 文字を 7 bit に圧縮することができるわけです。この例では 1 文字が 1.4 bit に圧縮されましたが、記号の出現確率によっては 1 文字が 1 bit 未満に圧縮できる場合もあります。ちなみに、ハフマン符号では 1 文字が 1 bit 未満になることはありえません。
次は復号を説明します。ここでは説明の都合上、符号語は下限値の 0.11296 とします。0.11296 は [0, 0.2) の間にあるので、最初の記号は a であることがわかります。次に、a を表す区間 [0, 0.2) を [0, 1.0) になるように拡大すると、符号語は次のように変換できます。
新しい符号語 = (符号語 - 記号の下限値) / 記号の区間幅
= (0.11296 - 0) / 0.2
= 0.5648
新しい符号語 0.5648 は [0.2, 0.8) の間にあるので、次の記号は b であることがわかります。このような操作を繰り返し行うことで、次のように記号列 abbbc を求めることができます。
| 符号語 | 記号 | 区間 |
|---|---|---|
| 0.11296 | a | [0, 0.2) |
| 0.5648 | b | [0.2, 0.8) |
| 0.608 | b | [0.2, 0.8) |
| 0.68 | b | [0.2, 0.8) |
| 0.8 | c | [0.8, 1) |
| 0 |
ところで、算術符号には 2 つの問題点があります。ひとつは記号列の最後を判定できないことです。さきほどの復号の例では最後に符号語が 0 になりましたが、このあとも記号 a を復号することができます。つまり、符号語 0.11296 は記号列 abbbc だけではなく、abbbca, abbbcaa, abbcaaa... などにも復号することができるのです。
この問題は、終端を表す記号を用意して終端記号を復号したら終了する、または、記号の総数をファイルの先頭に書き込んでおく、などといった方法で解決することができます。
もうひとつは、入力する記号列が長くなるほど、より多くの桁数が必要になることです。また、浮動小数点演算の誤差も考慮しなければいけません。これはとても大きな問題点で、解決するまでに 10 年以上の時間がかかりました。
1981 年に C. B. Jones によって発表された算術符号 (Jones 符号) は、実数のかわりに整数で演算するように工夫されています。Jones 符号に興味のある方は参考文献 1 『文書データ圧縮アルゴリズム入門』をお読みください。参考文献 2 『Ｃ言語による最新アルゴリズム事典』にも算術圧縮のプログラム (Ｃ言語) があります。
それでは、符号化のプログラムを作りましょう。最初に記号の出現確率を求めます。記号と記号列はシンボルとリストで表します。記号の出現確率は連想リストに格納して返します。
((記号 出現確率 区間の下限 区間の上限) ... )
先頭要素が記号、2 番目の要素が出現確率、3, 4 番目の要素が区間を表します。出現頻度表を作成する関数 make-frequency は次のようになります。
リスト : 出現頻度表の作成
;;; 記号
(defconstant syms '(a b c d e f g h))
;;; 記号のカウント
(defun count-symbol (xs)
(let ((table (mapcar (lambda (x) (cons x 0)) syms)))
(dolist (x xs table)
(let ((code (assoc x table)))
(unless code
(error "invalid symbol ~a~%" x))
(incf (cdr code))))))
;;; 記号の出現確率と区間を求める
(defun make-cumulative (xs size &optional (a 0))
(if (null xs)
nil
(cons (list (caar xs) ; 記号
(/ (cdar xs) size) ; 出現確率
(/ a size) ; 下限
(/ (+ a (cdar xs)) size)) ; 上限
(make-cumulative (cdr xs) size (+ a (cdar xs))))))
;;; 出現頻度表の作成
(defun make-frequency (xs)
(make-cumulative (count-symbol xs) (length xs)))
使用する記号はリスト SYMS に登録しておきます。最初に記号をカウントします。この処理を関数 count-symbol で行います。返り値は連想リストです。次に、関数 make-cumulative で記号の出現確率と区間を求めます。引数 XS には count-symbol の返り値を渡します。SIZE は記号の総数です。引数 A は記号の累積度数を表します。あとは、出現確率と区間の下限と上限を計算し、それを連想リストに格納して返すだけです。
次は記号列を符号化する関数 encode を作ります。Common Lisp の分数を使うと、encode はとても簡単にプログラムできます。次のリストを見てください。
リスト : 符号化
(defun encode (xs)
(let ((table (make-frequency xs))
(low 0)
(high 1))
(dolist (x xs (values table (make-code low high)))
(let ((c (assoc x table))
(w (- high low)))
(setq low (+ (* w (third c)) low)
high (+ (* w (second c)) low))))))
encode の引数 XS には記号列を渡します。最初に make-frequency で出現頻度表を作成して変数 TABLE にセットし、変数 LOW と HIGH を 0 と 1 に初期化します。次に、dolist で XS から記号を一つずつ取り出して変数 X にセットし、assoc で TABLE から X を検索します。
あとは算術符号の原理そのままで、区間 [LOW, HIGH) を記号の出現確率で分割していくだけです。分数を使っているので桁数や計算誤差を気にする必要はありません。最後に values で出現頻度表と符号語を返します。関数 make-code は符号語をビット列 (要素が 0 と 1 だけのリスト) で出力します。make-code は次のリストを見てください。
リスト : 符号をビット列で出力
(defun make-code (low high &optional (value 0) (n 1/2))
(cond
((and (<= low value) (< value high)) nil)
((<= high (+ value n))
(cons 0 (make-code low high value (/ n 2))))
(t
(cons 1 (make-code low high (+ value n) (/ n 2))))))
1/2, 1/4, 1/8, 1/16, ... (1/2)k の値を使って、LOW 以上 HIGH 未満の値 VALUE を求めます。VALUE に (1/2)k を加えた場合は返り値のリストに 1 を追加し、そうでない場合は 0 を追加します。これで符号語をビット列で表すことができます。
それでは実際に試してみましょう。
* (encode '(a b b b c)) ((A 1/5 0 1/5) (B 3/5 1/5 4/5) (C 1/5 4/5 1) (D 0 1 1) (E 0 1 1) (F 0 1 1) (G 0 1 1) (H 0 1 1)) (0 0 0 1 1 1 1)
low = 0, high = 1
code = A, low = 0, high = 1/5
code = B, low = 1/25, high = 4/25
code = B, low = 8/125, high = 17/125
code = B, low = 49/625, high = 76/625
code = C, low = 353/3125, high = 76/625
記号列 abbbc を 7 bit で表すことができました。符号化の様子を表示すると、区間を分割していく様子がよくわかると思います。
もうひとつ例を示しましょう。
* (encode '(a a a a a a a a a b)) ((A 9/10 0 9/10) (B 1/10 9/10 1) (C 0 1 1) (D 0 1 1) (E 0 1 1) (F 0 1 1) (G 0 1 1) (H 0 1 1)) (0 1 1)
low = 0, high = 1
code = A, low = 0, high = 9/10
code = A, low = 0, high = 81/100
code = A, low = 0, high = 729/1000
code = A, low = 0, high = 6561/10000
code = A, low = 0, high = 59049/100000
code = A, low = 0, high = 531441/1000000
code = A, low = 0, high = 4782969/10000000
code = A, low = 0, high = 43046721/100000000
code = A, low = 0, high = 387420489/1000000000
code = B, low = 3486784401/10000000000, high = 387420489/1000000000
記号列 aaaaaaaaab を 011 の 3 bit に圧縮することができました。1 文字あたり 0.3 bit になります。これがハフマン符号にはない算術符号の特徴です。
次は、復号のプログラムを作ってみましょう。最初に、符号語 (ビット列) を数値 (分数) に変換するプログラムを作ります。次のリストを見てください。
リスト : 符号語を数値に変換
(defun code-to-number (xs)
(loop for x in xs
for n = 1/2 then (/ n 2)
when (plusp x) sum n))
要素が k 個のビット列は (1/2 1/4 1/8 1/16 ... (1/2)k) に対応しています。この値を変数 N で表しています。あとはビット列の要素 X が正ならば、loop マクロの sum で N を加算していくだけです。たとえば、記号列 abbbc の符号語 (0 0 0 1 1 1 1) は 15/128 になります。
数値から記号を求める処理は簡単です。出現頻度表から「第 3 要素 (下限値) <= 数値 < 第 4 要素 (上限値)」を満たす記号を関数 find-if で探すだけです。たとえば、記号列 abbbc の出現頻度表の値は次のようになっています。
((A 1/5 0 1/5) (B 3/5 1/5 4/5) (C 1/5 4/5 1) ... )
15/128 は 0 <= 15/128 < 1/5 を満たすので記号は A となります。あとは、A を表す区間を [0, 1) になるように拡大すると、符号語は次の式で変換されます。
新しい符号語 = (符号語 - 記号の下限値) / 記号の区間幅
そして、新しい符号語に対応する記号を探します。あとはこの処理を繰り返すだけです。プログラムは次のようになります。
リスト : 復号
(defun decode (n code table)
(let ((value (code-to-number code)) buffer)
(dotimes (i n (reverse buffer))
(let ((data (find-if (lambda (xs)
(and (<= (third xs) value) (< value (fourth xs))))
table)))
(push (first data) buffer)
(setq value (/ (- value (third data)) (second data)))))))
;;; 復号のテスト
(defun test (xs)
(multiple-value-bind
(table code)
(encode xs)
(decode (length xs) code table)))
関数 decode の引数 N が復号する記号数、CODE が符号語 (ビット列)、TABLE が出現頻度表です。最初に code-to-number で符号語を数値に変換します。find-if で求めた記号は BUFFER に追加し、それから数値 VALUE を変換します。あとはこれを N 回繰り返すだけです。
それでは実行してみましょう。
* (test '(a b b b c)) (A B B B C)
value: 15/128 => symbol: A value: 75/128 => symbol: B value: 247/384 => symbol: B value: 851/1152 => symbol: B value: 3103/3456 => symbol: C
きちんと復号されましたね。もうひとつ例を示しましょう。
* (test '(a a a a a a a a a b)) (A A A A A A A A A B)
value: 3/8 => symbol: A value: 5/12 => symbol: A value: 25/54 => symbol: A value: 125/243 => symbol: A value: 1250/2187 => symbol: A value: 12500/19683 => symbol: A value: 125000/177147 => symbol: A value: 1250000/1594323 => symbol: A value: 12500000/14348907 => symbol: A value: 125000000/129140163 => symbol: B
今回のプログラムでは復号で記号の総数が必要になりますが、終端を表す記号を用意して終端記号を復号したら終了する方法もあります。プログラムは簡単に改造できるので、興味のある方は試してみてください。
今までに説明した算術符号は「静的符号化」といい、あらかじめ記号の出現確率を調べておいて、それに基づいて入力記号列を符号化していく方法です。この方法では、ハフマン符号がもっとも有名でしょう。
これに対し「動的符号化」は、入力記号列の符号化を行いながら記号の出現確率を変化させる方法で、「適応型符号化」とも呼ばれています。最初は、どの記号も同じ確率で出現すると仮定して、記号列を読み込みながら記号の出現確率を修正し、その時点での出現確率に基づいて記号の符号化を行います。
動的符号化の特徴は、入力記号列の性質 (出現確率) の変化に適応できることですが、このほかにも長所があります。静的符号化の場合、復号するときに符号化で用いた記号の出現確率が必要になります。
このため、実際にファイルを圧縮するプログラムでは、記号の出現頻度表を出力ファイルの先頭に付加する方法が用いられます。ところが、動的符号化の場合、記号の出現確率は復号しながら求めることができるので、出現頻度表を付加する必要はありません。
また、静的符号化でファイルを圧縮する場合、記号の出現頻度を求めるときにファイルからデータを読み込み、符号化を行うときに再度ファイルからデータを読み込む必要があります。このようにデータの入力が 2 回必要な圧縮アルゴリズムを「2 パスの圧縮アルゴリズム」といいます。動的符号化は 1 パスで済むので、オンラインでのデータ圧縮にも対応することができます。
このように、動的符号化には有利な点があるため、ハフマン符号を動的符号化に対応させた「適応型ハフマン符号」が考案されています。しかしながら、適応型ハフマン符号は実装方法が難しく、処理速度も遅いという欠点があります。
これに対し、適応型算術符号 (レンジコーダ) は簡単な方法で実装することができ、適応型レンジコーダは処理速度も遅くありません。優れた実装方法なのです。
そこで、今回は適応型算術符号のプログラムを Common Lisp で作ります。Common Lisp の分数を使うと、算術符号の原理そのままにプログラムを作ることができます。そして、適応型算術符号にも簡単に対応することができます。
これから作成するプログラムは学習用なので実用性はありませんが、実際にプログラムを作ることで適応型算術符号の理解は深まると思います。
それでは、符号化のプログラムから作りましょう。最初に記号の出現頻度表を作成します。適応型算術符号の場合、最初はどの記号も同じ確率で出現すると仮定するのが一般的なので、出現する可能性がある記号はすべて出現頻度を 1 に初期化します。今回は学習用のプログラムなので、記号は A, B, C, D, E, F, G, H の 8 種類に限定します。プログラムは次のようになります。
リスト : 出現頻度表の作成
(defun make-frequency-ad (xs &optional (sum 0))
(if (null xs)
nil
(cons (list (car xs) 1 sum (1+ sum))
(make-frequency-ad (cdr xs) (1+ sum)))))
出現頻度表は連想リストを使って表していて、要素はリスト (記号 出現頻度 区間の下限 区間の上限) です。区間の上限と下限は累積度数で表します。したがって、出現頻度表を更新するときは、記号の出現頻度を増やすだけではなく、それ以降の記号の累積度数も増やす必要があります。出現頻度表は次のように初期化されます。
((A 1 0 1) (B 1 1 2) (C 1 2 3) (D 1 3 4) (E 1 4 5) (F 1 5 6) (G 1 6 7) (H 1 7 8))
次は記号列を符号化する関数 encode-ad を作ります。次のリストを見てください。
リスト : 符号化
(defun encode-ad (xs)
(let ((table (make-frequency-ad syms))
(sum (length syms))
(low 0)
(high 1))
(dolist (x xs (make-code low high))
(let ((c (find x table :key #'car))
(w (- high low)))
(setq low (+ (* w (/ (third c) sum)) low)
high (+ (* w (/ (second c) sum)) low))
;; 出現頻度表の更新
(incf sum)
(update-frequency x table)))))
関数 encode-ad の引数 XS には記号列を渡します。TABLE から記号 X を検索するため、関数 find のキーワード :key に car を指定しています。あとのプログラムは算術符号の原理そのままで、区間 [LOW, HIGH) を記号の出現確率で分割していくだけです。適応型算術符号の場合、出現確率は記号を読み込むたびに変化するので、その時点での出現確率を記号のデータと記号の総数 SUM で計算していることに注意してください。そして、関数 update-frequency で出現頻度表を更新します。
リスト : 出現頻度表の更新
(defun update-frequency (code freq)
(let ((xs (member code freq :key #'car)))
;; 読み込んだ記号を更新
(incf (second (car xs)))
(incf (fourth (car xs)))
;; それ以降の記号の累積度数を更新
(dolist (x (cdr xs))
(incf (third x))
(incf (fourth x)))))
最初に、関数 member で出現頻度表 FREQ から記号 CODE を探索します。CODE の出現頻度とその区間の上限値を incf で増やします。あとは dolist でそれ以降の記号の区間の下限値と上限値を incf で増やすだけです。
これでプログラムは完成です。それでは実行してみましょう。
* (encode-ad '(a b b b c)) (0 0 0 0 1 0 0 0 0 0 0 1)
code = A, low = 0, high = 1/8, sum = 9 ((A 2 0 2) (B 1 2 3) (C 1 3 4) (D 1 4 5) (E 1 5 6) (F 1 6 7) (G 1 7 8) (H 1 8 9)) code = B, low = 1/36, high = 1/24, sum = 10 ((A 2 0 2) (B 2 2 4) (C 1 4 5) (D 1 5 6) (E 1 6 7) (F 1 7 8) (G 1 8 9) (H 1 9 10)) code = B, low = 11/360, high = 1/30, sum = 11 ((A 2 0 2) (B 3 2 5) (C 1 5 6) (D 1 6 7) (E 1 7 8) (F 1 8 9) (G 1 9 10) (H 1 10 11)) code = B, low = 41/1320, high = 7/220, sum = 12 ((A 2 0 2) (B 4 2 6) (C 1 6 7) (D 1 7 8) (E 1 8 9) (F 1 9 10) (G 1 10 11) (H 1 11 12)) code = C, low = 83/2640, high = 499/15840, sum = 13 ((A 2 0 2) (B 4 2 6) (C 2 6 8) (D 1 8 9) (E 1 9 10) (F 1 10 11) (G 1 11 12) (H 1 12 13))
記号列が (a b b b c) の場合、静的符号化では 7 bit に符号化できましたが、適応型算術符号では 12 bit になりました。このように、出現する記号の種類が少なく記号列の長さが短いデータでは、適応型算術符号の圧縮率が低下するのはしょうがないでしょう。ですが、記号列が長くて記号の出現確率の違いが大きくなると、適応型算術符号でも圧縮率は向上するので大丈夫です。
次は復号のプログラムを作ります。次のリストを見てください。
リスト : 適応型算術符号の復号
(defun decode-ad (n code)
(let ((table (make-frequency-ad syms))
(sum (length syms))
(value (code-to-number code)) buffer)
(dotimes (i n (reverse buffer))
(let ((data (find-if (lambda (xs)
(and (<= (/ (third xs) sum) value)
(< value (/ (fourth xs) sum))))
table)))
(push (car data) buffer)
(setq value (/ (- value (/ (third data) sum))
(/ (second data) sum)))
(incf sum)
(update-frequency (car data) table)))))
関数 decode の引数 N が復号する記号の総数、CODE が符号語 (ビット列) を表します。最初に make-frequency-ad で出現頻度表を作成し、code-to-number で符号語を数値に変換します。記号の復号処理は TABLE から「第 3 要素 (下限値)<= 数値 < 第 4 要素 (上限値)」を満たす記号を find-if で探すだけです。上限値と下限値は累積度数で表されているので、記号の総数 SUM で割り算していることに注意してください。
find-if で求めた記号は BUFFER に追加し、それから数値 VALUE を変換します。適応型算術符号の場合、出現確率は記号を読み込むたびに変化するので、その時点での出現確率を記号のデータと SUM で計算していることに注意してください。そして、関数 update-frequency で出現頻度表を更新します。あとはこれを N 回繰り返すだけです。
それでは実行してみましょう。
* (decode-ad 5 '(0 0 0 0 1 0 0 0 0 0 0 1)) (A B B B C)
value = 129/4096 => symbol: A, sum = 9 ((A 2 0 2) (B 1 2 3) (C 1 3 4) (D 1 4 5) (E 1 5 6) (F 1 6 7) (G 1 7 8) (H 1 8 9)) value = 129/512 => symbol: B, sum = 10 ((A 2 0 2) (B 2 2 4) (C 1 4 5) (D 1 5 6) (E 1 6 7) (F 1 7 8) (G 1 8 9) (H 1 9 10)) value = 137/512 => symbol: B, sum = 11 ((A 2 0 2) (B 3 2 5) (C 1 5 6) (D 1 6 7) (E 1 7 8) (F 1 8 9) (G 1 9 10) (H 1 10 11)) value = 173/512 => symbol: B, sum = 12 ((A 2 0 2) (B 4 2 6) (C 1 6 7) (D 1 7 8) (E 1 8 9) (F 1 9 10) (G 1 10 11) (H 1 11 12)) value = 293/512 => symbol: C, sum = 13 ((A 2 0 2) (B 4 2 6) (C 2 6 8) (D 1 8 9) (E 1 9 10) (F 1 10 11) (G 1 11 12) (H 1 12 13))
きちんと復号されましたね。このように、算術符号は出現頻度表を更新するだけで簡単に適応型符号を実現することができます。
;;;
;;; ari.lisp : 算術符号
;;;
;;; Copyright (C) 2003-2020 Makoto Hiroi
;;;
;;; 記号
(defconstant syms '(a b c d e f g h))
;;; 静的符号化
;;; 記号のカウント
(defun count-symbol (xs)
(let ((table (mapcar (lambda (x) (cons x 0)) syms)))
(dolist (x xs table)
(let ((code (assoc x table)))
(unless code
(error "invalid symbol ~a~%" x))
(incf (cdr code))))))
;;; 累積度数表
(defun make-cumulative (xs size &optional (a 0))
(if (null xs)
nil
(cons (list (caar xs) ; 記号
(/ (cdar xs) size) ; 出現確率
(/ a size) ; 下限
(/ (+ a (cdar xs)) size)) ; 上限
(make-cumulative (cdr xs) size (+ a (cdar xs))))))
;;; 出現頻度表の作成
(defun make-frequency (xs)
(make-cumulative (count-symbol xs) (length xs)))
;;; 符号語の生成
(defun make-code (low high &optional (value 0) (n 1/2))
(cond
((and (<= low value) (< value high)) nil)
((<= high (+ value n))
(cons 0 (make-code low high value (/ n 2))))
(t
(cons 1 (make-code low high (+ value n) (/ n 2))))))
;;; 符号化
(defun encode (xs)
(let ((table (make-frequency xs))
(low 0)
(high 1))
(dolist (x xs (values table (make-code low high)))
(let ((c (assoc x table))
(w (- high low)))
(setq low (+ (* w (third c)) low)
high (+ (* w (second c)) low))))))
;;; 符号を数値に変換
(defun code-to-number (xs)
(loop for x in xs
for n = 1/2 then (/ n 2)
when (plusp x) sum n))
;;; 復号
(defun decode (n code table)
(let ((value (code-to-number code)) buffer)
(dotimes (i n (reverse buffer))
(let ((data (find-if (lambda (xs) (and (<= (third xs) value) (< value (fourth xs)))) table)))
(push (first data) buffer)
(setq value (/ (- value (third data)) (second data)))))))
;;; 復号のテスト
(defun test (xs)
(multiple-value-bind
(table code)
(encode xs)
(decode (length xs) code table)))
;;; 適応型符号化
;;; 出現頻度表の作成
(defun make-frequency-ad (xs &optional (sum 0))
(if (null xs)
nil
(cons (list (car xs) 1 sum (1+ sum))
(make-frequency-ad (cdr xs) (1+ sum)))))
;;; 出現頻度表の更新
(defun update-frequency (code freq)
(let ((xs (member code freq :key #'car)))
(incf (second (car xs)))
(incf (fourth (car xs)))
(dolist (x (cdr xs))
(incf (third x))
(incf (fourth x)))))
;;; 符号化
(defun encode-ad (xs)
(let ((table (make-frequency-ad syms))
(sum (length syms))
(low 0)
(high 1))
(dolist (x xs (make-code low high))
(let ((c (find x table :key #'car))
(w (- high low)))
(setq low (+ (* w (/ (third c) sum)) low)
high (+ (* w (/ (second c) sum)) low))
;; 出現頻度表の更新
(incf sum)
(update-frequency x table)))))
;;; 復号
(defun decode-ad (n code)
(let ((table (make-frequency-ad syms))
(sum (length syms))
(value (code-to-number code)) buffer)
(dotimes (i n (reverse buffer))
(let ((data (find-if (lambda (xs)
(and (<= (/ (third xs) sum) value)
(< value (/ (fourth xs) sum))))
table)))
(push (car data) buffer)
(setq value (/ (- value (/ (third data) sum))
(/ (second data) sum)))
(incf sum)
(update-frequency (car data) table)))))