數據壓縮之ROLZ字典編碼

　　在字典編碼中，最常用的恐怕要算LZ77編碼了。LZ77的思想很簡單，就是用一個<offset, length>元組來表示當前位置的字節串在前offset個字節中出現過。正是由於這個簡單的思想，所有基於LZ77實現的實用算法都有着不錯的解壓速度。經典的使用LZ77編碼的壓縮算法有zip/gz的deflate算法，7z的lzma算法等。

　　在對LZ77算法研究中，我們也發現算法中的一些不足之處，LZ77最明顯的不足是offset值的過度零散導致對<offset, length>元組的后續處理效果不好。例如處理一個16MB的數據塊，一個<offset, length>中offset的取值就有16777216種之多，雖然可以對offset進行分段處理，但還是或多或少會影響到壓縮率的提升。

　　於是我們有了ROLZ算法，ROLZ全稱是Reduced Offset Lempel Ziv，即減少了offset的LZ編碼。在ROLZ中，重復串出現的位置不再是用相對當前的偏移量offset，而是用一個表項編號index來表示。

　　在詳細介紹ROLZ算法之前，我們先介紹一個壓縮算法中常用的概念——上下文（context），上下文就是當前編碼位置之前的內容，在實際應用中，我們只取當前編碼位置前k個字符做為上下文，稱為k階上下文。上下文是提高壓縮率的一個重要工具，例如字符串"queen"，如果用一般的統計模型來處理，由於英語中字母"u"出現頻率不高，我們可能會給"u"分配一個長的前綴編碼。但考慮一階上下文"q"，在英語中"q"后面出現"u"的概率非常高，我們就會給"u"分配相對較短的前綴編碼，從而利用上下文提高了編碼效率。

　　在編碼中，我們總會保存每個上下文的一些信息，但是當上下文階數k >= 3時，上下文的個數會急劇升高，當k=4時就已經有256^4 = 4GB的上下文各數。所以對於較長的上下文，我們往往取一個對前k個字符的一個Hash值作為上下文，這樣的結果是可能幾個完全不同的上下文會占用相同的存儲空間，導致上下文的預測准確率有所降低，但這種做法滿足了程序對內存的要求，實際中也被廣泛使用。

　　現在進入正題，在ROLZ算法中，我們會建立一張二維表rolz_table[context][index]，在處理完每個位置后，我們將這個位置存入對應上下文的桶中，以下就字符串"construct- destruct"舉個例子（考慮一階上下文，假設當前編碼到位置12）：

0 1 2 3 4 5 6 7 8 9 10 11 $ 13 14 15 16 17
c o n s t r u c t - d  e  s t  r  u  c  t

table["c"] = {8, 1}
table["o"] = {2}
table["n"] = {3}
table["s"] = {4}
table["t"] = {9, 5}
table["r"] = {6}
table["u"] = {7}
table["-"] = {10}
table["d"] = {11}

　　如果是傳統的LZ77算法，這時候我們應該輸出一個匹配對<offset=9, length=6>表示重復串"struct"。但在ROLZ算法中，我們只在對應上下文的桶中尋找匹配，在當前上下文"e"的桶中找不到匹配（因為是一個空的桶），只好將"s"原樣輸出，然后將當前位置加入上下文"e"的桶中：

 

0 1 2 3 4 5 6 7 8 9 10 11 12 $ 14 15 16 17
c o n s t r u c t - d  e  s  t r  u  c  t

table["c"] = {8, 1}
table["o"] = {2}
table["n"] = {3}
table["s"] = {4}
table["t"] = {9, 5}
table["r"] = {6}
table["u"] = {7}
table["-"] = {10}
table["d"] = {11}
table["e"] = {12}

　　現在編碼位置13的"t"，我們在上下文"s"的桶中尋找匹配，發現table["s"][0]就是我們要找的匹配，於是匹配成功，輸出一個ROLZ匹配項<index=0, length=5>表示當前位置的重復串"truct"。

　　所以，對字符串"construct-destruct"，LZ77和ROLZ的編碼結果分別如下：

LZ77: construct-de<9,6>
ROLZ: construct-des<0,5>

　　可以看到ROLZ的輸出多了一個字符"s"，但是offset變成了更容易處理的index，實踐表明對大量輸入的時候ROLZ通常能得到更好的壓縮率。

　　ROLZ的編碼實現比LZ77要簡單得多，因為rolz_table本身就是一張很好的散列表，我們不需要像LZ77那樣另外創建一個查找結構。而解壓相對LZ77要復雜一些，我們必須和壓縮過程一樣構建起整個rolz_table，然后通過查表將index轉換為重復串的位置，所以ROLZ解壓相對LZ77要慢一些。

　　我個人有一個ROLZ+算術編碼的完整實現comprolz，目前壓縮率已經和7z相當。下面也附上一個簡單的ROLZ+Polar編碼實現（壓縮率略高於gz），使用3階Hash上下文，每個上下文有15個桶項：

/*******************************************************************************
 * RichSelian's nooblike compressor: ROLZ + Polar coding
 ******************************************************************************/
#include <stdio.h>
#include <string.h>

/*******************************************************************************
 * POLAR Coder
 ******************************************************************************/
#define POLAR_SYMBOLS   512 /* should be even */
#define POLAR_MAXLEN    15  /* should be less than 16, so we can pack two length values into a byte */

#define M_round_down(x)     while((x)&(-(x)^(x))) { (x) &= (-(x)^(x)); }
#define M_round_up(x)       while((x)&(-(x)^(x))) { (x) &= (-(x)^(x)); } (x) <<= 1;
#define M_int_swap(x, y)    {int (_)=(x); (x)=(y); (y)=(_);}

int polar_make_leng_table(const int* freq_table, int* leng_table) {
    int symbols[POLAR_SYMBOLS];
    int i;
    int s;
    int total;
    int shift = 0;

    memcpy(leng_table, freq_table, POLAR_SYMBOLS * sizeof(int));

MakeTablePass:
    /* sort symbols */
    for(i = 0; i < POLAR_SYMBOLS; i++) {
        symbols[i] = i;
    }
    for(i = 0; i < POLAR_SYMBOLS; i++) {
        if(i > 0 && leng_table[symbols[i - 1]] < leng_table[symbols[i]]) {
            M_int_swap(symbols[i - 1], symbols[i]);
            i -= 2;
        }
    }

    /* calculate total frequency */
    total = 0;
    for(i = 0; i < POLAR_SYMBOLS; i++) {
        total += leng_table[i];
    }

    /* run */
    M_round_up(total);
    s = 0;
    for(i = 0; i < POLAR_SYMBOLS; i++) {
        M_round_down(leng_table[i]);
        s += leng_table[i];
    }
    while(s < total) {
        for(i = 0; i < POLAR_SYMBOLS; i++) {
            if(s + leng_table[symbols[i]] <= total) {
                s += leng_table[symbols[i]];
                leng_table[symbols[i]] *= 2;
            }
        }
    }

    /* get code length */
    for(i = 0; i < POLAR_SYMBOLS; i++) {
        s = 2;
        if(leng_table[i] > 0) {
            while((total / leng_table[i]) >> s != 0) {
                s += 1;
            }
            leng_table[i] = s - 1;
        } else {
            leng_table[i] = 0;
        }

        /* code length too long -- scale and rebuild table */
        if(leng_table[i] > POLAR_MAXLEN) {
            shift += 1;
            for(i = 0; i < POLAR_SYMBOLS; i++) {
                if((leng_table[i] = freq_table[i] >> shift) == 0 && freq_table[i] > 0) {
                    leng_table[i] = 1;
                }
            }
            goto MakeTablePass;
        }
    }
    return 0;
}

int polar_make_code_table(const int* leng_table, int* code_table) {
    int i;
    int s;
    int t1;
    int t2;
    int code = 0;

    memset(code_table, 0, POLAR_SYMBOLS * sizeof(int));

    /* make code for each symbol */
    for(s = 1; s <= POLAR_MAXLEN; s++) {
        for(i = 0; i < POLAR_SYMBOLS; i++) {
            if(leng_table[i] == s) {
                code_table[i] = code;
                code += 1;
            }
        }
        code *= 2;
    }

    /* reverse each code */
    for(i = 0; i < POLAR_SYMBOLS; i++) {
        t1 = 0;
        t2 = leng_table[i] - 1;
        while(t1 < t2) {
            code_table[i] ^= (1 & (code_table[i] >> t1)) << t2;
            code_table[i] ^= (1 & (code_table[i] >> t2)) << t1;
            code_table[i] ^= (1 & (code_table[i] >> t1)) << t2;
            t1++;
            t2--;
        }
    }
    return 0;
}

int polar_make_decode_table(const int* leng_table, const int* code_table, int* decode_table) {
    int i;
    int c;

    for(c = 0; c < POLAR_SYMBOLS; c++) {
        if(leng_table[c] > 0) {
            for(i = 0; i + code_table[c] < 65536; i += (1 << leng_table[c])) {
                decode_table[i + code_table[c]] = c;
            }
        }
    }
    return 0;
}

/*******************************************************************************
 * ROLZ
 ******************************************************************************/
#define ROLZ_BUCKET_SIZE    65536
#define MATCH_IDX_SIZE      15  /* make element of rolz_table[] 64 Bytes */
#define MATCH_LEN_MIN       2
#define MATCH_LEN_MAX       17  /* MATCH_LEN_MAX < MATCH_LEN_MIN + (POLAR_SYMBOLS-256) / MATCH_IDX_SIZE */

static int ch1 = 0;
static int ch2 = 0;
static int ch3 = 0;

#define M_rolz_item(x, n)  (rolz_table[(x)].m_item[(rolz_table[(x)].m_head + MATCH_IDX_SIZE - (n)) % MATCH_IDX_SIZE])

static struct {
    unsigned int  m_item[MATCH_IDX_SIZE];
    unsigned char m_head;
} rolz_table[ROLZ_BUCKET_SIZE];

static inline unsigned int rolz_context() {
    return (unsigned)(ch1 * 131313131 + ch2 * 131313 + ch3 * 131) % ROLZ_BUCKET_SIZE;
}

static inline void rolz_update_context(unsigned char* buf, int pos, int cache) {
    rolz_table[rolz_context()].m_head = (rolz_table[rolz_context()].m_head + 1) % MATCH_IDX_SIZE;
    if(cache) {
        M_rolz_item(rolz_context(), 0) = pos | (buf[pos] << 24);
    } else {
        M_rolz_item(rolz_context(), 0) = pos;
    }
    ch3 = ch2;
    ch2 = ch1;
    ch1 = buf[pos];
    return;
}

int rolz_encode(unsigned char* ibuf, unsigned short* obuf, int ilen) {
    int olen = 0;
    int pos = 0;
    int i;
    int j;
    int match_idx;
    int match_len;

    memset(rolz_table, 0, sizeof(rolz_table));
    ch3 = 0;
    ch2 = 0;
    ch1 = 0;
    while(pos < ilen) {
        match_len = MATCH_LEN_MIN - 1;
        match_idx = -1;
        if(pos + MATCH_LEN_MAX < ilen) { /* find match */
            for(i = 0; i < MATCH_IDX_SIZE; i++) {
                if(M_rolz_item(rolz_context(), i) == 0 || M_rolz_item(rolz_context(), i) >> 24 != ibuf[pos]) {
                    continue;
                }

                j = 1;
                while(j < MATCH_LEN_MAX && ibuf[pos + j] == ibuf[(M_rolz_item(rolz_context(), i) & 0x00ffffff) + j]) {
                    j++;
                }
                if(j > match_len) {
                    match_len = j;
                    match_idx = i;
                }
            }
        }
        if(match_len < MATCH_LEN_MIN) {
            match_len = 1;
            match_idx = -1;
        }

        if(match_idx == -1) { /* encode */
            obuf[olen++] = ibuf[pos];
        } else {
            obuf[olen++] = 256 + (match_len - MATCH_LEN_MIN) * MATCH_IDX_SIZE + match_idx;
        }

        for(i = 0; i < match_len; i++) { /* update context */
            rolz_update_context(ibuf, pos, 1);
            pos += 1;
        }
    }
    return olen;
}

int rolz_decode(unsigned short* ibuf, unsigned char* obuf, int ilen) {
    int olen = 0;
    int pos = 0;
    int i;
    int match_pos;
    int match_idx;
    int match_len;

    ch3 = 0;
    ch2 = 0;
    ch1 = 0;
    for(pos = 0; pos < ilen; pos++) {
        if(ibuf[pos] < 256) { /* decode */
            match_idx = -1;
            match_len = 1;
            obuf[olen++] = ibuf[pos];
        } else {
            match_idx = (ibuf[pos] - 256) % MATCH_IDX_SIZE;
            match_len = (ibuf[pos] - 256) / MATCH_IDX_SIZE + MATCH_LEN_MIN;
        }

        if(match_idx != -1) { /* expand match */
            match_pos = M_rolz_item(rolz_context(), match_idx);
            for(i =-0; i < match_len; i++) {
                obuf[olen++] = obuf[match_pos + i];
            }
        }
        for(i = 0; i < match_len; i++) { /* update context */
            rolz_update_context(obuf, olen - match_len + i, 0);
        }
    }
    return olen;
}

/*******************************************************************************
 * MAIN
 ******************************************************************************/
int main(int argc, char** argv) {
    static unsigned char  ibuf[10000000]; /* blocksize = 10MB */
    static unsigned short rbuf[10000000];
    static unsigned char  obuf[12000000];
    int ilen;
    int rlen;
    int olen;
    int rpos;
    int opos;
    int i;
    int freq_table[POLAR_SYMBOLS];
    int leng_table[POLAR_SYMBOLS];
    int code_table[POLAR_SYMBOLS];
    int decode_table[1 << (POLAR_MAXLEN + 1)];
    int code_buf;
    int code_len;

    if(argc == 2 && strcmp(argv[1], "e") == 0) {
        while((ilen = fread(ibuf, 1, sizeof(ibuf), stdin)) > 0) {
            rlen = rolz_encode(ibuf, rbuf, ilen);
            olen = 0;

            memset(freq_table, 0, sizeof(freq_table));
            code_buf = 0;
            code_len = 0;

            for(i = 0; i < rlen; i++) {
                freq_table[rbuf[i]] += 1;
            }
            polar_make_leng_table(freq_table, leng_table);
            polar_make_code_table(leng_table, code_table);

            /* write length table */
            for(i = 0; i < POLAR_SYMBOLS; i += 2) {
                obuf[olen++] = leng_table[i] * 16 + leng_table[i + 1];
            }

            /* encode */
            for(i = 0; i < rlen; i++) {
                code_buf += code_table[rbuf[i]] << code_len;
                code_len += leng_table[rbuf[i]];
                while(code_len > 8) {
                    obuf[olen++] = code_buf % 256;
                    code_buf /= 256;
                    code_len -= 8;
                }
            }
            if(code_len > 0) {
                obuf[olen++] = code_buf;
                code_buf = 0;
                code_len = 0;
            }
            fwrite(&rlen, sizeof(rlen), 1, stdout);
            fwrite(&olen, sizeof(olen), 1, stdout);
            fwrite(obuf, 1, olen, stdout);
        }
        return 0;
    }
    if(argc == 2 && strcmp(argv[1], "d") == 0) {
        while(fread(&rlen, sizeof(rlen), 1, stdin) == 1 && fread(&olen, sizeof(olen), 1, stdin) == 1) {
            olen = fread(obuf, 1, olen, stdin);
            rpos = 0;
            opos = 0;
            code_buf = 0;
            code_len = 0;

            /* read length table */
            for(i = 0; i < POLAR_SYMBOLS; i += 2) {
                leng_table[i] =     obuf[opos] / 16;
                leng_table[i + 1] = obuf[opos] % 16;
                opos++;
            }

            /* decode */
            polar_make_code_table(leng_table, code_table);
            polar_make_decode_table(leng_table, code_table, decode_table);

            while(rpos < rlen) {
                while(opos < olen && code_len < POLAR_MAXLEN) {
                    code_buf += obuf[opos++] << code_len;
                    code_len += 8;
                }
                i = decode_table[code_buf % 65536];

                rbuf[rpos++] = i;
                code_buf >>= leng_table[i];
                code_len -=  leng_table[i];
            }

            ilen = rolz_decode(rbuf, ibuf, rlen);
            fwrite(ibuf, 1, ilen, stdout);
        }
    }
    return -1;
}