目录

有限域可能是RS-FEC中最难懂的部分了，其实就数学要求而言，高中知识也许就足矣，但是现有的很多介绍RS-FEC的文章对这块知识的介绍比较简略，很难有比较系统的理解。为了学习有限域，是有必要找本数学书来看看，比如冯克勤的<<走向数学从书17 - 有限域>>，看完第一章有限域就可以了，45页纸，虽然不比看小说那么轻松，但应该不存在哪里有什么看不懂的地方，因为确实只需要高中的数学知识。这篇文章的目的不是去重复书本的内容，而是综述性的介绍一些来龙去脉和知识点的关联，类似一个读书笔记，另外还会结合”A Tutorial on Reed-Solomon Coding for Fault-Tolerance in RAID-like Systems”这篇论文来介绍一下具体的实现。

有限域从名称来来看，一个是有限的元素，另外一个是域，域可以理解为在元素构成的集合上可以做加减乘除的运算，并且满足结合律，交换律和分配律，还有满足0元素和1元素的相关特性（具体定义见书14页）。FEC如果是在无限域上运算，运算可能会溢出，精度也会有损失，有限域上因为元素有限，数值的长度可以固定，不存在溢出的情况，而且不存在会有浮点运算之类的精度问题，另外还有个好处在后面具体的知识介绍和剧透实现可以看到，有限域上的运算性能相对无限域在理论上（乘法除法可以转换为加减法，本质上类似在实数域我们也有指数和对数运算，但这个表相比元素无限个数的实数域会小很多）和实现上都可以有很大的优化（因为个数有限，可以预先计算好，运算直接查表来的到结果）。

P元有限域

对于每个素数p，可以构作出一个p元有限域Fp

p元有限域的构建

为什么要是素数? 证明之前先举个直观的例子。取模操作大家都比较熟悉了，对于m为>=2的正整数，a是整数，定义如下的一个集合，集合里面的元素[a]表示模m与a同余的所有整数构成的集合(集合里面的元素也是集合), [a]叫做模m的一个同余类。那这样就一共存在m个同余类: [0], [1], … [m-1]。

然后定义如下操作（对，是定义，不要去纠结操作的具体的现实意义，很多理论是到了几百年后才有应用场景的，有限域理论在三百多年前就出现了，到了信息时代才得到应用。不是说没有，而是说重点不在那。其实这些定义是可以根据上面的同余类的定义推导出来的）:

1
2
3
[a] + [b] = [a+b]
[a] - [b] = [a-b]
[a] . [b] = [a.b]

可以证明基于上面的定义，加法和乘法满足交换律，结合律，分配律，而且[0]和[1]的作用和我们熟悉的有限域下一样。

其次对于加法和减法，两者互为逆运算

1
2
3
4
[a] + [b] = [c] <=> [c] - [a]
                  = [[a] + [b]] - [a]
                  = [[a]] + [[b]] - [a]
                  = [a] + [b] - [a] = [b]

剩下最后一点有待确认: 是否存在乘法的逆运算，即除法运算。上面没有类似加减乘法来定义[a]/[b] = [a/b]，因为首先a/b这个运算对于整数而言就没有意义，因为结果不一定是整数。根据同余类概念的定义[x]=[a]/[b] ([b]!=[0])相当于求解同余方程b.x=a(mod)m。对于乘法的逆运算，应该是对每个同余类[b]，存在同余类[c]，是的[b].[c] = [1]

1
[a]/[b] = [a].[1]/[b] = [a].[b].[c]/[b] = [a].[c]

例如对于m=4，对于[2]，就找不到对应的[a]使得[a].[2]=[1]。

1
2
3
4
[0].[2] = [0]
[1].[2] = [2]
[2].[2] = [0]
[3].[2] = [2]

具体的证明见书中第7页引理1.1.6：设m>=2, [a]为模m的同余类，则存在同余类[b]使得[a].[b] = [1]的充要条件是a,m互素。并且当a,m互素时，只有一个同余类[b]使得[a].[b]=[1]。

由引理知，对于非素数的m，总可以找到一个正因子d>1(d和m不是互素关系)，也就是说[d]的逆元素不存在。

所以上面也就证明了最开始的结论: 对于每个素数p，可以构作出一个p元有限域Fp

对于p=5，可以有如下的F5的乘法表:

1
2
3
4
5
6
7
8
9
10
11
      [0]  [1]  [2]  [3]  [4]
     -------------------------
[0] | [0]  [0]  [0]  [0]  [0]
    |
[1] | [0]  [1]  [2]  [3]  [4]
    |
[2] | [0]  [2]  [4]  [1]  [4]
    |
[3] | [0]  [3]  [1]  [4]  [2]
    |
[4] | [0]  [4]  [4]  [2]  [1]

p元有限域的特性

后面为了书写方便，有限域内的数会省略掉”[]”。下面是一些p元有限域的一些重要的特性:

Fp中的每个元素a，均有a^^p = a，如果a!=0，则a^^(p-1) = 1

例如: 3^^4 = 3.3.3.3 = 4.3.3 = 2.3 = 1

证明: 因为共有p个不同的元素，所以当a!=0时，a, 2a, 3a, (p-1)a也分别是p个不同的元素, 所以1.2 … p-1 = a. 2a. 3a. (p-1)a，消去1.2 … p-1，得到a^^(p-1) = 1。

Fp中非0元素a的阶定义为最小的正整数n，有a^^n = 1。因为上面已经证明了a^^(p-1)=1，所以n一定是存在的。a的阶n一定是p-1的因子。

另外，如果Fp中存在元素g，且g的阶是p-1，那么1=g^^0, g^^1, g^^(p-2)分别对于p-1个不同的元素。（见11页“注纪”里的证明），也就是说非零元素都可以由g来生成，这个就称为本原元素，有了这个，Fp的数学结构就简单了。

对于上面p=5的例子, 3就是一个本原元素:

1
2
3
4
5
3^^0 = 1
3^^1 = 3
3^^2 = 3.3 = 4
3^^3 = 3.3.3 = 4.3 = 2
3^^4 = 3.3.3.3 = 4.3.3 = 2.3 = 1

有了本原元素指数方式来表示每一个元素，乘法运算和除法运算就可以转换为加法和减法运算了。

问题到这里并没有结束，虽然我们知道对于素数p存在p元有限域，但给定任意一个素数，除了穷举法，我们没办法快速找出其本原元素。如果大家看过RS-FEC相关的文献，就会发现里面有用到“生成多项式”，后面将介绍怎么会和多项式有关系的。

任意有限域

元素个数

有限域F的元素个数必为p^^n，其中p为素数, n>=1。由前面可知，Fp是F的一个子域，利用最小生成集，可以将F中的元素表示为x = a1.x1 + ... + an.xn，可以证明对于任意的x，ai系数是唯一的。ai的组合共有p^^n个，所以共有p^^n个元素。

存在性

对于每个素数幂p^^n，均存在p^^n元有限域(多项式x^^(p^^n)-x的在Fp的代数闭包中的所有根形成p^^n元有限域)。这里的证明用到了重根判别法，然后证明根构成一个域。

唯一性

对于Fp的代数闭包F，其中恰好由一个p^^n元有限域。

令q=p^^n，F中非零元素为x[1], …, x[q-1]（这里借用数组的[]来表示下标），每个元素乘以非零元素a，由x[1]. … . x[q-1] = a^^(q-1).x[1]. … .x[q-1]。则由a^^(q-1)=1，即a^^q=a。从而每个元素都是f(x) = x^^q-x = x^^(p^^n)-x的根，且只有p^^n个根，从而证明了唯一性。

本原元素的存在性

q=p^^n，q-1阶元素叫做q元域Fq的本原元素。有了本原元素后，做乘法操作就很方便了。但加法应该怎么做呢？

构建有限域

Fp中的n次不可约首一多项式f(x)，也是RS-FEC文献中的生成多项式。设a为f(x)的一个根，那么集合

1
2
S = {c[0].a^^(n-1) + c[1].a^^(n-2) + ... + c[n-1]
其中c[0], ..., c[n-1] 都是Fp中的元素

是一个q元域，q=p^^n。
这里的证明思路是先证明S中有p^^n个不同的元素，其次再根据域的定义来证明S中的元素上的操作满足定义。

有了这个结论后，加法操作就是多项式的加法操作（系数属于Fp）。

例子

F3[x]中的不可约多项式f(x)=x^^2 + x + 2（如果可约，那一定会有一次多项式因子，使得其在F3中有根，但f(0)=2, f(1)=1, f(2)=2，都不为0，所以f(x)不可约）。根据上面的定理，设a为f(x)在扩域的一个根，则有a^^2 + a + 2 = 0 <=> a^^2 = 2a + 1，那么对于n=2, p=3, q=p^^n=3^^2=9元有限域F9有如下9个元素:

1
S = {c[0].a + c[1]} = {0, 1, 2, a, a+1, a+2, 2a, 2a+1, 2a+2}

此外可以看到a是8阶元素，即是F9的本原元素:

1
2
3
4
5
6
7
a^^2 = 2a + 1
a^^3 = (2a + 1).a = 2a^^2 + a = 2(2a + 1) + a = 2a + 2
a^^4 = (2a + 2).a = 2a^^2 + 2a = 2(2a + 1) + 2a = 2
a^^5 = 2a
a^^6 = 2a^^2 = 2(2a + 1) = a + 2
a^^7 = (a + 2).a = a^^2 + 2a = 2a + 1 + 2a = a + 1
a^^8 = (a + 1).a = a^^2 + a = 2a + 1 + a = 1

所以F9中的元素有如下两种表示方式，纪c[0].a + c[1]为(c[0], c[1]):

1
2
3
   0 = (0, 0)       1 = (0, 1)       a = (1, 0)
a^^2 = (2, 1)    a^^3 = (2, 2)    a^^4 = (0, 2)
a^^5 = (2, 0)    a^^6 = (1, 2)    a^^7 = (1, 1)

在计算乘法的时候可以用左边的表达方式操作:

1
(2a + 1).(a + 2) = (2, 1).(1, 2) = a^^2 . a^^6 = a^^8 = 1 = (0, 1)

加法可以用右边的表达方式来操作:

1
a^^3 + a^^5 = (2, 2) + (2, 0) = (1, 2) = a^^6

代码实现

实现中都会取p=2, 对应的n次不可约多项式分别为:

1
2
3
4
5
6
n = 2:    x^^2 + x + 1
n = 4:    x^^4 + x + 1
n = 8:    x^^8 + x^^4 + x^^3 + x^^2 + 1
n = 16:   x^^16 + x^^12 + x^^3 + x + 1
n = 32:   x^^32 + x^^22 + x^^2 + x + 1
n = 64:   x^^64 + x^^4 + x^^3 + x + 1

下面以n=4为例，GF(p^^n) = GF(2^^4) = 16，即16元有限域，设a是x^^4 + x + 1的一个解，即有a^^4 = a + 1，下面可用看的a也为本原元素（a^^15＝1，只有本原元素才能生成所有的非零元素，也才能用于指数方式的表达），16元有限域的元素就仅仅限于前16项，a^^15和a^^0是同一个元素。

指数方式 多项式方式 二进制
0 0 0000
a^^0 1 0001
a^^1 a 0010
a^^2 a^^2 0100
a^^3 a^^3 1000
a^^4 a + 1 0011
a^^5 a^^2 + a 0110
a^^6 a^^3 + a^^2 1100
a^^7 a^^3 + a + 1 1011
a^^8 a^^2 + 1 0101
a^^9 a^^3 + a 1010
a^^10 a^^2 + a +1 0111
a^^11 a^^3 + a^^2 + a 1110
a^^12 a^^3 + a^^2 + a + 1 1111
a^^13 a^^3 + a^^2 + 1 1101
a^^14 a^^3 + 1 1001
a^^15 1 0001

上面可以很清楚的看到二进制方式只是多项式方式略去a后的一种缩写的表示方式而已。而指数方式则是根据a^^4 = a + 1带入运算的到的:

1
2
3
4
a^^4 = a + 1
a^^5 = a^^2 + a
a^^6 = a^^3 + a^^2
a^^7 = a^^4 + a^^3 = a^^3 + a + 1

如果要讲上面的计算过程转换为代码:

a^^7实际上可以a^^6 . a，也就是a^^6的二进制方式左移一位，从1100(a^^3 + a^^2)变成11000(a^^4 + a^^3)，因为a^^4 + a + 1 = 0，任何元素加上0不变，所以a^^4 + a^^3 + a^^4 + a + 1，转换为XOR操作(11000 ^ 10011 = 01011)，即a^^3 + a + 1。

可以看到加法操作其实也就是系数为F2的多项式加法，也就是系数的XOR运算，乘法可以先转换为指数方式，然后做（模2^^n-1=15）加法，然后转换为多项式方式。

上面的过程转换为c语言的实现，首先是要生成两个表:

  • gflog用来将“多项式方式 => 指数方式”（根据二进制方式，范围从1到2^^n-1，查找对应的指数），
  • gfilog用来将“指数方式 => 多项式方式”（根据指数方式，范围从0到2^^n-2，查找对应的二进制方式）。
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
unsigned int prim_poly_4 = 023; //x^^4 + x + 1 => 0b10011 = 023
unsigned int prim_poly_8 = 0435;
unsigned int prim_poly_16 = 0210013;
unsigned short * gflog, * gfilog;

int setup_tables(int w)
{            
    unsigned int b, log, x_to_w, prim_poly;
         
    switch (w) {
        case 4:
            prim_poly = prim_poly_4;
            break;
         
        case 8:
            prim_poly = prim_poly_8;
            break;
            
        case 16:
            prim_poly = prim_poly_16;
            break;
            
        default:
            return -1;
    }
     
    x_to_w = 1 << w;
    gflog = (unsigned short *) malloc (sizeof(unsigned short) * x_to_w);
    gfilog = (unsigned short *) malloc (sizeof(unsigned short) * x_to_w);
         
    b = 1; 
    // b - poly form, log - exp form
    for (log = 0; log < x_to_w - 1; ++log) {
        gflog[b] = (unsigned short) log;
        gfilog[log] = (unsigned short) b;
        b = b << 1;
        if ( b & x_to_w ) b = b ^ prim_poly;
    }

    return 0;
}

然后是乘法和除法操作，加法减法对于p=2的情况等同于XOR操作。

1
2
#define NW (1 << w)   /* In other words, NW equals 2 to the w-th power */

int mult(int a, int b)
{
    int sum_log;
    if (a == 0 || b == 0) return 0;
    sum_log = gflog[a] + gflog[b];
    if (sum_log >= NW-1) sum_log -= NW-1;
    return gfilog[sum_log];
}

int div(int a, int b)
{
    int diff_log;
    if (a == 0) return 0;
    if (b == 0) return -1;
    diff_log = gflog[a] - gflog[b];
    if (diff_log < 0) diff_log += NW-1;
    return gfilog[diff_log];
}