第25章
整数算術符号化
（やさしい版）
Pearls of Functional Algorithm Design（関数プログラミングによるアルゴリズム設計の真珠）
どんな問題？
前章では、算術符号化を「区間を分数（有理数）でどんどん狭めていく」やり方で説明しました。しかし実際のコンピュータで、際限なく大きくなる分数をそのまま扱うのは現実的ではありません。
そこで本章では、分数の代わりに「有限のけた数の整数」だけで算術符号化を実装する方法を組み立てます。ゴールは、Haskell の Int 型のような「桁が限られた整数」でも動く encode と decode を作ることです。
この章で扱うこと
- 区間を整数の対 (ℓ, r) で表す新しい表現方法
- 整数だけを使うと発生する「区間がつぶれる」問題と、その回避策
- 逐次符号化（少しずつビットを出す）と区間拡張（区間を広げて表現しなおす）という 2 つのアイデア
- 復号側での「ストリームの反転」と、それを効率化した最終形
新しい定義
区間を整数の対で表す
今回は区間を、有限精度の整数の対 (ℓ, r) で表現します。ただし 0 ≤ ℓ < r ≤ 2e であり、e はあらかじめ決めた固定の整数です。この対は、単位区間 (0, 1) 上の部分区間 (ℓ/2e, r/2e) を表していると考えます。
モデルから来る区間も同じ形式ですが、こちらは別の整数 d を使い、(p/2d, q/2d) の形で表します。
大事な前提
あとで示すように、d ≤ e−2 でなくてはなりません。つまり d と e は同じ値にはできず、モデル側の精度 d は符号化側の精度 e より少なくとも 2 小さくする必要があります。e, d は、すべての計算が Int の範囲に収まる程度に十分小さく選んでおきます。
狭める操作 ▶ の定義
前章の「区間を狭める」操作 ▷ を、整数版の ▶ に置き換えます。定義は次のとおり。
(▶) :: Interval → Interval → Interval
(ℓ, r) ▶ (p, q) = (ℓ+⌊(r−ℓ)∗p/2d⌋, ℓ+⌊(r−ℓ)∗q/2d⌋)
Dart
// 精度パラメータ e, d は const で持つ（本文では d ≤ e-2 が必要）
const int e = 5;
const int d = 3;
const int e1 = 1 << (e - 2); // 2^(e-2)
const int e2 = 1 << (e - 1); // 2^(e-1)
const int e3 = e1 + e2; // 3∗2^(e-2)
const int e4 = 1 << e; // 2^e
typedef Interval = (int, int);
// (ℓ, r) ▶ (p, q)
Interval narrow(Interval i, Interval j) {
final (l, r) = i;
final (p, q) = j;
final w = r - l;
return (l + ((w * p) >> d), l + ((w * q) >> d));
}
この計算の途中に出てくる最大の整数は 2e+d 程度です（(r−ℓ)∗q がそのくらいまで大きくなり得るため）。これが Int の範囲に収まるように e, d を選んでおけば、区間計算はすべて Int でできます。
ビット取り出しの整数版 ibit
前章では区間を小数として扱い、bit という関数で「小数点以下」を 1 ビットずつ取り出していました。整数版でも同じ考え方で、bit を ibit に置き換えます。
ibit (ℓ, r) | r ≤ 2e−1 = Just (0, (2∗ℓ, 2∗r))
| 2e−1 ≤ ℓ = Just (1, (2∗ℓ−2e, 2∗r−2e))
| otherwise = Nothing
Dart
// Haskell の Maybe (Bit, Interval) を (int, Interval)? で表現。null が Nothing。
(int, Interval)? ibit(Interval i) {
final (l, r) = i;
if (r <= e2) return (0, (2 * l, 2 * r));
if (e2 <= l) return (1, (2 * l - e4, 2 * r - e4));
return null; // まん中 e2 をまたぐ → 何も出せない
}
ibit は、区間を 2e 倍にスケールしたもの上で bit と同じ働きをする関数です。
新しい encode の最初の案
この ibit を使って、encode を次のように書き直せます。
encode1 m = unfoldr ibit · foldl (▶) (0, 2e) · intervals m (25.1)
Dart
// モデル m は「今の状態から記号 x に対応する部分区間 (p, q)/2^d」を返すもの、と抽象化。
// intervals m xs は記号列 xs をそのモデルに通して部分区間の列に変える。
typedef Model<X> = ({
Interval Function(X x) interval, // 記号 x → 部分区間 (p, q)
Model<X> Function(X x) adapt, // 適応（本章では不変モデルを想定）
});
Iterable<Interval> intervals<X>(Model<X> m, List<X> xs) sync* {
var cur = m;
for (final x in xs) {
yield cur.interval(x);
cur = cur.adapt(x);
}
}
// unfoldr ibit
Iterable<int> unfoldIbit(Interval start) sync* {
var i = start;
while (true) {
final r = ibit(i);
if (r == null) return;
yield r.$1;
i = r.$2;
}
}
List<int> encode1<X>(Model<X> m, List<X> xs) {
var i = (0, e4);
for (final j in intervals(m, xs)) {
i = narrow(i, j); // foldl (▶)
}
return unfoldIbit(i).toList(); // → ゼロ幅につぶれると無限ループ
}
手順としては、テキストの記号を区間の列に変え、それを ▶ で最終区間まで狭め、そこから ibit でビット列を取り出します。
しかし、これは動かない！
▷ では起こらなかった現象ですが、▶ で狭め続けると、そのうち区間が空になってしまう（つぶれる）ことがあります。
つぶれる例
e = 5, d = 3 として、モデル m が文字「a」に (3, 5) を、文字「b」に (5, 6) を割り当てるとします。文字列 "bba" を符号化すると:
encode1 m "bba"
= foldl (▶) (0, 32) [(5, 6), (5, 6), (3, 5)]
= foldl (▶) (20, 24) [(5, 6), (3, 5)]
= foldl (▶) (22, 23) [(3, 5)]
= (22, 22)
最後は (22, 22) という幅ゼロの区間になってしまい、そこから unfoldr ibit を走らせると、無限にゴミのようなビット列を吐き続けてしまいます。困りました。
逐次符号化と区間拡張
アイデア 1: 逐次符号化
1 つ目のアイデアは、ビットを最後にまとめて出すのではなく、途中でどんどん出してしまうという「逐次符号化」です。foldl を stream に置き換えます。
encode2 m = stream ibit (▶) (0, 2e) · intervals m (25.2)
Dart
// stream f g s = 「f s が Just(y, s') なら y を出して stream f g s' に続く。
// Nothing なら入力から要素 x を取って s を g s x に進める」
// ここでは f = ibit, g = narrow (▶), s = 区間, 入力 = 部分区間の列。
Iterable<int> encode2<X>(Model<X> m, List<X> xs) sync* {
var s = (0, e4);
final it = intervals(m, xs).iterator;
while (true) {
final r = ibit(s);
if (r != null) {
yield r.$1;
s = r.$2;
} else if (it.moveNext()) {
s = narrow(s, it.current);
} else {
return; // 入力が尽きて ibit も Nothing なら終了
}
}
}
▶ は結合的でないため、これは前章のようなストリーミング条件を満たしていません。それでも、先ほどの "bba" の例を追いなおすと、途中でビット 101 や 101101 が出てくることで区間が (0, 32) にリセットされ、うまくつぶれずに済みます。
encode2 m "bba"
= stream ibit (▶) (0, 32) [(5, 6), (5, 6), (3, 5)]
= stream ibit (▶) (20, 24) [(5, 6), (3, 5)]
= 101 : stream ibit (▶) (0, 32) [(5, 6), (3, 5)]
= 101101 : stream ibit (▶) (0, 32) [(3, 5)]
= 101101 : stream ibit (▶) (12, 20) [ ]
= 101101
ところが、"aab" では話が違います。
encode2 m "aab"
= stream ibit (▶) (0, 32) [(3, 5), (3, 5), (5, 6)]
= stream ibit (▶) (12, 20) [(3, 5), (5, 6)]
= stream ibit (▶) (15, 17) [(5, 6)]
= stream ibit (▶) (16, 16) [ ]
= 0111...
ここでは各中間区間が中点 16 をまたいでおり、ibit は「まん中をまたぐ区間」については何も出せない（Nothing を返す）ので、区間が (16, 16) までつぶれてしまいます。
アイデア 2: 区間拡張
2 つ目のアイデアは、まん中 2e−1 をまたぐような細い区間を、いったん広げてから狭めることです。狭めても区間がつぶれないための条件を計算すると:
⌊(r−ℓ)∗p/2d⌋ < ⌊(r−ℓ)∗q/2d⌋
これがすべての p < q について成り立つ条件は、実は簡単で、2d ≤ r−ℓ であることです。
重要な観察
区間の幅を最低 2e−2 以上に保てば、つぶれは起こりません。そのためには d ≤ e−2 が必要です。これが最初に予告した条件の理由です。
区間拡張
区間の三つ組表現
区間 (ℓ, r) を、三つ組 (n, (ℓ′, r′)) の形で表します。ここで n は「何回広げたか」を表し、ℓ′ = widen n ℓ, r′ = widen n r とします。
widen n x = 2n(x − 2e−1) + 2e−1
この widen は、まん中 2e−1 を中心にして幅を 2n 倍に引き伸ばす写像です。完全に拡張された区間とは、0 ≤ ℓ′ < r′ ≤ 2e の範囲を守りつつ、n をできるだけ大きく取ったもののことです。
具体例
e = 5 のとき、幅 4 の区間 (13, 17) は、幅 16 の完全拡張版 (2, (4, 20)) と等価に表せます。
拡張と縮小の関数
記述を簡単にするため、四つの整数 e1, e2, e3, e4 を ei = (i/4)·2e と定めます。すると e2 がちょうどまん中で、e4 = 2e です。
もう一段拡張できる条件は e1 ≤ ℓ かつ r ≤ e3（区間が中央よりだいぶ内側にあるとき）と表せて、expand は次のように書けます。
expand i = extend (0, i)
extend (n, (ℓ, r))
| e1 ≤ ℓ ∧ r ≤ e3 = extend (n+1, 2∗ℓ−e2, 2∗r−e2)
| otherwise = (n, (ℓ, r))
Dart
typedef Extended = (int, Interval); // (n, (ℓ, r))
Extended expand(Interval i) => extend((0, i));
Extended extend(Extended ei) {
var (n, i) = ei;
var (l, r) = i;
while (e1 <= l && r <= e3) { // まん中を含みつつ内側にあるなら拡張
n += 1;
l = 2 * l - e2;
r = 2 * r - e2;
}
return (n, (l, r));
}
逆の縮小操作 contract は次のとおり。
contract (n, (ℓ, r)) = (shorten n ℓ, shorten n r)
Dart
// shorten n x = (x − e2) / 2^n + e2 （負値もあるので ~/ を使う）
int shorten(int n, int x) => (x - e2) ~/ (1 << n) + e2;
Interval contract(Extended ei) {
final (n, (l, r)) = ei;
return (shorten(n, l), shorten(n, r));
}
ここで shorten n x = (x − e2)/2n + e2 です。contract · expand = id ですが、逆の expand · contract は一般には恒等ではありません（データ精緻化でよくある構図）。
enarrow: 拡張してから狭める
enarrow（extend + narrow の略）は、部分拡張された区間をいったん完全拡張してから ▶ で狭めます。
enarrow :: (Int, Interval) → Interval → (Int, Interval)
enarrow ei j = (n, i ▶ j) where (n, i) = extend ei
Dart
Extended enarrow(Extended ei, Interval j) {
final (n, i) = extend(ei);
return (n, narrow(i, j)); // 拡張してから ▶ で狭める
}
この工夫のおかげで、▶ で狭められるのは幅の十分ある区間だけになり、つぶれが避けられます。
新しい定義
encode の三度目の書き直し
ここまでのアイデアを合わせて、encode の新しい定義に置き換えます。
encode3 m = stream ebit enarrow (0, (0, 2e)) · intervals m (25.3)
Dart
// stream 化した encode。状態は Extended（拡張回数付き区間）。
Iterable<int> encode3Bits<X>(Model<X> m, List<X> xs) sync* {
Extended s = (0, (0, e4));
final it = intervals(m, xs).iterator;
while (true) {
final r = ebit(s);
if (r != null) {
yield r.$1;
s = r.$2;
} else if (it.moveNext()) {
s = enarrow(s, it.current);
} else {
return;
}
}
}
ここで ebit は、拡張された区間の上で動作する ibit の対応物です。次の性質を満たすように設計します。
unfoldr ebit = unfoldr ibit · contract (25.4)
ebit は「e2 をまたぐ区間」については Nothing を返すので、▶ で狭めるのは幅が e1 以上のときだけになります。これで d ≤ e−2 の下でのつぶれ回避が保証されます。
ebit を融合則で導く
(25.4) は unfoldr の融合則の適用を示唆しています。融合則とは、次の 2 条件が満たされれば unfoldr h = unfoldr f · g となる、という定理です。
h x = Nothing ⇒ f (g x) = Nothing
h x = Just (y, x′) ⇒ f (g x) = Just (y, g x′)
h = ebit, f = ibit, g = contract と置いて、条件を満たすように ebit を定めると次のようになります。
ebit (0, (ℓ, r))
| r ≤ e2 = Just (0, (0, (2∗ℓ, 2∗r)))
| e2 ≤ ℓ = Just (1, (0, (2∗ℓ−e4, 2∗r−e4)))
| otherwise = Nothing
ebit (n+1, (ℓ, r))
| r ≤ e2 = Just (0, (n, (ℓ+2n∗e2, r+2n∗e2)))
| e2 ≤ ℓ = Just (1, (n, (ℓ−2n∗e2, r−2n∗e2)))
| otherwise = Nothing
Dart
// 融合則から導かれた ebit。1 ビットずつ取り出す版。
(int, Extended)? ebit(Extended ei) {
final (n, (l, r)) = ei;
if (n == 0) {
if (r <= e2) return (0, (0, (2 * l, 2 * r)));
if (e2 <= l) return (1, (0, (2 * l - e4, 2 * r - e4)));
return null;
} else {
final k = 1 << (n - 1); // 2^(n-1) （n' = n-1 として 2^n' ∗ e2 の代わりに k*e2）
if (r <= e2) return (0, (n - 1, (l + k * e2, r + k * e2)));
if (e2 <= l) return (1, (n - 1, (l - k * e2, r - k * e2)));
return null;
}
}
この定義でも動きますが、少しごちゃごちゃしていて非効率です。
ebits: 1 ステップでまとめて出す
r ≤ e2 の場合をよく観察すると、unfoldr ebit は「まず 0 が 1 個」「次に 1 が n 個連続」という形の出力を続けます。だったら、これを 1 ステップでまとめて出せば効率的です。
ebits (n, (ℓ, r))
| r ≤ e2 = Just (bits n 0, (0, (2∗ℓ, 2∗r)))
| e2 ≤ ℓ = Just (bits n 1, (0, (2∗ℓ−e4, 2∗r−e4)))
| otherwise = Nothing
Dart
// bits n b = [b, 1-b, 1-b, ..., 1-b] （b のあと 1-b が n 個）
List<int> bits(int n, int b) => [b, ...List.filled(n, 1 - b)];
// 1 ステップで n+1 ビットまとめて出す版。
(List<int>, Extended)? ebits(Extended ei) {
final (n, (l, r)) = ei;
if (r <= e2) return (bits(n, 0), (0, (2 * l, 2 * r)));
if (e2 <= l) return (bits(n, 1), (0, (2 * l - e4, 2 * r - e4)));
return null;
}
ここで bits n b = b : replicate n (1−b) は「b のあとに 1−b を n 個」を返す関数です。これによって、(25.3) は同値でより効率的な次の版になります。
encode3 m
= concat · stream ebits enarrow (0, (0, 2e)) · intervals m (25.5)
Dart
// concat · stream ebits enarrow ... の最終形。
Iterable<int> encode3<X>(Model<X> m, List<X> xs) sync* {
Extended s = (0, (0, e4));
final it = intervals(m, xs).iterator;
while (true) {
final r = ebits(s);
if (r != null) {
yield* r.$1; // まとめて出したビット列を平坦化（= concat）
s = r.$2;
} else if (it.moveNext()) {
s = enarrow(s, it.current);
} else {
return;
}
}
}
これが encode の最終形です。
肝心な問い
encode3 は結局なにをしている？
前章の encode は「出力ビット列を分数に戻すと、記号列から得られる区間の中に入る」という性質を満たしていました。encode3 ではこの性質はそのままの形では成り立ちません。
そこで、開始区間を追加の引数に持つ変種 encode′3 を定義します。
encode′3 m ei xs = concat · stream ebits enarrow ei · intervals m xs
Dart
// 開始区間 ei を引数に取れるようにした版（性質 (25.6) の主役）。
Iterable<int> encode3Prime<X>(Model<X> m, Extended ei, List<X> xs) sync* {
var s = ei;
final it = intervals(m, xs).iterator;
while (true) {
final r = ebits(s);
if (r != null) {
yield* r.$1;
s = r.$2;
} else if (it.moveNext()) {
s = enarrow(s, it.current);
} else {
return;
}
}
}
すると、次の性質が任意のモデル m、拡張区間 ei、記号列 xs について成り立ちます。
2e ∗ toFrac (encode′3 m ei xs) ∈ contract ei (25.6)
なぜ大事？
この性質 (25.6) は、decode を作るときの土台になります。「符号化された結果を分数に戻したものが、開始区間を縮めた領域に入っている」ことを言っているのです。証明は付録にあります。
最後の難点
拡張回数 n の暴走
実は (25.5) にも弱点があります。完全拡張された区間 (n, i) の n（拡張回数）がとても大きくなり、有限精度整数の範囲を超えてしまうことが理論的にはあり得るのです。
たとえば、(0, e4) を (3/8, 5/8) のような「まん中をまたぐ細い区間」で何度も何度も狭めていく状況を想像してみましょう。ebit は何も出さないまま、拡張だけがどんどん進み、n がふくれあがります。
現実的なまとめ
この状況は実際にはほとんど起こりませんが、論理的には起こり得ます。要は、算術符号化のいかなる版も、有限精度算術での完全な動作は保証できないのです。もし起きたら、エラーで中断するか有理数版に切り替える、という選択になります。
ストリームの反転
decode の方針
decode は xs ⊑ decode m (encode m xs) を満たすように作ります。encode が stream を使って作られている以上、decode を作るには「ストリームを反転する」道具が要ります。それが destream です。
destream f g h s ys = unfoldr step (s, ys)
where step (s, ys) = case f s of
Just (y, s′) → step (s′, ys ↓ [y])
Nothing → Just (x, (g s x, ys))
where x = h s ys
Dart
// stream の双対。f/g/h は下で ebit/enarrow/助手 h を渡すために抽象化。
// ys ↓ [y] は「ys の先頭が y なら 1 つ取り除く」に相当。
Iterable<X> destream<S, Y, X>(
(Y, S)? Function(S) f,
S Function(S, X) g,
X Function(S, List<Y>) h,
S s0,
List<Y> ys0,
) sync* {
var s = s0;
var ys = ys0;
while (true) {
final r = f(s);
if (r != null) {
// f が y を出す → ys の先頭から y を 1 つ削る
s = r.$2;
ys = ys.sublist(1);
} else {
// f が Nothing → h で出力の要素 x を作る
final x = h(s, ys);
yield x;
s = g(s, x);
}
}
}
ここで演算 ↓ は (us ⧺ vs) ↓ us = vs（先頭の us を取り除く）で定義されます。動作イメージは stream の双対です。
- f s が値 y を返すとき: 入力ストリーム ys の先頭から y を取り除く。
- f s が Nothing のとき: 「補助関数」 h を使って、出力の要素を 1 つ生み出す。
反ストリーミング定理
stream と destream の関係は、次の定理で結ばれます。
定理 25.1（反ストリーミング定理）
stream f g s xs が有限リストを返し、f s = Nothing のとき補助関数 h が h s (stream f g s (x : xs)) = x を満たすと仮定する。このとき
xs ⊑ destream f g h s (stream f g s xs)
この定理を f = ebit, g = enarrow, s = ei と置いて適用すると、次の decode が得られます。
decode m = destream ebit enarrow h (m, (0, (0, e4)))
Dart
// 状態は (モデル m, Extended)。ここではモデル型 X 固定で示す。
// destream の型に合わせるため、状態を (Model, Extended) にまとめる。
typedef DState<X> = (Model<X>, Extended);
Iterable<X> decodeDestream<X>(
Model<X> m,
List<int> bs,
X Function(DState<X>, List<int>) hAux,
) {
return destream<DState<X>, int, X>(
(s) {
final r = ebit(s.$2);
if (r == null) return null;
return (r.$1, (s.$1, r.$2));
},
(s, x) => (s.$1.adapt(x), enarrow(s.$2, s.$1.interval(x))),
hAux,
(m, (0, (0, e4))),
bs,
);
}
ただし、補助関数 h が次を満たすことが条件です（ei が e2 をまたぐすべての場合について）。
h (m, ei) (encode′3 m ei (x : xs)) = x (25.7)
ebits 版に書き直す
ebit を ebits に置き換えたのと同じ要領で、decode も次の形になります。
decode m bs = unfoldr step (m, (0, (0, e4)), bs)
step (m, (n, (ℓ, r)), bs)
| r ≤ e2 = step (m, (0, (2∗ℓ, 2∗r)), bs ↓ bits n 0)
| e2 ≤ ℓ = step (m, (0, (2∗ℓ−e4, 2∗r−e4)), bs ↓ bits n 1)
| otherwise = Just (x, (adapt m x,
enarrow (n, (ℓ, r)) (interval m x), bs))
where x = h (m, (n, (ℓ, r))) bs
Dart
// ebits 版に書き直した decode。1 ステップで n+1 ビットまとめて剥がす。
Iterable<X> decodeEbits<X>(
Model<X> m0,
List<int> bs0,
X Function(Model<X>, Extended, List<int>) hAux,
) sync* {
var m = m0;
var ei = (0, (0, e4));
var bs = bs0;
while (true) {
final (n, (l, r)) = ei;
if (r <= e2) {
ei = (0, (2 * l, 2 * r));
bs = bs.sublist(n + 1); // bs ↓ bits n 0
} else if (e2 <= l) {
ei = (0, (2 * l - e4, 2 * r - e4));
bs = bs.sublist(n + 1); // bs ↓ bits n 1
} else {
final x = hAux(m, ei, bs);
yield x;
ei = (n, narrow((l, r), m.interval(x))); // enarrow の展開結果を再構成
m = m.adapt(x);
}
}
}
あとは補助関数 h を実際に作ることが残された仕事です。
補助関数
◀ の導出
まず、▷ に対する ◁ と同じ役割を、▶ に対して果たす演算 ◀ を作ります。使うのは床の規則: 整数 n と実数 f について n ≤ f ≡ n ≤ ⌊f⌋ という性質です。
計算をたどると次に到達します。
(◀) :: Int → Interval → Int
k ◀ (ℓ, r) = ((k − ℓ + 1)∗2d − 1) div (r − ℓ)
Dart
// k ∈ (i ▶ j) ≡ (k ◀ i) ∈ j という関係を満たす逆演算。
int nnarrow(int k, Interval i) {
final (l, r) = i;
return ((k - l + 1) * (1 << d) - 1) ~/ (r - l);
}
これで k ∈ (i ▶ j) ≡ (k ◀ i) ∈ j が成立します。
h の定義
性質 (25.6) を widen/contract の言葉で書き直し、床の規則を使ってさらに整理すると、記号 x は次のように取り出せることがわかります。
h (m, ei) bs = symbol m (⌊widen n (2e ∗ toFrac bs)⌋ ◀ i)
where (n, i) = extend ei
Dart
// symbol はモデル上で「区間内の整数 k を、それを含む記号に写す」関数
// （ここでは Model に symbol を含めず、探索は呼び出し側で行う想定）。
typedef ModelSym<X> = ({
Interval Function(X x) interval,
ModelSym<X> Function(X x) adapt,
X Function(int k) symbol,
});
// toFrac bs = 0.b0 b1 b2... を有理数として見る（e4 倍して整数化）
// widen n x = 2^n(x − e2) + e2
int widen(int n, int x) => (1 << n) * (x - e2) + e2;
// bs から先頭 e ビットを整数化した値をベースに、e4∗toFrac(bs) の床を近似。
int toIntBits(List<int> bs, int k) {
var v = 0;
for (var i = 0; i < k; i++) {
v = (v << 1) | (i < bs.length ? bs[i] : 0);
}
return v;
}
X hAux<X>(ModelSym<X> m, Extended ei, List<int> bs) {
final (n, i) = extend(ei);
final f = widen(n, toIntBits(bs, e)); // ⌊widen n (e4∗toFrac bs)⌋
return m.symbol(nnarrow(f, i));
}
逐次復号
いまの decode の弱点
今の decode にはたくさんの弱点があります。
- 有理数算術を使う（toFrac bs は分数）
- 逐次的ではない
- extend の計算が h と enarrow で重複している
- toFrac を毎回計算しなおしている
- widen n はコストの高い指数計算を含む
これらを解決するために、decode を3 段階で逐次化します。
段階 1: extend を step に取り込む
step の中で拡張の場合分けもいっしょに扱うことで、外側の extend を追い出します。
step (m, (n, (ℓ, r)), bs)
| r ≤ e2 = step (m, (0, (2∗ℓ, 2∗r)), bs ↓ bits n 0)
| e2 ≤ ℓ = step (m, (0, (2∗ℓ−e4, 2∗r−e4)), bs ↓ bits n 1)
| e1 ≤ ℓ ∧ r ≤ e3 = step (m, (n+1, (2∗ℓ−e2, 2∗r−e2)), bs)
| otherwise = Just (x,
(adapt m x, (n, (ℓ, r) ▶ interval m x), bs))
where x = symbol m (⌊widen n (e4∗toFrac bs)⌋ ◀ (ℓ, r))
Dart
// 段階 1: 拡張の場合分けを step 内に取り込む。extend の呼び出しを排除。
Iterable<X> decodeStage1<X>(ModelSym<X> m0, List<int> bs0) sync* {
var m = m0;
var ei = (0, (0, e4));
var bs = bs0;
while (true) {
final (n, (l, r)) = ei;
if (r <= e2) {
ei = (0, (2 * l, 2 * r));
bs = bs.sublist(n + 1);
} else if (e2 <= l) {
ei = (0, (2 * l - e4, 2 * r - e4));
bs = bs.sublist(n + 1);
} else if (e1 <= l && r <= e3) {
ei = (n + 1, (2 * l - e2, 2 * r - e2)); // 拡張を 1 段進める
} else {
final f = widen(n, toIntBits(bs, e));
final x = m.symbol(nnarrow(f, (l, r)));
yield x;
ei = (n, narrow((l, r), m.interval(x)));
m = m.adapt(x);
}
}
}
段階 2: toFrac の繰り返しをやめる
「毎回計算しなおしている f = widen n (e4∗toFrac bs) の値そのもの」を状態として引きずり回します。次の関係が成り立つのがカギです。
f 0 (bs ↓ bits n b) = 2∗f n bs − e4∗b
f (n+1) bs = 2∗f n bs − e2
結果は次のとおり。
decode m bs = unfoldr step (m, (0, (0, e4)), e4∗toFrac bs) (25.8)
step (m, (n, (ℓ, r)), f)
| r ≤ e2 = step (m, (0, (2∗ℓ, 2∗r)), 2∗f)
| e2 ≤ ℓ = step (m, (0, (2∗ℓ−e4, 2∗r−e4)), 2∗f−e4)
| e1 ≤ ℓ ∧ r ≤ e3 = step (m, (n+1, (2∗ℓ−e2, 2∗r−e2)), 2∗f−e2)
| otherwise = Just (x,
(adapt m x, (n, (ℓ, r) ▶ interval m x), f))
where x = symbol m (⌊f⌋ ◀ (ℓ, r))
Dart
// 段階 2: 状態に f = ⌊widen n (e4∗toFrac bs)⌋ を持たせ、更新は加減算とシフトで。
// （ここではまだ f を「有理数の床値」扱いだが、値は int で保持できる。）
Iterable<X> decodeStage2<X>(ModelSym<X> m0, List<int> bs) sync* {
var m = m0;
var ei = (0, (0, e4));
var f = widen(0, toIntBits(bs, e)); // 初期値: e4∗toFrac bs の床
while (true) {
final (n, (l, r)) = ei;
if (r <= e2) {
ei = (0, (2 * l, 2 * r));
f = 2 * f;
} else if (e2 <= l) {
ei = (0, (2 * l - e4, 2 * r - e4));
f = 2 * f - e4;
} else if (e1 <= l && r <= e3) {
ei = (n + 1, (2 * l - e2, 2 * r - e2));
f = 2 * f - e2;
} else {
final x = m.symbol(nnarrow(f, (l, r)));
yield x;
ei = (n, narrow((l, r), m.interval(x)));
m = m.adapt(x);
}
}
}
さらに、n は使わなくてよい冗長な変数と分かるので落とせます。
step (m, (ℓ, r), f)
| r ≤ e2 = step (m, (2∗ℓ, 2∗r), 2∗f)
| e2 ≤ ℓ = step (m, (2∗ℓ−e4, 2∗r−e4), 2∗f−e4)
| e1 ≤ ℓ ∧ r ≤ e3 = step (m, (2∗ℓ−e2, 2∗r−e2), 2∗f−e2)
| otherwise = Just (x,
(adapt m x, (ℓ, r) ▶ interval m x, f))
where x = symbol m (⌊f⌋ ◀ (ℓ, r))
Dart
// n は使わない冗長変数と分かるので落とす。区間 (ℓ, r) と f だけで回る版。
Iterable<X> decodeStage2NoN<X>(ModelSym<X> m0, List<int> bs) sync* {
var m = m0;
var i = (0, e4);
var f = widen(0, toIntBits(bs, e));
while (true) {
final (l, r) = i;
if (r <= e2) {
i = (2 * l, 2 * r);
f = 2 * f;
} else if (e2 <= l) {
i = (2 * l - e4, 2 * r - e4);
f = 2 * f - e4;
} else if (e1 <= l && r <= e3) {
i = (2 * l - e2, 2 * r - e2);
f = 2 * f - e2;
} else {
final x = m.symbol(nnarrow(f, i));
yield x;
i = narrow(i, m.interval(x));
m = m.adapt(x);
}
}
}
段階 3: 逐次計算にする
e4 = 2e なので、⌊e4∗toFrac bs⌋ は bs の先頭 e ビットだけから決まります。実際:
⌊e4∗toFrac bs⌋ = toInt (take e (bs ⧺ repeat 0))
ここで toInt は「ビット列を整数に読み替える」関数です。さらに bs′ = bs ⧺ repeat 0 とすれば:
⌊2∗e4∗toFrac bs⌋ = 2∗toInt (take e bs′) + head (drop e bs′)
これを使うと、状態 f を「整数 n（先頭 e ビットを整数化したもの）＋残りビット列 ds」に置き換えて、完全に有限精度整数だけで逐次計算できます。
decode m bs = unfoldr step (m, (0, e4), toInt (take e bs′), drop e bs′)
where bs′ = bs ⧺ repeat 0
step (m, (ℓ, r), n, b : bs)
| r ≤ e2 = step (m, (2∗ℓ, 2∗r), 2∗n+b, bs)
| e2 ≤ ℓ = step (m, (2∗ℓ−e4, 2∗r−e4), 2∗n−e4 + b, bs)
| e1 ≤ ℓ ∧ r ≤ e3 = step (m, (2∗ℓ−e2, 2∗r−e2), 2∗n−e2 + b, bs)
| otherwise = Just (x,
(adapt m x, (ℓ, r) ▶ interval m x, n, b : bs))
where x = symbol m (n ◀ (ℓ, r))
Dart
// 段階 3: 状態 f を「先頭 e ビットを整数化した n」＋「残りビット列 bs」に置換。
// これで完全に有限精度整数だけで逐次計算できる（図 25.1 の最終形）。
Iterable<X> decodeFinal<X>(ModelSym<X> m0, List<int> bs0) sync* {
var m = m0;
var i = (0, e4);
// bs' = bs ⧺ repeat 0 の代わりに、必要になったら 0 を補うヘルパで代用。
var n = toIntBits(bs0, e);
var idx = e; // bs0 内の次に読むビット位置
int nextBit() {
final b = idx < bs0.length ? bs0[idx] : 0;
idx += 1;
return b;
}
while (true) {
final (l, r) = i;
if (r <= e2) {
i = (2 * l, 2 * r);
n = 2 * n + nextBit();
} else if (e2 <= l) {
i = (2 * l - e4, 2 * r - e4);
n = 2 * n - e4 + nextBit();
} else if (e1 <= l && r <= e3) {
i = (2 * l - e2, 2 * r - e2);
n = 2 * n - e2 + nextBit();
} else {
final x = m.symbol(nnarrow(n, i));
yield x;
i = narrow(i, m.interval(x));
m = m.adapt(x);
}
}
}
図 25.1 decode の最終形
Dart
// 動作確認: 本文の例に沿って、記号 'a'→(3,5), 'b'→(5,6) を割り当てた
// 静的モデル（adapt は自身を返す）で encode3 を走らせる。
void main() {
final staticModel = <X>(Map<X, Interval> table) {
late Model<X> m;
m = (
interval: (X x) => table[x]!,
adapt: (X _) => m,
);
return m;
};
final m = staticModel<String>({'a': (3, 5), 'b': (5, 6)});
final encoded = encode3(m, ['b', 'b', 'a']).toList();
print(encoded); // 例: [1, 0, 1, 1, 0, 1]
}
結びの覚え書き
ここまでお付き合いいただいた読者は、「算術符号化」には本当に色々な算術が絡んでいることを実感されたと思います。数値の計算だけでなく、foldl や unfoldr といった折りたたみ・展開の算術も含めた、二重の意味での「arithmetic」でした。
まとめ
算術符号化のアイデア自体は単純だが、有限精度の整数で実装しようとすると、区間のつぶれを避ける工夫（逐次符号化と区間拡張）や、復号側の逐次化のための丁寧な変形が必要になる。
付録
反ストリーミング定理の証明の骨格
証明は、stream と destream の定義から出る次の 2 つの性質に基づきます。
f s = Nothing ⇒
stream f g s (x : xs) = stream f g (g s x) xs ∧
destream f g h s ys = x : destream f g h (g s x) ys
f s = Just (y, s′) ⇒
stream f g s xs = y : stream f g s′ xs ∧
destream f g h s (y : ys) = destream f g h s′ ys
示したい主張は次のとおり。
xs ⊑ destream f g h s (stream f g s xs)
証明は xs と n（stream f g s xs の長さ）についての二重帰納法です。
- [ ] の場合: [ ] は任意のリストの接頭辞なので、直ちに成立。
- x : xs の場合 かつ f s = Nothing: 上の第一性質から、両辺の先頭が x となり、残りの部分に帰納法を適用できる。
- x : xs の場合 かつ f s = Just (y, s′): 上の第二性質から、stream の長さが 1 減るので、こちらも帰納法で決着する。
性質 (25.6) の証明の骨格
示したいのは e4∗toFrac (encode′3 m ei xs) ∈ contract ei です。xs と n（encode′3 m ei xs の長さ）についての二重帰納法で示します。
[ ] の場合: encode′3 m ei [ ] = concat (unfoldr ebits ei) に (25.4) と ibit の定義を順に適用すると true に落ちる。
x : xs の場合: encode′3 の場合分けごとに帰納法を回す。狭める節では、次の包含関係がカギとなる。
contract (enarrow ei j) ⊆ contract (extend ei) = contract ei
ビットを出す節では、対応する開始区間の長さが厳密に短くなることを使って帰納法を適用する。最後に toFrac (bits n b ⧺ bs) = (2n + (b−1) + toFrac bs)/2n+1 という簡単な等式で締めくくれば、証明が完了する。