Constant-weight code的压缩

问题

假设我们的信源都是constant-weight code,或者叫m-of-n code,也就是说它们有nn个bit,且其中有mm个为1。那么我们可以对信源进行压缩,最少只需log(nm)\log \binom{n}{m}个bit,即可表示所有的m-of-n code。那么,应该如何设计有效的编码和解码方法,能将m-of-n code压缩到log(nm)\log \binom{n}{m}个bit呢?StackExchange[1]上已经有人提出了该问题(另一种表述)。

简单方法

上面的链接中已经有人给出了一个回答。回答将所有的m-of-n code按字典序排好然后列出来,那么在编码的时候,只要把code对应的编号找到,将编号二进制表示即可满足要求:

ENCODE(S,n,k):1. If k=0, return 0.2. If nS, return ENCODE(S,n1,k).3. If nS, return (n1k)+ENCODE(S,n1,k1).\begin{aligned} &\operatorname{ENCODE}(S,n,k):\\ & \text{1. If }k = 0,\text{ return }0.\\ & \text{2. If }n \notin S,\text{ return }\operatorname{ENCODE}(S,n-1,k).\\ & \text{3. If }n \in S,\text{ return }\binom{n-1}{k}+\operatorname{ENCODE}(S,n-1,k-1). \end{aligned}

解码也很简单:

DECODE(x,n,k):1. If k=0, return .2. If x<(n1k), return DECODE(x,n1,k).3. If x(n1k), return DECODE(x(n1k),n1,k1){n}.\begin{aligned} &\operatorname{DECODE}(x,n,k):\\ & \text{1. If } k = 0,\text{ return }\emptyset.\\ & \text{2. If } x < \binom{n-1}{k},\text{ return }\operatorname{DECODE}(x,n-1,k).\\ & \text{3. If } x \geq \binom{n-1}{k},\text{ return }\operatorname{DECODE}(x - \binom{n-1}{k},n-1,k-1) \cup \{n\}. \end{aligned}

显然这种方法的平均和最坏时间都是O((nk))O(\binom{n}{k}),或者也需要O((nk))O(\binom{n}{k})的预处理以及存储空间。通常来说,这种指数级别的结果是无法接受的。

算术编码

算术编码非常适合此情形。同学发给我一篇thesis[2],是对算术编码的一个详细介绍,如有不了解算术编码的读者可以参考,或者也可以看看我写过的一个slides。我们把一个m-of-n code看成一个流式编码的过程。最开始的时候,字符1和0出现的概率分别是m/nm/n(nm)/n(n-m)/n。假设我们已经编码了aa个1和bb个0,即剩余mam-a个1和nmbn-m-b个0,那么下一个字符为1和0的概率分别是(ma)/(nab)(m-a)/(n-a-b)(nmb)/(nab)(n-m-b)/(n-a-b)。这样整个编码/解码流程的用时就缩减为O(n)O(n)
算术编码在此情形会遇到精度问题,因为理论上输入字符串会被编码成一个浮点数,由于计算机一般采用固定精度表示,nn较大时计算误差会一步步累计,最终甚至会导致解码失败。但如果我们分配O(n)O(n)空间来表示这个浮点数,其运算时间则将不会是常数。Ramabadran于1990年发表的文章[3](以下简称为文章R)就这个问题进行了讨论,并给出了一个解决方案。

A coding scheme for m-out-of-n codes

文章R本质上是基于Jones在1981年发布的文章[4](以下简称为文章J)。文章J研究的是如何将算术编码方法在固定数量寄存器,也就是我们的计算机上实现,但它假设了信源中各个字符出现的概率是完全独立的,也就是Huffman编码能处理的那种情形,且给出了算术编码在计算机上实现的一个算法。算法与标准的算术编码方法相比,在编码同样的字符串之后,得到的bit数会稍微多一点。记这个差值为ϵ\epsilon,文章J还给出了ϵ\epsilon关于寄存器个数、字符串长度及每个字符取值数量的一个上界。
文章R则把文章J的方法作了一些小的修改,运用到m-of-n code的情形,同样给出了算法和ϵ\epsilon的上界。上界的表达式比较复杂,简单的来说ϵ=O(nlogn/2k)\epsilon=O(n\log n/2^k),这里kk是寄存器的个数。
两篇文章的算法实现其实都有个BUG。先看文章J的TABLE II,有一行是tt,表示"digit ejected from x plus a possible carry",即进位值,最大值写的是2d12d-1。在0-1字符串中d=2d=2,所以最大值就是3。但是这篇文章在进制中表示数字是从1数起的,所以此处的3对应于文章R中应该为2。再看文章R的TABLE II,第一行的ee就是进位值,这里最大值写的是3。在我的实验中发现最大值确实是3,所以这里是文章J遗漏了一种情形,文章R发现了并且改了过来。然而这种遗漏的情形会导致算法出错,这一点文章R在修改算法的时候并没有发现。
文章R存储xx时,会存储其最低的w+1w+1位(以下简称为低位),但其余位(以下简称为高位)不能直接发送出去,因为低位在运算的时候可能会有进位,进位的可能值为0,1,2,3。为了解决这种情形,算法还会增加一个变量来存储高位中最后面有多少个连续的1。这样当低位进位时,比如说进位了1,那么算法就知道哪些连续的1都会因为进位而变成0,从而可以发送出去。但我们考虑高位的最后面是"111…101"的情形,这时算法不会记录这些连续的1,只会记录高位最后面是恰有一个连续的1,而前面这些连续的1都已经发送出去了。但如果低位运算后进位了3,那这些连续的1实际上都应该变成0,于是就导致了编码错误。为了解决这个问题,实际上我们可以再加一个变量来存储"01"前边有多少个1即可。在改动后我对大量数据进行了测试,算法都通过了。

Arithmetic Coding For Data Compression

之前提到的thesis[2:1]也给出了算术编码在计算机上实现的一个算法,算法来源于文章[5],而thesis的方法可以算的上是算术编码的一个主流实现了,因为它比上面提到的算法更加容易理解和实现。比较遗憾的是无论thesis还是文章[5:1]似乎都没有给出其算法与理论上算术编码差距,即ϵ\epsilon,的上界。
我们知道算术编码是一开始有一个[0,1]的区间,随着对字符的不断编码,这个区间会不断缩小,编码完成时我们用尽量少的bit输出一个区间内的数,用来代表这个区间。算法的整数实现,则是直接把区间[0,1]放大到[0,0x7fffffff],也就是考虑到现在机器上整型目前主要是32bit的情形。每当区间被分割时,分割点会取整,也就是我们用一个整数来近似替代分割点,于是得到一个新的整数区间。问题在于,随着区间的减小,分割点的近似会对精度产生越来越大的影响,比如到最后这个区间可能就是一个类似[x,x+1]这样的形式,那么算法将无法继续进行下去。为了避免区间过小,算法会设定一个阈值,每当区间长度小于这个阈值,区间就会被放大,也就是thesis里提到的E1,E2和E3 scaling。我没有进行细致的分析,仅从方法上的对比来看,它和前面提到的方法所带来的误差应该是一样的,即运用到m-of-n code上,应有ϵ=O(nlogn)/232\epsilon=O(n\log n)/2^{32}

另外两篇

前面StackExchange[1:1]中还给出了两个不同的方法[6] [7],但它们所针对的问题是如何将固定长度的字符串编码成一个m-of-n code。看起来我们只要互换编码函数和解码函数就能解决我们的问题,但事实并非如此。比如Knuth based文章[7:1]会给出一个映射

f: ABf:~A\rightarrow B

这里AA是固定长度字符串的集合,BB是m-of-n code集合。那么文章给出的ff一定是单射,又因为AABB的元素个数绝大多数情况不相同,因此ff的逆映射对于BB中的某些元素是无法定义的。如果我们互换编码函数和解码函数,编码函数就会对于某些m-of-n code无法编码,这样肯定是不行的。
对于Knuth based文章[7:2],我大致浏览了一下,感觉是无法适用于我们的问题。对于Geometric approach文章[6:1],我只知道它的时间复杂度是O(m2)O(m^2),是无关于nn的,因此只有m=o(n)m=o(\sqrt{n})时才会比算术编码更好。稍微读了下其方法,感觉实现起来会很麻烦,也不确定能不能改到适用我们的问题,就先搁置了。

一般化

我们把constant-weight code扩展为一个非负整数元素的向量,即m-of-n code表示向量长度为nn,元素之和为mm。这样的向量个数等于方程

x1+x2++xn=mx_1+x_2+\cdots+x_n=m

的非负整数解的个数,即(m+n1m)\binom{m+n-1}{m}。对于这个方程的解,我们按如下方式写一个0-1字符串:先写x1x_1个1,然后写一个0,再写x2x_2个1,再写一个0……最后写xnx_n个1,然后结束。这样一来,这个字符串长度必定为m+n1m+n-1,且有mm个1。用这种方法,我们很容易找到非负整数的m-of-n code与二元的m+n-1 -of- m code的一一对应,且方法用时为O(m+n)O(m+n)。再运用算术编码,我们就有了O(m+n)O(m+n)的压缩非负整数的m-of-n code的方法。
另外一个思路是修改算术编码的信源定义域,即将{0,1}改成{0,1,…,m},然后依然可以计算条件概率,最后的用时为O(nlogm)O(n\log m),不过条件概率涉及到复杂的组合数运算,精确计算会耗时多且会溢出,近似计算则精度可能会有较大影响,因而这种方法没有大的实际意义。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
# 文章1的Python实现
# The modifications of Elias scheme
# Application of arithmetic coding on m-out-of-n codes

import math
import itertools

# return ceil(log2(n))
def bits(n):
n-=1
counter=0
while n>0:
n>>=1
counter+=1
return counter

def enc(src):
n=len(src)
n0=n-sum(src)
w=20#bits(n)
x=0
y=v=1<<w
mask = v-1
rl=0
pb=True
dst=[]
for c in src:
z = (2*y*n0 + n) // (2*n)
n -= 1
if c == 1:
x += z
y -= z
else:
n0 -= 1
y = z
while y<v:
y<<=1
e=x>>w
x=(x&mask)<<1
if e>1:
if pb:
rl+=e-1
else:
dst.append(1)
dst+=[0]*(rl-1)
rl=e-2
else:
if pb:
dst.append(0)
dst+=[1]*rl
rl=e
else:
rl+=e

pb=((e&1)==0) or (rl==0)

e=x>>w
if e<2:
dst.append(0)
dst+=[1]*rl
if pb:
dst.append(0)
dst.append(e)
else:
if pb:
dst.append(0)
dst+=[1]*(rl+1)
else:
dst.append(1)
dst+=[0]*rl
dst.append(e-2)
return dst[3:]


def dec(dst, n, n0):
w = 20
dst += [1] * w
src=[]
x=0
y=1
v=1<<w
j=0
for i in range(n):
while y<v:
x=(x<<1)+dst[j]
y<<=1
j+=1
z=(2*y*n0 + n) // (2*n)
n-=1
if x<z:
src.append(0)
y=z
n0-=1
else:
src.append(1)
x-=z
y-=z

return src

def gen_code(d, c): # generate constant weight codes (c out of d)
for indices in itertools.combinations(range(d), c):
arr = [0] * d
for i in indices:
arr[i] = 1
yield arr

def test_encode(n, m):
for arr in gen_code(n, m):
dst = enc(arr)
print(dst)

def test_decode(n, m):
for arr in gen_code(n, m):
dst = enc(arr)
src = dec(dst, n, n-m)
if src != arr:
print(src)
break

#test_decode(10,3)
#test_decode(10,3)
#enc([0, 1, 0, 1, 1, 0, 0, 0, 0])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
// 文章3运用到m-of-n code的C++实现
// 这里用到一个bitstream类,其实现略去,可通过调用方法名来猜测其含义 =w=
#include "bitstream.h"
#include <math.h>
#include <stdint.h>
#include <iostream>

#define FIRST_QTR 0x40000000
#define THIRD_QTR 0xc0000000
#define HALF 0x80000000


// Return an approx upper bound of log2(binom(n,k))
inline double log_binom(int n, int k)
{
int j = n - k;
int min = j < k ? j : k;
double sum = 0;
if (min <= 5) // Small case, where Stirling's approximation works not well
for (int i = 0; i < min; i++)
sum += log2((double)(n - i) / (min - i));
else // Use Stirling's approximation
sum = (1./k + 1./j - 1./(n+1)) / 12. + (n+0.5) * log2(n) - (k+0.5) * log2(k) - (j+0.5) * log2(j) - 1.325747;
// Here 1.325747 = log2(sqrt(2pi))
return sum;
}

// Since the length of encoded string is not known by the decoder,
// we calculate an upper bound of it, and fill it by 0 until the length of it is exact the upper bound.
inline int code_bits_bound(int n, int n0)
{
int n1 = n - n0;
if (n0 == 0 || n1 == 0)
return 0;
/*
Overhead bits because of integer implementation.
Details: "A Coding Scheme for rn-out-of-n Codes" - https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=58748
Here we assume 8 < n <= 2^26, so
double T = n0 * log(n1) + n1 * log(n0) + (n0 / n1 + n1 / n0) / 2 + 1.57722 * n;
< n*(log(n) + 2.1) < 2n*log(n)
double epsilon = log2(1+log(1/(1-T/2^32))) = log2(1+T/2^32)=T/2^32/log(2) < n*log(n)/log(2)/2^31=n*log2(n)/HALF
*/
double epsilon = n * log2(n) / HALF;

double total_length = log_binom(n, n0) + epsilon;
return (int)ceil(total_length);
}

void arithmetic_encode(bitstream *src, bitstream *dst, int n, int n0)
{
if (n > (1 << 26))
{
std::cout << "Arithmetic Encoding Error: source string is too long!" << std::endl;
return;
}
uint32_t low = 0, high = 0xffffffff;
int bits_to_follow = 0;
int counter = 0;
int code_bits = code_bits_bound(n, n0);
while (n0 < n && n0 > 0)
{
uint32_t sep = (uint32_t)((((uint64_t)high - low + 1) * n0 * 2 + n) / (n * 2)); // = round(step * n0 / n)
n--;
if (src->read())
low += sep;
else
{
n0--;
high = low + sep - 1;
}
while (1)
{
if (high < HALF || low >= HALF) // E1/E2 scaling
{
counter += bits_to_follow + 1;
dst->write(low >= HALF);
dst->write_repeatedly(high < HALF, bits_to_follow);
bits_to_follow = 0;
}
else if ((FIRST_QTR <= low) && (high < THIRD_QTR)) // E3 scaling
{
bits_to_follow++;
low ^= THIRD_QTR;
high ^= THIRD_QTR;
}
else
break;
low <<= 1; // Leftmost bit 1 should be discarded automatically
high = (high << 1) | 1;
}
}
counter++;
dst->write(1);
if (counter > code_bits)
{
std::cout << "Arithmetic Encoding Error: Invalid code_bits found!" << std::endl;
return;
}
dst->write_repeatedly(0, code_bits - counter); // so that the length of encoded string is exact code_bits_bound(n, n0).
}

inline unsigned char get_bit(bitstream *dst, int &counter, int code_bits)
{
if (counter >= code_bits)
return 0;
counter++;
return dst->read();
}

void arithmetic_decode(bitstream *src, bitstream *dst, int n, int n0)
{
uint32_t low = 0, high = 0xffffffff;
int counter = 0;
int code_bits = code_bits_bound(n, n0);
uint32_t value = 0;
for (int i = 0; i < 32; i++)
value = (value << 1) | get_bit(dst, counter, code_bits);
while (n0 < n && n0 > 0)
{
uint32_t sep = (uint32_t)((((uint64_t)high - low + 1) * n0 * 2 + n) / (n * 2)); // = round(delta * n0 / n)
n--;
if (low + sep <= value)
{
src->write(1);
low += sep;
}
else
{
src->write(0);
n0--;
high = low + sep - 1;
}
while (1)
{
if (high >= HALF && low < HALF)
{
if ((FIRST_QTR <= low) && (high < THIRD_QTR))
{
low ^= THIRD_QTR;
high ^= THIRD_QTR;
value ^= THIRD_QTR;
}
else
break;
}
low <<= 1; // Left bit 1 should be discarded automatically
high = (high << 1) | 1;
value = (value << 1) | get_bit(dst, counter, code_bits);
}
}
if (counter < code_bits)
{
std::cout << "Arithmetic Decoding Warning: Unusual case found!" << std::endl;
for (int i = 0; i < code_bits - counter; i++)
dst->read();
}
src->write_repeatedly(0, n0);
n -= n0;
src->write_repeatedly(1, n);
}

  1. https://cstheory.stackexchange.com/questions/19313/optimal-encoding-of-k-subsets-of-n ↩︎ ↩︎

  2. Bodden E, Clasen M, Kneis J. Arithmetic coding revealed[J]. McGill University. Sable Research Group www. sable. mcgill. ca, 2007. ↩︎ ↩︎

  3. Ramabadran T V. A coding scheme for m-out-of-n codes[J]. IEEE Transactions on Communications, 1990, 38(8): 1156-1163. ↩︎

  4. Jones C. An efficient coding system for long source sequences[J]. IEEE Transactions on Information Theory, 1981, 27(3): 280-291. ↩︎

  5. Witten I H, Neal R M, Cleary J G. Arithmetic coding for data compression[J]. Communications of the ACM, 1987, 30(6): 520-540. ↩︎ ↩︎

  6. Tian C, Vaishampayan V A, Sloane N J A. A coding algorithm for constant weight vectors: A geometric approach based on dissections[J]. IEEE Transactions on Information Theory, 2009, 55(3): 1051-1060. ↩︎ ↩︎

  7. Skachek V, Immink K A S. Constant weight codes: an approach based on Knuth’s balancing method[J]. IEEE Journal on Selected Areas in Communications, 2014, 32(5): 909-918. ↩︎ ↩︎ ↩︎