算术编码与 ANS：超越 Huffman

一、从 Huffman 的天花板说起

1952 年，David Huffman 在 MIT 的课堂作业中提出了 Huffman 编码。七十多年过去了，它依然是最广为人知的熵编码方案。但 Huffman 编码有一个根本性的局限：每个符号至少编码为 1 比特。

这意味着什么？考虑一个极端但常见的情况。假设某个符号的出现概率为 0.95，它携带的信息量为：

\[ I = -\log_2 0.95 \approx 0.074 \text{ bit} \]

理论上只需要 0.074 比特就能编码这个符号。但 Huffman 编码最短的码字就是 1 比特，对这个符号而言，实际编码长度是理论最优值的约 \(1/0.074 \approx 13.5\) 倍。

这不是理论上的吹毛求疵。在现实中，这种情况非常普遍：

• LZ77 系列压缩：match/literal 标志位（概率极度不均衡）
• 图像编码：预测残差集中在零附近
• 二进制算术编码：上下文建模后的二值化符号
• 通用数据压缩：经过上下文建模后，许多符号概率极高

这里要分清两个口径，否则容易得出互相矛盾的数字：

• 单个符号的浪费：上面那个 \(p=0.95\) 的符号，理论 0.074 比特、实际 1 比特，浪费比例 \((1-0.074)/1 \approx 92.6\%\)。
• 整个二元信源的冗余：一个只有两个取值、\(p=0.95\) 的信源，每符号平均要付出的代价是

\[ \text{冗余} = 1 - H(p),\qquad H(p) = -p\log_2 p - (1-p)\log_2(1-p) \]

当 \(p=0.95\) 时 \(H(p)\approx 0.286\) 比特，Huffman 每符号固定花 1 比特，冗余约 \(1-0.286 = 0.714\) 比特/符号（71.4%）。这就是 Huffman 的 “1 比特粒度” 问题：无论概率多极端，每个符号至少 1 比特。

上图把两种口径画在一起：左边是单符号必须付满 1 比特；右边是算术编码 / ANS 把分数比特摊到整条消息上。底栏则是整个二元信源的平均代价——Huffman 钉死在 1.0，AC/ANS 贴住熵 0.286。

缓解但无法根治

有几种方法可以缓解这个问题：

分组编码（Block Coding）：将 \(n\) 个符号分为一组联合编码，字母表大小从 \(k\) 膨胀到 \(k^n\)，码表指数增长。实践中 \(n\) 很难取大。

游程编码（Run-Length Coding）：对连续重复符号先做游程统计，再 Huffman 编码。但这引入了额外的建模假设，且只对特定分布有效。

下表是 Huffman “离散码字 vs 连续信息量” 这一核心矛盾的直观体现（“浪费” 指单符号在 Huffman 最短码字下相对其信息量的多付比例）：

根本的解决方案是：彻底放弃为每个符号分配整数比特码字的思路。这就是算术编码的核心洞察。

二、算术编码：区间细分的艺术

算术编码的核心思想可以用一句话概括：将整个消息编码为 \([0, 1)\) 上的一个子区间，最终输出落在该子区间内的任意一个数（tag）的二进制表示。

消息越长、符号越稀有，区间就越窄；区间宽度的负对数，就是需要写出的比特数。于是”每个符号付整数比特”被换成了”整条消息共享一条越来越细的数轴”。

基本原理：跟着 “BAC” 走一遍

三符号信源 \(\{A, B, C\}\)，概率 \(P(A)=0.6,\ P(B)=0.3,\ P(C)=0.1\)。先把 \([0,1)\) 按概率切成三段：

编码 “BAC” 时，每读一个符号，就在当前区间里按同样比例再切一刀，只保留该符号对应的那一段：

对应的数值：

\[ \begin{aligned} \text{初始}&: [0.0,\ 1.0), &\text{宽度}&=1.0 \\ \text{编码 B}&: [\,0.0+1.0\cdot0.6,\ 0.0+1.0\cdot0.9\,) = [0.6,\ 0.9), &\text{宽度}&=0.3 \\ \text{编码 A}&: [\,0.6+0.3\cdot0.0,\ 0.6+0.3\cdot0.6\,) = [0.6,\ 0.78), &\text{宽度}&=0.18 \\ \text{编码 C}&: [\,0.6+0.18\cdot0.9,\ 0.6+0.18\cdot1.0\,) = [0.762,\ 0.78), &\text{宽度}&=0.018 \end{aligned} \]

最后一行在图里做了放大：当前区间已经缩到 \([0.6, 0.78)\)，再切出 C 的 \([0.762, 0.78)\)。任意落在最终区间内的 tag（例如 \(v=0.762\)）都能唯一代表 “BAC”。宽度 \(0.018 = 0.3\cdot0.6\cdot0.1\)，所需比特约 \(-\log_2 0.018 \approx 5.8\)。

数学公式化

设累积分布 \(C(s)\)，符号概率 \(P(s)\)。编码器维护区间 \([\text{low}, \text{high})\)，对每个输入符号 \(s\)：

\[ \begin{aligned} \text{range} &= \text{high} - \text{low} \\ \text{high} &= \text{low} + \text{range}\cdot\big(C(s)+P(s)\big) \\ \text{low} &= \text{low} + \text{range}\cdot C(s) \end{aligned} \]

最终宽度等于各符号概率之积，所需比特数等于消息的自信息之和：

\[ \text{width} = \prod_{i=1}^{n} P(s_i),\qquad \text{bits} = -\log_2 \text{width} = -\sum_{i=1}^{n}\log_2 P(s_i) = \sum_{i=1}^{n} I(s_i) \]

这就是算术编码能逼近香农极限的原因：代价按真实概率连续累加，不再被整数码长卡住。

解码：把 tag 逐步还原成符号

解码是编码的逆过程。已知 tag \(v=0.762\)，每一步看它落在哪一段，输出对应符号，再把该段拉伸回 \([0,1)\)：

\[ \begin{aligned} v=0.762 &\in [0.6, 0.9) \Rightarrow \text{B},& v &\leftarrow (0.762-0.6)/0.3 = 0.54 \\ v=0.54 &\in [0.0, 0.6) \Rightarrow \text{A},& v &\leftarrow (0.54-0.0)/0.6 = 0.9 \\ v=0.9 &\in [0.9, 1.0) \Rightarrow \text{C},& v &\leftarrow (0.9-0.9)/0.1 = 0.0 \end{aligned} \]

拉伸公式 \(v \leftarrow (v - C(s)) / P(s)\) 与编码时的区间收缩互逆：编码把世界缩小到子区间，解码把子区间重新放大成世界。输出 “BAC”，与输入一致。

三、二进制小数表示与重归一化

上一节的 “教科书版” 有一个致命问题：区间越来越窄，精度需求随消息长度线性增长。编码 1MB 文件，理论上要数百万位浮点精度——不可用。

重归一化（Renormalization）

工程解法是边编码边吐比特：一旦区间完全落在左半或右半，最高位就已经定了，立刻写出该比特，再把区间放大一倍，把精度”借回来”。

区间 \([\text{low}, \text{high})\) 只有三种局面：

1. 完全在左半（\(\text{high} \le 0.5\)）：最高位必为 0。输出 0，再加倍：\(\text{low} \leftarrow 2\,\text{low}\)，\(\text{high} \leftarrow 2\,\text{high}\)。
2. 完全在右半（\(\text{low} \ge 0.5\)）：最高位必为 1。输出 1，减掉 \(0.5\) 再加倍：\(\text{low} \leftarrow 2(\text{low}-0.5)\)，\(\text{high} \leftarrow 2(\text{high}-0.5)\)。
3. 卡在中间（\(\text{low} < 0.5 \le \text{high}\)，且落在 \([0.25, 0.75)\)）：左右都还可能，比特暂时不定。记一次 underflow（pending++），以 \(0.25\) 为中心放大：\(\text{low} \leftarrow 2(\text{low}-0.25)\)，\(\text{high} \leftarrow 2(\text{high}-0.25)\)。等以后情况一或二终于定下方向时，再把积压的 pending 比特按相反极性一并刷出。

直观理解：情况一、二是”已经知道下一个二进制小数位”；情况三是”还差一口气，先把区间拉开，欠着的比特稍后补”。没有重归一化，算术编码只能停在纸上；有了它，固定位宽的整数寄存器就够用了。

整数实现

实际实现中，我们使用整数算术来避免浮点精度问题。设精度为 N 位（通常 N=32 或 N=16）：

#define PRECISION  32
#define WHOLE      (1ULL << PRECISION)      // 1.0
#define HALF       (1ULL << (PRECISION-1))  // 0.5
#define QUARTER    (1ULL << (PRECISION-2))  // 0.25

// 编码过程（整数版）
void encode_symbol(uint32_t cum_freq, uint32_t sym_freq, uint32_t total_freq) {
    uint64_t range = (uint64_t)(high - low) + 1;
    high = low + (range * (cum_freq + sym_freq)) / total_freq - 1;
    low  = low + (range * cum_freq) / total_freq;

    // 重归一化循环
    for (;;) {
        if (high < HALF) {
            output_bit(0);
            while (pending > 0) { output_bit(1); pending--; }
        } else if (low >= HALF) {
            output_bit(1);
            while (pending > 0) { output_bit(0); pending--; }
            low  -= HALF;
            high -= HALF;
        } else if (low >= QUARTER && high < 3 * QUARTER) {
            pending++;
            low  -= QUARTER;
            high -= QUARTER;
        } else {
            break;
        }
        low  = 2 * low;
        high = 2 * high + 1;
    }
}

这段代码的关键在于 pending 计数器，它处理了 “快要收敛但还没收敛” 的尴尬情况。当区间横跨中点但集中在中间时，我们无法立即输出比特，但可以记录下 “一旦确定方向，需要输出多少个相反的比特”。

四、Range Coder：无专利的实用替代

算术编码的一个重大问题是专利。IBM 在 1970 年代末到 1980 年代对算术编码的核心技术申请了大量专利。这使得很多开源项目和商业软件不敢使用算术编码。

早在 1979 年，IBM 英国科学中心的 G. Nigel N. Martin 就在 Southampton 的 Video & Data Recording 会议上提出了 Range Coder（区间编码器）。它在数学上等价于算术编码（本质是 Pasco 1976 年提出的 FIFO 算术码的再发现），但以字节而非比特为单位进行重归一化，工程上更简单。这种字节级、可用任意进制的形式，后来成为开源项目规避 IBM 算术编码专利的常用路径。

和经典算术编码差在哪

把第三节的重归一化从”看最高比特”改成”看最高字节”，就是 Range Coder 的核心：

• 维护 low（区间下界）和 range（宽度），不再显式存 high。
• 每当 range < 2^{24}（高字节已确定），就把 low 的高字节写出，range 左移 8 位——对应第三节情况一/二的字节版。
• 进位不确定时，字节先塞进 cache，等方向明确再刷——对应第三节的 pending / underflow。

Range Coder 实现

Range Coder 的核心思想是维护两个变量：range（区间宽度）和 low（区间下界）。

typedef struct {
    uint32_t range;
    uint32_t low;
    uint32_t cache;
    uint32_t cache_size;
    uint8_t *output;
    size_t   out_pos;
} RangeEncoder;

void rc_init_encoder(RangeEncoder *rc, uint8_t *output) {
    rc->range      = 0xFFFFFFFF;
    rc->low        = 0;
    rc->cache      = 0;
    rc->cache_size = 1;
    rc->output     = output;
    rc->out_pos    = 0;
}

static void rc_shift_low(RangeEncoder *rc) {
    // 检查是否需要进位
    if ((uint32_t)rc->low < 0xFF000000 || (rc->low >> 32)) {
        uint32_t carry = (uint32_t)(rc->low >> 32);
        rc->output[rc->out_pos++] = (uint8_t)(rc->cache + carry);
        while (rc->cache_size > 1) {
            rc->output[rc->out_pos++] = (uint8_t)(carry + 0xFF);
            rc->cache_size--;
        }
        rc->cache = (uint32_t)((rc->low >> 24) & 0xFF);
    } else {
        rc->cache_size++;
    }
    rc->low = (uint32_t)(rc->low << 8);
}

void rc_encode(RangeEncoder *rc, uint32_t cum_freq, uint32_t freq, uint32_t total) {
    rc->range /= total;
    rc->low   += cum_freq * rc->range;
    rc->range *= freq;

    // 字节级重归一化
    while (rc->range < (1 << 24)) {
        rc->range <<= 8;
        rc_shift_low(rc);
    }
}

Range Coder 的进位处理机制比算术编码的 underflow 处理更简洁。当 low 可能产生进位时（即高位不确定），我们将字节暂存在 cache 中，等到确定不会进位时再一起输出。

工程取舍

Range Coder 相比经典算术编码的差异可以归纳为：

字节级操作带来的速度提升是实实在在的，而字节对齐造成的压缩率损失在实践中通常可以忽略。这正是 7-Zip 的 LZMA 选择 Range Coder 的原因之一。

五、ANS：Jarek Duda 的革命性突破

2009 年，波兰数学家 Jarek Duda 在 arXiv 上发表了 Asymmetric Numeral Systems（非对称数字系统，ANS）。这项工作几乎改写了现代熵编码的工程版图。

核心洞察：一个整数代替一个区间

算术编码维护区间 \([\text{low}, \text{high})\)——两个端点、一套 pending/进位逻辑。ANS 只维护一个整数状态 \(x\)：编码是 \(x \mapsto x'\)，解码是逆映射，同时吐出符号。

最简单的 uANS（均匀 ANS）就是普通 \(M\) 进制：

\[ \text{编码}:\ x' = x\cdot M + s,\qquad \text{解码}:\ s = x \bmod M,\ \ x = \lfloor x/M\rfloor \]

问题立刻变成：非均匀分布怎么办？

从均匀到非均匀：用槽位密度表达概率

Duda 的做法是：把模 \(M\) 的余数轴切成 \(M\) 个槽位（slot），按频率 \(f_s\) 把槽位分给各符号。频率高的符号占更多槽，编码时状态涨得更慢——代价更接近 \(-\log_2 P(s)\)。

设 \(f_s\) 为频率、\(c_s\) 为累积频率、总频率 \(M=\sum f_s\)：

\[ \text{编码}:\quad x' = \left\lfloor \frac{x}{f_s}\right\rfloor\cdot M + (x \bmod f_s) + c_s \]

\[ \text{解码}:\quad \text{slot} = x \bmod M,\quad \text{取 } s \text{ 使 } c_s \le \text{slot} < c_s + f_s,\quad x' = f_s\left\lfloor \frac{x}{M}\right\rfloor + \text{slot} - c_s \]

编码与解码互为逆：从 \(x'\) 一定能还原 \((x, s)\)。上图 \(M=8\)、\(f(A)=4\)、\(f(B)=3\)、\(f(C)=1\) 时，编 A 状态约乘 2，编 C 约乘 8——稀有符号更贵，与信息论一致。

信息论的优美

为什么 ANS 能逼近熵？编码频率为 \(f_s\) 的符号时，状态近似按 \(x' \approx x\cdot M/f_s\) 增长：

\[ \log_2 x' \approx \log_2 x + \log_2\frac{M}{f_s} = \log_2 x + \log_2\frac{1}{P(s)} \]

\(\log_2 x\) 增加的量恰好是 \(-\log_2 P(s)\) 比特。算术编码用区间宽度的收缩来记账；ANS 用整数状态的增长来记账——两者在信息论上等价，状态表示不同。

六、rANS 与 tANS：两种实现路径

ANS 有两种主要的实现变体：rANS（range ANS）和 tANS（table ANS）。

rANS：基于算术运算

rANS 直接使用上一节的数学公式进行编码和解码。它的特点是：

• 使用除法和取模运算
• 不需要大型查找表
• 适合频率表动态变化的场景
• 编码器和解码器使用相同的频率表

rANS 的状态空间管理（重归一化）需要把状态 \(x\) 始终约束在 \([L,\ b\cdot L)\) 内。其中 \(L\) 是下界（32 位实现常取 \(L=2^{23}\)），\(b\) 是输出基数（常取 \(b=2^8=256\)，即按字节输出）。编码符号 \(s\) 之前，先按下式压低状态再执行 ANS 编码公式：

\[ x_{\max} = \left\lfloor \frac{b\cdot L}{M}\right\rfloor f_s - 1 \]

while (x > x_max):
    输出 x 的低 8 位
    x >>= 8

tANS：基于查找表

tANS 将整个 ANS 编码和解码过程预计算为一张查找表（有限状态机）。它的特点是：

• 没有除法操作（用表查找替代）
• 编码和解码都是 O(1) 的表查找 + 位移操作
• 需要预计算查找表
• 查找表的大小为 L 项（通常 L = 1024 到 4096）
• 不适合频率表频繁变化的场景

tANS 的查找表构建过程：

1. 选择表大小 L（通常 L = 2^R，R 在 10-12 之间）
2. 将 [0, L) 的位置按照符号频率分配给各符号
3. 对每个状态 x in [L, 2L):
   - 确定 x 对应的符号 s
   - 计算编码后的新状态和需要输出的比特数
   - 存入查找表

rANS vs tANS 对比

实际工程中的选择：

• LZFSE（Apple）：使用 tANS，因为 LZ 解析后的符号分布相对固定
• zstd（Facebook/Meta）：使用 tANS（称为 FSE），为每个压缩块预计算查找表
• JPEG XL：使用 rANS，因为需要在图像块之间快速切换上下文

七、Streaming rANS：逆向编码，正向解码

ANS 有一个反直觉的硬约束：编码必须逆序，解码才是正序。根因是栈（LIFO）语义，不是实现疏忽。

为什么需要逆向编码

编码相当于把符号 push 进状态 \(x\)，解码相当于 pop。后进先出：若希望解码器按 \(s_1,s_2,s_3\) 吐出，编码器就必须按 \(s_3,s_2,s_1\) 压入。

若正向 push s1→s2→s3，则 pop 得到 s3→s2→s1（颠倒）
要让 pop 得到 s1→s2→s3，必须先 push s3，再 s2，再 s1

算术编码是 FIFO：编解码同向。ANS 用单个整数换来了更简单的状态与更好的交错并行，代价就是这块必须缓冲、逆序扫一遍。

上图把重归一化也画进去了：编码侧状态过大就吐低字节；解码侧状态过小就从比特流补字节。整条路径是镜像的。

实际实现策略

在实际实现中，有几种策略来处理这个问题：

策略一：缓冲整个块

最简单的方法是将输入数据按块处理，将每个块的符号读入缓冲区，然后逆序编码：

// 将符号缓冲到数组中
for (i = 0; i < block_size; i++) {
    symbols[i] = read_next_symbol();
}

// 逆序编码
for (i = block_size - 1; i >= 0; i--) {
    rans_encode(&state, symbols[i], &output);
}

策略二：交错多个 rANS 流

这是一种高性能技巧，使用多个独立的 rANS 状态并行处理：

// 4 路交错 rANS
uint32_t state[4];
for (i = block_size - 4; i >= 0; i -= 4) {
    rans_encode(&state[3], symbols[i+3], &output);
    rans_encode(&state[2], symbols[i+2], &output);
    rans_encode(&state[1], symbols[i+1], &output);
    rans_encode(&state[0], symbols[i+0], &output);
}

// 写入 4 个最终状态
write_state(state[0]); write_state(state[1]);
write_state(state[2]); write_state(state[3]);

交错的好处是：

• 指令级并行：4 个独立的编码操作没有数据依赖
• 解码端同样可以并行：4 个解码器独立工作
• 隐藏延迟：除法指令的延迟被其他流的操作覆盖

这种技巧在 SIMD 优化中尤为重要：由于单条除法指令延迟较高，多路交错让 CPU 在等待一路除法结果时去推进其它路，能显著提升吞吐量（具体倍数取决于微架构与实现）。

逆向编码的缓冲区管理

编码器生成的比特流也是逆序的，这需要特殊的缓冲区管理：

typedef struct {
    uint8_t *buf;        // 缓冲区起始
    uint8_t *buf_end;    // 缓冲区末尾
    uint8_t *ptr;        // 当前写入位置（从末尾向前）
} RansOutBuf;

void rans_out_init(RansOutBuf *ob, uint8_t *buf, size_t size) {
    ob->buf     = buf;
    ob->buf_end = buf + size;
    ob->ptr     = buf + size;  // 从末尾开始写
}

void rans_out_byte(RansOutBuf *ob, uint8_t byte) {
    *(--ob->ptr) = byte;  // 向前写入
}

八、完整 rANS 编解码器实现

下面是一个完整的、可编译运行的 rANS 编解码器 C 实现。

/* rans.c -- 完整的 rANS 编解码器
 *
 * 编译: gcc -O2 -o rans rans.c
 * 运行: ./rans
 *
 * 参数说明:
 *   RANS_L     = 状态下界 (2^23)
 *   RANS_SCALE = 频率精度 (2^12 = 4096)
 *   RANS_BYT   = 输出基数 (2^8 = 256，字节对齐)
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <assert.h>
#include <math.h>
#include <time.h>

/* ------------------------------------------------------------------ */
/*  基本参数                                                           */
/* ------------------------------------------------------------------ */

#define RANS_L          (1u << 23)          /* 状态下界 L              */
#define RANS_SCALE_BITS 12                  /* 频率表精度比特数         */
#define RANS_SCALE      (1u << RANS_SCALE_BITS) /* M = 4096           */

typedef uint32_t RansState;

/* ------------------------------------------------------------------ */
/*  符号统计结构                                                       */
/* ------------------------------------------------------------------ */

typedef struct {
    uint16_t freq;       /* 量化后的频率 f_s               */
    uint16_t cum_freq;   /* 累积频率 c_s                   */
} RansSymbol;

typedef struct {
    RansSymbol sym[256];         /* 每个字节值的频率信息   */
    uint8_t    cum2sym[RANS_SCALE]; /* 累积频率 -> 符号 查找表 */
} RansStats;

/* ------------------------------------------------------------------ */
/*  构建频率表                                                         */
/* ------------------------------------------------------------------ */

static void stats_build(RansStats *st, const uint8_t *data, size_t len) {
    uint32_t counts[256] = {0};
    size_t i;

    for (i = 0; i < len; i++)
        counts[data[i]]++;

    /* 将原始频率量化到 [0, RANS_SCALE] 范围 */
    /* 确保每个出现过的符号至少有频率 1         */
    uint32_t total = 0;
    uint16_t freqs[256] = {0};

    for (i = 0; i < 256; i++) {
        if (counts[i] > 0) {
            freqs[i] = (uint16_t)((uint64_t)counts[i] * RANS_SCALE / len);
            if (freqs[i] == 0)
                freqs[i] = 1;
            total += freqs[i];
        }
    }

    /* 调整总频率使其精确等于 RANS_SCALE */
    /* 找到频率最大的符号来吸收误差       */
    int32_t diff = (int32_t)RANS_SCALE - (int32_t)total;
    if (diff != 0) {
        uint32_t max_freq = 0;
        int max_idx = 0;
        for (i = 0; i < 256; i++) {
            if (freqs[i] > max_freq) {
                max_freq = freqs[i];
                max_idx = (int)i;
            }
        }
        freqs[max_idx] = (uint16_t)((int32_t)freqs[max_idx] + diff);
    }

    /* 构建累积频率表和逆查找表 */
    uint16_t cum = 0;
    for (i = 0; i < 256; i++) {
        st->sym[i].freq     = freqs[i];
        st->sym[i].cum_freq = cum;
        /* 填充逆查找表 */
        uint16_t j;
        for (j = 0; j < freqs[i]; j++)
            st->cum2sym[cum + j] = (uint8_t)i;
        cum += freqs[i];
    }
    assert(cum == RANS_SCALE);
}

/* ------------------------------------------------------------------ */
/*  rANS 编码器                                                        */
/* ------------------------------------------------------------------ */

static inline void rans_enc_init(RansState *state) {
    *state = RANS_L;
}

static inline void rans_enc_flush(RansState state, uint8_t **pptr) {
    /* 将最终状态以 4 字节写入 */
    uint8_t *ptr = *pptr;
    ptr -= 4;
    ptr[0] = (uint8_t)(state >>  0);
    ptr[1] = (uint8_t)(state >>  8);
    ptr[2] = (uint8_t)(state >> 16);
    ptr[3] = (uint8_t)(state >> 24);
    *pptr = ptr;
}

static inline void rans_enc_put(RansState *state, uint8_t **pptr,
                                uint16_t freq, uint16_t cum_freq) {
    RansState x = *state;

    /* 重归一化: 确保编码后状态不会溢出 */
    /* x_max = floor((RANS_L >> RANS_SCALE_BITS) << 8) * freq - 1 */
    uint32_t x_max = ((RANS_L >> RANS_SCALE_BITS) << 8) * freq;
    while (x >= x_max) {
        *(--(*pptr)) = (uint8_t)(x & 0xFF);
        x >>= 8;
    }

    /* ANS 编码公式: x' = (x / f) * M + (x % f) + c */
    *state = ((x / freq) << RANS_SCALE_BITS) + (x % freq) + cum_freq;
}

/* ------------------------------------------------------------------ */
/*  rANS 解码器                                                        */
/* ------------------------------------------------------------------ */

static inline void rans_dec_init(RansState *state, const uint8_t **pptr) {
    const uint8_t *ptr = *pptr;
    uint32_t x;
    x  = (uint32_t)ptr[0] <<  0;
    x |= (uint32_t)ptr[1] <<  8;
    x |= (uint32_t)ptr[2] << 16;
    x |= (uint32_t)ptr[3] << 24;
    *pptr += 4;
    *state = x;
}

static inline uint32_t rans_dec_get(RansState state) {
    return state & (RANS_SCALE - 1);
}

static inline void rans_dec_advance(RansState *state, const uint8_t **pptr,
                                    uint16_t freq, uint16_t cum_freq) {
    RansState x = *state;

    /* ANS 解码状态更新: x' = f * (x / M) + (x % M) - c */
    x = freq * (x >> RANS_SCALE_BITS) + (x & (RANS_SCALE - 1)) - cum_freq;

    /* 重归一化: 从比特流中读入字节 */
    while (x < RANS_L) {
        x = (x << 8) | **pptr;
        (*pptr)++;
    }

    *state = x;
}

/* ------------------------------------------------------------------ */
/*  编码入口                                                           */
/* ------------------------------------------------------------------ */

static size_t rans_compress(const uint8_t *input, size_t in_len,
                            uint8_t *output, size_t out_cap,
                            const RansStats *st) {
    uint8_t *ptr = output + out_cap;   /* 从尾部向前写 */
    RansState state;
    size_t i;

    rans_enc_init(&state);

    /* 关键: 逆序编码 */
    for (i = in_len; i > 0; i--) {
        uint8_t s = input[i - 1];
        rans_enc_put(&state, &ptr, st->sym[s].freq, st->sym[s].cum_freq);
    }

    rans_enc_flush(state, &ptr);

    /* 将压缩数据移动到 output 开头 */
    size_t compressed_size = (size_t)(output + out_cap - ptr);
    memmove(output, ptr, compressed_size);
    return compressed_size;
}

/* ------------------------------------------------------------------ */
/*  解码入口                                                           */
/* ------------------------------------------------------------------ */

static void rans_decompress(const uint8_t *compressed, size_t comp_len,
                            uint8_t *output, size_t out_len,
                            const RansStats *st) {
    const uint8_t *ptr = compressed;
    RansState state;
    size_t i;
    (void)comp_len;

    rans_dec_init(&state, &ptr);

    /* 正序解码 */
    for (i = 0; i < out_len; i++) {
        uint32_t slot = rans_dec_get(state);
        uint8_t  s    = st->cum2sym[slot];
        output[i] = s;
        rans_dec_advance(&state, &ptr, st->sym[s].freq, st->sym[s].cum_freq);
    }
}

/* ------------------------------------------------------------------ */
/*  测试与基准                                                         */
/* ------------------------------------------------------------------ */

static double get_time_ms(void) {
    struct timespec ts;
    clock_gettime(CLOCK_MONOTONIC, &ts);
    return ts.tv_sec * 1000.0 + ts.tv_nsec / 1e6;
}

static double calc_entropy(const uint8_t *data, size_t len) {
    uint32_t counts[256] = {0};
    double entropy = 0.0;
    size_t i;
    for (i = 0; i < len; i++) counts[data[i]]++;
    for (i = 0; i < 256; i++) {
        if (counts[i] == 0) continue;
        double p = (double)counts[i] / (double)len;
        entropy -= p * (log(p) / log(2.0));
    }
    return entropy;
}

int main(void) {
    /* 生成测试数据: 模拟非均匀分布 */
    size_t data_size = 1 << 20;  /* 1 MB */
    uint8_t *original   = (uint8_t *)malloc(data_size);
    uint8_t *compressed = (uint8_t *)malloc(data_size * 2);
    uint8_t *decoded    = (uint8_t *)malloc(data_size);

    srand(42);
    size_t i;
    for (i = 0; i < data_size; i++) {
        /* 偏斜分布: 大部分值集中在 0-15 */
        int r = rand() % 100;
        if (r < 60)      original[i] = (uint8_t)(rand() % 4);
        else if (r < 85)  original[i] = (uint8_t)(4 + rand() % 12);
        else if (r < 95)  original[i] = (uint8_t)(16 + rand() % 48);
        else               original[i] = (uint8_t)(64 + rand() % 192);
    }

    double entropy = calc_entropy(original, data_size);
    printf("=== rANS 编解码器测试 ===\n\n");
    printf("原始数据大小 : %zu bytes\n", data_size);
    printf("零阶信息熵   : %.4f bits/byte\n", entropy);
    printf("理论最小大小 : %.0f bytes\n\n", entropy * data_size / 8.0);

    /* 构建频率表 */
    RansStats stats;
    stats_build(&stats, original, data_size);

    /* 编码 */
    double t0 = get_time_ms();
    size_t comp_size = rans_compress(original, data_size,
                                     compressed, data_size * 2, &stats);
    double t1 = get_time_ms();
    double enc_time = t1 - t0;

    printf("压缩后大小   : %zu bytes (%.2f bits/byte)\n",
           comp_size, 8.0 * comp_size / data_size);
    printf("压缩率       : %.2f%%\n", 100.0 * comp_size / data_size);
    printf("编码时间     : %.2f ms (%.1f MB/s)\n",
           enc_time, data_size / (enc_time * 1000.0));

    /* 解码 */
    double t2 = get_time_ms();
    rans_decompress(compressed, comp_size, decoded, data_size, &stats);
    double t3 = get_time_ms();
    double dec_time = t3 - t2;

    printf("解码时间     : %.2f ms (%.1f MB/s)\n",
           dec_time, data_size / (dec_time * 1000.0));

    /* 验证 */
    if (memcmp(original, decoded, data_size) == 0) {
        printf("\n验证通过: 解码数据与原始数据完全一致\n");
    } else {
        printf("\n错误: 解码数据与原始数据不一致!\n");
        free(original); free(compressed); free(decoded);
        return 1;
    }

    free(original);
    free(compressed);
    free(decoded);
    return 0;
}

代码要点解析

频率表量化：将真实频率量化到 RANS_SCALE（4096）范围。最关键的约束是所有频率之和必须精确等于 M，否则编解码器会失步。代码中通过将误差分配给频率最大的符号来保证这一点。

逆向编码：编码循环从 in_len - 1 到 0，同时输出缓冲区也从末尾向前写入。这是 rANS 的核心特征。

重归一化条件：编码时检查 x >= x_max，当状态过大时输出低字节并右移；解码时检查 x < RANS_L，当状态过小时读入字节并左移。这保证了状态始终在 [L, b*L) 范围内。

逆查找表：cum2sym 表将累积频率直接映射到符号，使得解码只需要一次数组访问就能确定符号，而不需要二分查找。

九、基准测试：Huffman vs 算术编码 vs rANS

下面用同一套实现实测比较 Huffman、Range Coder 与 rANS。三个编解码器共用同一个零阶静态模型（频率量化到 \(M=4096\)，与第八节一致），都是单路（非交错）、未做 SIMD / 表优化的参考实现，全部源码见本文同目录的 bench.c，可直接 gcc -O2 -o bench bench.c -lm && ./bench 复现。

测试方法

数据集为可复现的合成分布：uniform-256（256 符号均匀）、english-synth（按英文字母频率合成）、image-residual（零附近的几何分布，模拟图像残差）、binary-p0.95 / binary-p0.99（两符号偏斜源）。

压缩率对比

实测压缩率（bits per input byte，越接近 entropy 越好）：

几个观察：

1. 偏斜二元源是 Huffman 的 “阿喀琉斯之踵”：无论 \(p=0.95\) 还是 \(0.99\)，Huffman 都被钉死在 1.0 bit/符号，而 rANS 与 Range Coder 直接贴住熵（0.2866、0.0809），相差一个数量级。这正是前面”1 比特粒度”问题的实测体现。

2. 接近均匀、字母表又大时，Huffman 反而略胜：uniform-256 上 Huffman 给出 8.0000，而 rANS / Range 是 8.0318——把 256 个几乎相等的频率量化到 \(M=4096\) 反而引入了量化误差，加上每个 rANS 流末尾要刷 4 字节状态。这说明 ANS 并非在所有分布上都占优。

3. 一般偏斜分布（english、image-residual）：rANS 与 Range Coder 几乎贴住熵（差 0.001 以内），Huffman 稳定高出一截。rANS 与 Range Coder 在压缩率上基本无差别，因为它们用的是同一个量化模型。

速度对比

在 image-residual（16 MiB，7 次中位数）上的吞吐量：

这些数字只代表上述朴素参考实现，不能用来给三种算法的”最快速度”排座次，几个要点：

1. rANS 编解码基本对称（203 / 206 MB/s），因为编解码路径几乎是镜像的乘除法。
2. Range Coder 解码快于编码：本实现编码端的进位传播比解码端多做事。
3. 本文 Huffman 解码很慢（102 MB/s）是实现问题：这里用的是逐比特走规范码的朴素解码；换成查表式 Huffman 解码会快出数倍。所以不要据此得出”Huffman 解码慢”的结论。
4. 真实高性能实现远快于此：生产级 rANS 普遍采用 4–8 路交错 + SIMD 来掩盖除法延迟，tANS（FSE）则把整条编码路径预计算成查找表、完全消除除法，解码尤其快——这些都不是本文单路参考实现所能代表的。

定性总结

把压缩率（实测）与工程属性（推理与文献）合在一起：

十、专利历史：算术编码的枷锁与 ANS 的自由

专利问题深刻地影响了熵编码技术的采用和发展。这段历史值得详细回顾。

算术编码的专利困局

专利的影响是深远的。JPEG 标准本身就是一个活生生的例子：算术编码的压缩率通常比 Huffman 高约 5%–10%，但由于专利顾虑，几乎所有 JPEG 编码器都只用 Huffman。即使专利早已过期，绝大多数 JPEG 实现至今仍只支持 Huffman——生态一旦固化就很难回头。

ANS 的开放之路

Jarek Duda 从一开始就坚持 ANS 应该是自由的、不受专利限制的技术。

Duda 的立场非常明确：ANS 是数学，数学不应该被专利化。他公开表示，若有公司试图在 ANS 上设置专利壁垒，他将提供充分的现有技术来挑战这些专利。

这个故事的教训是：开放的技术标准和专利自由的实现，对于推动技术普及至关重要。ANS 之所以能在短短十年内被 Apple、Facebook/Meta、Google（JPEG XL）等巨头广泛采用，与其无专利负担密切相关。

十一、工程实践中的陷阱与细节

在实现和使用 ANS/算术编码时，有许多看似微小但足以导致灾难的工程细节。以下是经过实战检验的经验总结。

关键陷阱速查表

频率表量化的艺术

频率表量化是 ANS 实现中最容易出错的环节。核心要求：

1. 所有频率之和精确等于 M
2. 每个出现过的符号频率至少为 1
3. 频率比例尽可能接近真实概率

/* 频率量化: 三轮策略 */
void quantize_freqs(uint32_t *raw, uint16_t *qf,
                    int nsym, uint32_t total, uint32_t M) {
    uint32_t assigned = 0;
    int i;

    /* 第一轮: 非零符号至少分配 1 */
    int nz = 0;
    for (i = 0; i < nsym; i++) {
        if (raw[i] > 0) { qf[i] = 1; assigned++; nz++; }
        else             { qf[i] = 0; }
    }

    /* 第二轮: 按比例分配剩余预算 */
    uint32_t rem = M - nz;
    for (i = 0; i < nsym; i++) {
        if (raw[i] == 0) continue;
        uint32_t extra = (uint32_t)((uint64_t)raw[i] * rem / total);
        qf[i] += (uint16_t)extra;
        assigned += extra;
    }

    /* 第三轮: 贪心微调使总和精确等于 M */
    /* 每次调整选择使 KL 散度变化最小的符号 */
    while (assigned != M) {
        int dir = (assigned < M) ? 1 : -1;
        double best = (dir > 0) ? -1e30 : 1e30;
        int idx = -1;
        for (i = 0; i < nsym; i++) {
            if (raw[i] == 0 || (dir < 0 && qf[i] <= 1)) continue;
            double p = (double)raw[i] / total;
            double q0 = (double)qf[i] / M;
            double q1 = (double)(qf[i] + dir) / M;
            double v = p * log(q1 / q0);
            if ((dir > 0 && v > best) || (dir < 0 && v < best))
                { best = v; idx = i; }
        }
        qf[idx] += dir;
        assigned += dir;
    }
}

状态空间的选择

rANS 实现中，状态空间参数的选择影响着压缩率和速度的平衡：

常见配置：SCALE_BITS=12, L=2^23 适合大多数场景；SCALE_BITS=14, L=2^31 追求压缩率；SCALE_BITS=10, L=2^23 适合小型嵌入式系统。

十二、现实世界中的 ANS 应用

ANS 已经被广泛应用于现代压缩系统中。让我们看看几个重要的实际案例。

JPEG XL (2020-)

JPEG XL 是新一代图像压缩标准，由 Google 和 Cloudinary 等联合推动。它使用 rANS（而非 tANS）作为核心熵编码引擎，并把多个概率上下文聚类以减少频率表的传输开销。

之所以选 rANS 而非 tANS，核心原因在于上下文切换的代价：图像编码需要在大量不同的概率上下文之间频繁切换，rANS 切换上下文只需更换频率表指针，而 tANS 每次都要重建整个查找表，代价过高。此外，ANS 无专利负担这一点对标准化也至关重要。

LZFSE (Apple, 2016)

LZFSE 是 Apple 在 2016 年开源的压缩算法，用于 iOS 和 macOS。它先做 LZ77 解析，再用多个独立的 tANS（FSE）流分别编码 literal 值、literal 长度、match 长度和 match 距离这几类符号，输出按字节对齐。

LZFSE 的设计目标是”接近 zlib 的压缩率，但在 Apple 自家硬件上更快、更省电”，因此针对 ARM64 做了深度优化。它在通用文本/二进制数据上通常能取得优于 zlib 的压缩率与速度，但对极小数据或已压缩数据并无优势。

zstd / FSE (Facebook/Meta, 2015-)

zstd 是目前最成功的通用压缩算法之一。它的熵编码引擎 FSE（Finite State Entropy）是 tANS 的一种实现，同时 zstd 也用一个独立的 Huffman 编码器处理 literals。

zstd 的大致流程：数据分块后逐块独立处理，LZ77 解析生成 literals、match length、match offset 等符号流，其中 literals 用 Huffman、长度与偏移类符号用 FSE 编码，频率表写入块头。FSE 的表大小可配置，并支持预定义频率表以省去传表开销——这对小块数据尤其重要。相比纯 Huffman，FSE 让 zstd 在符号分布偏斜时拿到更接近熵的压缩率，而解码仍是纯查表，速度与 Huffman 相当。

其他应用

十三、ANS 与算术编码的深度对比

让我们从更本质的层面比较这两种方法。

信息论视角

两者在信息论意义上等价：都能逼近信息熵极限。但它们的”编码路径”完全不同：

操作复杂度

为什么 ANS 更快

ANS 的速度优势来自几个方面：

1. 单一状态变量：只有一个整数 x 需要跟踪，而算术编码需要 low 和 high 两个变量。

2. 无分支重归一化：rANS 的重归一化可以写成无分支代码（对于已知的频率表精度）。

3. 多路交错：由于每个 rANS 状态是独立的，可以轻松实现 4 路甚至 8 路交错编解码。

4. tANS 消除除法：tANS 将所有运算预计算到查找表中，完全消除了除法操作。

/* tANS 解码: 极致简洁，零除法 */
typedef struct {
    uint16_t new_x;      /* 基础新状态                    */
    uint8_t  symbol;     /* 输出符号                      */
    uint8_t  nb_bits;    /* 需要从比特流中读取的比特数     */
} TansDecEntry;

static inline uint8_t tans_decode(uint32_t *state,
                                  const TansDecEntry *table,
                                  const uint8_t **bitstream) {
    TansDecEntry e = table[*state];
    /* 读取 nb_bits 个比特并更新状态 */
    *state = e.new_x;
    /* 从比特流读入比特（此处简化） */
    *state = (*state << e.nb_bits) | read_bits(bitstream, e.nb_bits);
    return e.symbol;
}

什么时候仍然选择算术编码

尽管 ANS 在很多方面优于算术编码，但算术编码在某些场景下仍有优势：

1. 自适应编码：概率模型频繁更新时，算术编码不需要预计算任何表。
2. 二进制算术编码：当符号只有两个值时（如 CABAC），专门优化的二进制算术编码器非常高效。
3. 硬件实现：算术编码器的流水线结构在某些硬件加速器中更自然。
4. 正向编码：不需要缓冲数据来逆向编码。

十四、个人思考

写这篇文章的过程中，我反复感叹信息论的深度和熵编码技术发展的曲折。

关于 Huffman 的 “够用论”。在很多工程场景中，Huffman 编码确实 “够用” 了。当符号概率接近 2 的负整数次幂时，Huffman 几乎是最优的，而且实现简单、速度快。不要因为学了 ANS 就觉得到处都该用 ANS。选择编码方案应该基于实际的概率分布和性能需求。

关于 ANS 的美。从纯数学的角度看，ANS 比算术编码更优雅。算术编码需要维护一个区间的两个端点，而 ANS 只需要一个整数。这种简约不仅仅是审美上的，它直接导致了更好的并行性和更高的吞吐量。Duda 的工作展示了一个深刻的道理：数学上更简洁的形式往往也导致更高效的实现。

关于专利与创新。算术编码的专利史是一个反面教材。IBM 的专利可能短期内保护了其商业利益，但长期来看，它阻碍了更好的压缩技术的普及，最终导致整个行业在 JPEG 基线中 “锁定” 在 Huffman 编码上长达三十年。相比之下，Duda 坚持 ANS 的开放性，使其在十年内就被多个行业标准和主流软件采用。这个对比值得每一个技术工作者深思。

关于工程与理论的距离。从理论上说，算术编码在 1976 年就已经 “解决” 了信息编码的最优性问题。但从工程上说，直到 ANS 出现（2009 年），再到 zstd（2016 年）和 JPEG XL（2020 年）的普及，这条路走了将近半个世纪。信息论的完美理论和可部署的高性能实现之间，隔着的不只是代码量，还有专利壁垒、硬件限制、生态惯性，以及无数工程细节的积累。

关于未来。随着数据量的爆炸性增长，高效的熵编码变得越来越重要。ANS/FSE 目前是通用压缩的最佳选择，但我认为未来可能出现针对特定硬件（如 GPU、专用加速器）优化的新型熵编码方案。同时，机器学习驱动的概率建模正在快速发展，更好的概率模型与高效的 ANS 编码器的结合，可能带来新一轮的压缩率突破。

十五、参考资料

1. Huffman, D.A. “A Method for the Construction of Minimum-Redundancy Codes.” Proceedings of the IRE, 1952.
2. Rissanen, J.J. “Generalized Kraft inequality and arithmetic coding.” IBM Journal of R&D, 1976.
3. Witten, I.H., Neal, R.M., and Cleary, J.G. “Arithmetic Coding for Data Compression.” CACM, 1987.
4. Martin, G.N.N. “Range encoding: An algorithm for removing redundancy.” Video & Data Recording Conf, 1979.
5. Duda, J. “Asymmetric numeral systems.” arXiv:0902.0271, 2009.
6. Duda, J. “ANS: entropy coding combining speed of Huffman with compression rate of arithmetic coding.” arXiv:1311.2540, 2013.
7. Collet, Y. “Finite State Entropy - A new breed of entropy coder.” 2013.
8. Alakuijala, J. et al. “JPEG XL next-generation image compression architecture.” Proc. SPIE 11137, 2019.
9. Apple Inc. “LZFSE - LZ style compression using Finite State Entropy coding.” GitHub, 2016.
10. RFC 8478: Zstandard Compression and the application/zstd Media Type, 2018.

上一篇: zstd 深度解剖 下一篇: 整数压缩 相关阅读: - Huffman 编码与 DEFLATE - zstd 深度解剖