1 实验环境

语言：C++

2 实验内容

2.1 编程实现 $(m,n)$ 群码

2.1.1 输入输出

输入：$m,n$ 和一致性校验矩阵 $H$ 中的 $H_{m \times r}$ 部分，其中 $r=n-m$.

输出：所有的码字，按原码对应的十进制整数升序排列.

2.1.2 设计思路及核心算法

  1. 首先输入 $m,n$：

    std::cin >> m >> n;

  2. 调用 inputH(m,n) 函数输入一致性校验矩阵 $H$：

    BitMatrix H = inputH(m, n);

    inputH(m,n) 函数代码如下：

    uint64_t inputBit() {
    std::string line;
    std::cin >> line;
    uint64_t x = 0;
    int k = 0;
    for (size_t i = 0; i < line.size(); i++)
        if (line[i] == '1' || line[i] == '0')
            x |= (line[i] - '0') << k++;
    return x;
}
BitMatrix inputH(int m, int n) {
    BitMatrix H(n, n - m);
    for (int i = 0; i < m; i++) H.d[i] = inputBit();
    for (int i = m; i < n; i++) H.d[i] = 1 << (i - m);
    return H;
}

    inputH 函数传入行列参数m,n，声明一个 BitMatrix 类，用来存放一致性校验矩阵的 $m \times r$部分.

  3. 调用 HtoH2(H) 函数生成 $m \times n$ 的编码矩阵： $$ \left[ \begin{matrix} I_m & H_{m\times r} \ \end{matrix} \right] $$

    BitMatrix HtoH2(BitMatrix &H) {
    BitMatrix H2(H.n - H.m, H.n);
    for (int i = 0; i < H.n - H.m; i++) H2.d[i] = 1 << i;
    for (int i = 0; i < H.n - H.m; i++)
        for (int j = 0; j < H.m; j++)
            H2.d[i] |= (H.d[i] >> j & 1) << (j + H.n - H.m);
    return H2;
}

  4. 定义一个 Bm 矩阵用于存放原码，矩阵每行从上到下按原码对应的十进制递增顺序排列，并进行运算：

    $$ e_H(B^m)=B^m\times
    \left[ \begin{matrix} I_m & H_{m\times r} \ \end{matrix} \right] $$

    求出群码对应的矩阵并输出.

    BitMatrix H2 = HtoH2(H), Bm(1 << m, m);
for (int i = 0; i < (1 << m); i++)
    for (int j = 0; j < m; j++)
        Bm.d[i] |= (i >> j & 1) << (m - j - 1);
BitMatrix N = Bm * H2;
std::cout << "1. e(B^m):" << std::endl;
for (int i = 0; i < (1 << m); i++) {
    std::cout << "  e(";
    printBit(Bm.d[i], m);
    std::cout << ") = ";
    printBit(N.d[i], n);
    std::cout << std::endl;
}

2.2 构造陪集表

2.2.1 设计思路及核心算法

生成 coset table 的算法原理如下：

  1. 首先，定义了一个空的哈希表 cosets 用于存储 coset 的信息，以及一个空的哈希表 visited 用于标记已经访问过的向量；

  2. 创建一个向量 numbers，其中存储了所有可能的向量（共 2^n 个向量），并按照向量中 1 的个数进行排序，从少到多；

  3. 遍历 numbers 中的每个向量 i，对于每个向量，执行以下步骤：

    a. 创建一个空的 coset 对象 coset，用于存储属于同一个 coset 的向量；

    b. 对于矩阵 N 的每一行（表示每个 coset），将当前向量 i 与该行的向量 x 进行异或运算，得到新的向量 x' = N.d[j] ^ uint64_t(i)；

    c. 检查 visited 哈希表，如果在 visited 中找到了向量 x'，说明该向量已经属于某个 coset，跳过当前循环；

    d. 如果在 visited 中没有找到向量 x'，则将 x' 添加到当前 coset 中，并将 x' 添加到 visited 哈希表中，值为当前 coset 的索引；

    e. 将当前 coset 的陪集头 e 更新为 coset 中具有最少 1 的个数的向量；

    f. 如果当前 coset 的陪集头 e 在 cosets 哈希表中不存在，则将当前 coset 添加到 cosets 中，并将 e 对应的特征值（通过一致性校验矩阵 H 计算）添加到 syncdrome 哈希表中，值为当前 coset 的索引；

    g. 检查 cosetLeaders 的大小，如果已经达到了 $2^{n-m}$（即余类的数量），则停止循环.

  4. 最后，得到了所有的 coset 信息，并按照 coset 陪集头对应的十进制递增的顺序输出.

通过上述算法，可以生成所有的 coset，并找到每个 coset 的陪集头和特征值.

2.2.2 BitMatrix 类

该程序将译码时所需的比特矩阵进行了封装，能够实现乘法 (元素对应异或) 操作.

class BitMatrix {
public:
    int n, m;
    uint64_t d[100];
    BitMatrix(int n, int m): n(n), m(m) {
        memset(d, 0, sizeof(uint64_t) * n);
    }
    BitMatrix operator*(const BitMatrix& rhs) const {
        assert (m == rhs.n);
        BitMatrix res(n, rhs.m);
        for (int i = 0; i < n; i++)
            for (int k = 0; k < m; k++)
                if (d[i] >> k & 1)
                    res.d[i] ^= rhs.d[k];
        return res;
    }
    void print() const {
        for (int i = 0; i < n; i++) {
            printBit(d[i], m);
            std::cout << std::endl;
        }
    }
};

2.2.3 Coset 类

Coset 类由陪集、陪集头及添加操作组成，陪集中的元素使用一维向量进行存储，陪集头用整型变量单独存储，方便调用. 添加操作在每一次执行时更新一次陪集头.

class Coset {
public:
    std::vector<uint64_t> d; // coset
    uint64_t e; // coset leader
    Coset() : e(0) {}
    void add(uint64_t x) {
        if (d.empty())
            d.push_back(e = x);
        else {
            d.push_back(x);
            if (bitcount(x) < bitcount(e)) e = x;
        }
    }
};

2.3 极大似然译码

基于 2.2 构造陪集表 求出的左陪集表，用和函数 $e_H$ 相关联的极大似然译码函数 $d$ ，对 $x_t\in B^n$ 进行译码.

2.3.1 输入输出格式

  • 输入：首行是译码个数 $N$ ，整数格式；
  • 然后是 $N$ 行比特串 $x_t \in B^n$ ，字符串格式；
  • 输出：$N$行译码结果，格式为 $d(x_t) = b_1 b_2 … b_m$ . 如：d(001011)=001.

2.3.2 设计思路及核心算法

  1. 输入 $N$ ，然后输入待解码的 $N$ 个比特串. 输入后需要进行字符串转整数操作；

    std::cin >> k;
for (int i = 0; i < k; i++) {
    std::cerr << "Enter the word #" << i + 1 << ": ";
    uint64_t x = inputBit();
/*
    ......................
*/

    inputBit() 函数如下：

    uint64_t inputBit() {
    std::string line;
    std::cin >> line;
    uint64_t x = 0;
    int k = 0;
    for (size_t i = 0; i < line.size(); i++)
        if (line[i] == '1' || line[i] == '0')
                x |= (line[i] - '0') << k++;
    return x;
}

  2. 确定输入的 $x_t$ 位于哪个 coset 中，可通过 visited[x] 得到所属陪集编号；

  3. 取出所属陪集的陪集头 $\epsilon$ ，与 $x_t$ 进行异或运算，得到的结果必然在群码对应的矩阵 $N$ 中找到，由此得到解码后结果.

    核心代码如下：

    std::vector<uint64_t > cosetLeaders;
    for (auto i : numbers) {
        Coset coset;
        bool notFound = true;
        for (int j = 0; j < N.n; j++) {
            uint64_t x = N.d[j] ^ uint64_t(i);
            if (visited.find(x) != visited.end()) {
                notFound = false;
                break;
            }
            visited[x] = cosetLeaders.size();
            coset.add(x);
        }
        if (notFound) {
            cosets[coset.e] = coset;
            syndrome[getSyndrome(coset.e, H)] = cosetLeaders.size();
            cosetLeaders.push_back(coset.e);
        }
        // if cosetLeaders.size() = 2^(n-m), then stop. 
        // Because there are only 2^(n-m) cosets.
        if (cosetLeaders.size() == (1 << (n - m))) break;
    }
/*
............
*/
printBit(Bm.d[NtoBm[cosetLeaders[visited[x]] ^ x]], m);

2.4 基于特征值的译码

基于 2.2 构造陪集表 求出的左陪集表，利用函数 $f_H:B^n \rightarrow B^r，f_H(x) = x\cdot H$ 求出每个陪集头对应的特征值 (Syndrome) ，并利用特征值，对$x_t \in B^n$ 进行译码.

2.4.1 输入输出格式

  • 输入：首行是译码个数N，整数格式；
  • 然后是N行比特串 $x_t \in B^n$ ，字符串格式；
  • 输出：N行译码所用的特征值、陪集头和译码结果，格式为 $f(x_t)$ = *** ，$\epsilon$ = ***，$d(xt) = b_1b_2…b_m$.

注：本程序采用的输出方式，对于每一个待解码二进制串输入，输出结果分为两行，第一行为采用极大似然译码计算的结果，结果前用序号 (a) 进行区分，第二行为基于特征值译码的结果，结果前用序号 (b) 区分，例如：

(a) d(10111) = 10
(b) f(10111) = 100, ε = 00100, d(10111)

2.4.2 设计思路及核心算法

  1. 输入 N ，然后输入待解码的 N 个比特串. 输入后需要进行字符串转整数操作；

  2. 通过函数 $f_H(x_t) = x_t \cdot H$ 计算 $x_t$ 特征值；

    std::cout << " (b) f("; //特征值
printBit(x, n);
uint64_t s = getSyndrome(x, H);
std::cout << ") = ";
printBit(s, n - m);

    getSyndrome() 函数细节：

    uint64_t getSyndrome(uint64_t x, BitMatrix &H) {
    BitMatrix X(1, H.n);
    X.d[0] = x;
    return (X * H).d[0];
}

  3. 根据所得特征值在哈希表中对应到陪集头；

    std::cout << ", ε = "; //陪集头
printBit(cosetLeaders[syndrome[s]], n);

  4. 陪集头与待译码比特串进行异或操作，得到的结果必然在群码对应的矩阵 $N$ 中找到，由此得到解码后结果.

    std::cout << ", d("; //译码结果
printBit(x, n);
std::cout << ") = ";
printBit(Bm.d[NtoBm[cosetLeaders[syndrome[s]] ^ x]], m);
std::cout << std::endl;

3 代码优化

我们在四个地方优化了原本的代码，使得程序更加高效.

  1. 由于实验中满足 $1<m<n<7$，所以使用 uint64_t 替代数组，并使用位运算进行相应计算，同时降低了空间复杂度和时间复杂度；

    其中，矩阵乘法的核心代码如下：

    // main.cpp, Line 30-33
for (int i = 0; i < n; i++)
    for (int k = 0; k < m; k++)
        if (d[i] >> k & 1)
            res.d[i] ^= rhs.d[k];

    普通的矩阵乘法的时间复杂度是 $O(n^3)$，而使用位运算后，时间复杂度降低到了 $O(n^2)$.

  2. 当已经找到了所有的陪集时，就可以停止循环，不需要再继续遍历；

    // main.cpp, Line 132
if (cosetLeaders.size() == (1 << (n - m))) break;

  3. 使用哈希表记录每一个码字所属的陪集，在解码的时候可以直接查找，而不需要再次遍历；

    // main.cpp, Line 123, 记录每一个码字所属的陪集
visited[x] = cosetLeaders.size();

    // main.cpp, Line 155, 直接解码
printBit(Bm.d[NtoBm[cosetLeaders[visited[x]] ^ x]], m);

    该哈希表还可以快速地判断当前陪集是否重复，减少了计算量；

    // main.cpp, Line 119-122, 判断当前陪集是否重复
if (visited.find(x) != visited.end()) {
    notFound = false;
    break;
}

    考虑到 $1<m<n<7$，所以码字的范围最大不超过 $2^8 = 256$，所以实际上可以使用数组代替哈希表，降低部分常数.

  4. 特征值解码方式也使用了类似的方式；

    // main.cpp, Line 128, 记录每一个特征值对应的陪集(头)
syndrome[getSyndrome(coset.e, H)] = cosetLeaders.size();

    // main.cpp, Line 155, 直接获取特征值对应的陪集头并进行异或运算
printBit(Bm.d[NtoBm[cosetLeaders[syndrome[s]] ^ x]], m);

  5. 使用了 lowbit 快速计算最低位的 1 的位置，加快了 bitcount 的速度.

    // main.cpp, Line 17
while (n) n -= n & -n, s++;

4 实验成果示例

$ ./main # 运行
2 5 # 输入m n
011
101 #这两行是一致性校验矩阵中的 H_{m*r} 部分
1. e(B^m):
  e(00) = 00000
  e(01) = 01101
  e(10) = 10011
  e(11) = 11110
2. Cosets:
  00000 | 00000 01101 10011 11110
  00001 | 00001 01100 10010 11111
  00010 | 00010 01111 10001 11100
  00100 | 00100 01001 10111 11010
  01000 | 01000 00101 11011 10110
  10000 | 10000 11101 00011 01110
  00110 | 00110 01011 10101 11000
  01010 | 01010 00111 11001 10100
Enter the number of words to decode: 3 # 在输出码字和陪集之后，在这里输入你要测试的次数
Enter the word #1: 10111
 (a) d(10111) = 10
 (b) f(10111) = 100, ε = 00100, d(10111) = 10 #输出译码的结果
······

附录：源代码

#include <iostream>
#include <string>
#include <vector>
#include <unordered_map>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <cstdint>
#include <cassert>
void printBit(uint64_t x, int n) {
    for (int i = 0; i < n; i++)
        std::cout << (x >> i & 1);
}
template <typename T>
int bitcount(T n) {
    int s = 0;
    while (n) n -= n & -n, s++;
    return s;
}
class BitMatrix {
public:
    int n, m;
    uint64_t d[100];
    BitMatrix(int n, int m): n(n), m(m) {
        memset(d, 0, sizeof(uint64_t) * n);
    }
    BitMatrix operator*(const BitMatrix& rhs) const {
        assert (m == rhs.n);
        BitMatrix res(n, rhs.m);
        for (int i = 0; i < n; i++)
            for (int k = 0; k < m; k++)
                if (d[i] >> k & 1)
                    res.d[i] ^= rhs.d[k];
        return res;
    }
    void print() const {
        for (int i = 0; i < n; i++) {
            printBit(d[i], m);
            std::cout << std::endl;
        }
    }
};
class Coset {
public:
    std::vector<uint64_t> d; // coset
    uint64_t e; // coset leader
    Coset() : e(0) {}
    void add(uint64_t x) {
        if (d.empty())
            d.push_back(e = x);
        else {
            d.push_back(x);
            if (bitcount(x) < bitcount(e)) e = x;
        }
    }
};
uint64_t inputBit() {
    std::string line;
    std::cin >> line;
    uint64_t x = 0;
    int k = 0;
    for (size_t i = 0; i < line.size(); i++)
        if (line[i] == '1' || line[i] == '0')
            x |= (line[i] - '0') << k++;
    return x;
}
BitMatrix inputH(int m, int n) {
    BitMatrix H(n, n - m);
    for (int i = 0; i < m; i++) H.d[i] = inputBit();
    for (int i = m; i < n; i++) H.d[i] = 1 << (i - m);
    return H;
}
BitMatrix HtoH2(BitMatrix &H) {
    BitMatrix H2(H.n - H.m, H.n);
    for (int i = 0; i < H.n - H.m; i++) H2.d[i] = 1 << i;
    for (int i = 0; i < H.n - H.m; i++)
        for (int j = 0; j < H.m; j++)
            H2.d[i] |= (H.d[i] >> j & 1) << (j + H.n - H.m);
    return H2;
}
uint64_t getSyndrome(uint64_t x, BitMatrix &H) {
    BitMatrix X(1, H.n);
    X.d[0] = x;
    return (X * H).d[0];
}
int main() {
    int m, n;
    std::cin >> m >> n;
    assert (1 < m && m < n && n < 64);
    BitMatrix H = inputH(m, n);
    BitMatrix H2 = HtoH2(H), Bm(1 << m, m);
    for (int i = 0; i < (1 << m); i++)
        for (int j = 0; j < m; j++)
            Bm.d[i] |= (i >> j & 1) << (m - j - 1);
    BitMatrix N = Bm * H2;
    std::cout << "1. e(B^m):" << std::endl;
    for (int i = 0; i < (1 << m); i++) {
        std::cout << "  e(";
        printBit(Bm.d[i], m);
        std::cout << ") = ";
        printBit(N.d[i], n);
        std::cout << std::endl;
    }
    std::unordered_map<uint64_t, Coset> cosets;
    std::unordered_map<uint64_t, int> visited, NtoBm, syndrome;
    for (int i = 0; i < N.n; i++) NtoBm[N.d[i]] = i;
    std::vector<int> numbers;
    for (int i = 0; i < (1 << n); i++) numbers.push_back(i);
    std::sort(numbers.begin(), numbers.end(), [&](int a, int b) {
        int ca = bitcount(a), cb = bitcount(b);
        return ca == cb ? a > b : ca < cb;
    });
    std::vector<uint64_t > cosetLeaders;
    for (auto i : numbers) {
        Coset coset;
        bool notFound = true;
        for (int j = 0; j < N.n; j++) {
            uint64_t x = N.d[j] ^ uint64_t(i);
            if (visited.find(x) != visited.end()) {
                notFound = false;
                break;
            }
            visited[x] = cosetLeaders.size();
            coset.add(x);
        }
        if (notFound) {
            cosets[coset.e] = coset;
            syndrome[getSyndrome(coset.e, H)] = cosetLeaders.size();
            cosetLeaders.push_back(coset.e);
        }
        // if cosetLeaders.size() = 2^(n-m), then stop.
        // Because there are only 2^(n-m) cosets.
        if (cosetLeaders.size() == (1 << (n - m))) break;
    }
    std::cout << "2. Cosets:" << std::endl;
    for (auto i: cosetLeaders) {
        std::cout << "  ";
        printBit(i, n);
        Coset coset = cosets[i];
        std::cout << " |";
        for (size_t i = 0; i < coset.d.size(); i++) {
            std::cout << " ";
            printBit(coset.d[i], n);
        }
        std::cout << std::endl;
    }
    std::cerr << "Enter the number of words to decode: ";
    int k;
    std::cin >> k;
    for (int i = 0; i < k; i++) {
        std::cerr << "Enter the word #" << i + 1 << ": ";
        uint64_t x = inputBit();
        std::cout << " (a) d(";
        printBit(x, n);
        std::cout << ") = ";
        printBit(Bm.d[NtoBm[cosetLeaders[visited[x]] ^ x]], m);
        std::cout << std::endl;
        std::cout << " (b) f(";
        printBit(x, n);
        uint64_t s = getSyndrome(x, H);
        std::cout << ") = ";
        printBit(s, n - m);
        std::cout << ", ε = ";
        printBit(cosetLeaders[syndrome[s]], n);
        std::cout << ", d(";
        printBit(x, n);
        std::cout << ") = ";
        printBit(Bm.d[NtoBm[cosetLeaders[syndrome[s]] ^ x]], m);
        std::cout << std::endl;
    }
    return 0;
}