JErasure库提供一般的RS码和CRS码两种编码方式，是基于C的纠删码的类库。

JErasure库的相关学习

一、Galois Field
- （一）基本概念
- （二）GF(2w)GF(2^w)GF(2w)
二、CRS码
- （一）RS码简介
- （二）CRS Code简介
- - 2.2.1 Cauchy Matrix
  - 2.2.2 相关证明
三、JErasure库源码分析
- （一）Galois.c/Galois.h
- （二）JErasure.c/JErasure.h
- - 3.2.1 函数的参数介绍
  - 3.2.2 位矩阵编解码的时空优化方案
  - 3.2.3 无返回值(void)的函数
  - 3.2.4 返回整数(integer)的函数
  - 3.2.5 求点积（dot product）的函数
  - 3.2.6 矩阵转置操作函数
  - 3.2.7 一些基本的矩阵操作函数
  - 3.2.8 监测函数
- （三）cauchy.c/cauchy.h

一、Galois Field

（一）基本概念

1. 域（field）
域代表一组元素和元素间的四则运算的集合。其中，域内元素的四则运算必须满足封闭性，也就是说域内元素进行四则运算后得到的结果仍为域内元素。

2. 有限域（finite field）
有限域即为仅含有有限个元素的域。有限域又称作伽罗华域(Galois Field)。记为GFGFGF。例如，GF(7)GF(7)GF(7)的元素为0~6.

有限域中必须含有加法单位元和乘法单位元。加法单位元即域中的零元，即存在元素eee，使得域中任一元素aaa，有e+a=ae+a = ae+a=a。乘法单位元即域中的幺元，即存在元素e，使得域中任一元素aaa，有e∗a=ae*a = ae∗a=a。GF(p)GF(p)GF(p)的加法单位元和乘法单位元分别为0和1.

此外，有限域中的每个元素还应有对应的逆元。乘法逆元与我们常见的倒数类似。即存在元素b，使得b∗a=eb*a = eb∗a=e，则bbb为aaa的乘法逆元。其中eee为乘法单位元。加法逆元与我们常见的相反数类似。即存在元素bbb，使得b+a=eb+a = eb+a=e，则bbb为aaa的加法逆元。其中eee为加法单位元。

有限域中的减法和除法即为加法和乘法的逆过程。GF(p)GF(p)GF(p)的加法和乘法与实数域内的加法和乘法相似，只不过最后要模ppp。即GF(a+b)=(a+b)modp，GF(a×b)=(a×b)modpGF(a+b) = (a+b) \mod p， GF(a×b) = (a×b) \mod pGF(a+b)=(a+b)modp，GF(a×b)=(a×b)modp。

一般而言，ppp应为质数，因为如果ppp非质数，则可能有一些元素找不到乘法逆元。例如GF(9)GF(9)GF(9)中，元素3找不到GF(9)GF(9)GF(9)中的一个元素aaa使3×amodp=13\times a \mod p = 13×amodp=1，即3没有乘法逆元。

（二）GF(2w)GF(2^w)GF(2w)

1. GF(2w)GF(2^w)GF(2w)
GF(2w)GF(2^w)GF(2w)上的加法运算和乘法运算不使用一般的加法和乘法，而是使用多项式的计算。

2. 多项式计算
GF(2w)GF(2^w)GF(2w)上的多项式计算中，多项式的系数为1或0，即只能取GF(2)GF(2)GF(2)上的元素；当进行加法运算时，合并同类项时，不是一般的系数相加，而是进行异或计算，且加法和减法的操作一致。例如：
(x3+x2)+(x4+x2)=x4+x3(x^3+x^2)+(x^4+x^2)=x^4+x^3 (x3+x2)+(x4+x2)=x4+x3

GF(2w)GF(2^w)GF(2w)上的多项式的最高次不超过www。例如，在GF(23)GF(2^3)GF(23)中，f1(x)=x2+1f_1(x)=x^2+1f1(x)=x2+1为GF(23)GF(2^3)GF(23)的多项式，f2(x)=x3+x2+1f_2(x)=x^3+x^2+1f2(x)=x3+x2+1则不是GF(23)GF(2^3)GF(23)的多项式。

但是，多项式的加法计算后的结果一定仍为域中多项式，而对于多项式乘法，若直接按一般的多项式乘法计算的计算结果则可能不再是域中的多项式。与GF(p)GF(p)GF(p)类似，需要找到一个多项式对其取模，使之仍然为域中多项式。因此，提出本原多项式的概念。

3. 本原多项式/质多项式(primitive polynomial)

质多项式与质数的概念类似，质多项式即为不能表示成除幺元对应的多项式与该多项式以外的其他多项式的乘积，即不能再进行因式分解的多项式。例如，在GF(24)GF(2^4)GF(24)中，f(x)=x4+1f(x) = x^4 + 1f(x)=x4+1不是GF(24)GF(2^4)GF(24)上的质多项式，因为f(x)=(x2+1)×(x2+1)f(x) = (x^2 + 1) × (x^2 + 1)f(x)=(x2+1)×(x2+1)。g(x)=x4+x+1g(x) = x^4 + x + 1g(x)=x4+x+1是GF(24)GF(2^4)GF(24)上的质多项式，因为它不可再分解。

以GF(23)GF(2^3)GF(23)为例，指数小于3的多项式一共有8个，分别是0，1，x，x+1，x2，x2+1，x2+x，x2+x+10，1，x，x+1，x^2，x^2+1，x^2+x，x^2+x+10，1，x，x+1，x2，x2+1，x2+x，x2+x+1，一共有8个多项式。其对应的系数分别可以表示为000，001，010，011，100，101，110，111，若用十进制表示，则正好为0，1，2，3，4，5，6，7。这说明0~7与8个多项式之间有着必然的映射关系，每个多项式会对应一个值。对于GF(23)GF(2^3)GF(23)，取素多项式f(x)=x3+x+1f(x) = x^3+x+1f(x)=x3+x+1，则x+1x+1x+1的乘法逆元为x2+xx^2+xx2+x，因为((x+1)×(x2+x)mod(x3+x+1)=1)((x+1) × (x^2+x) mod (x^3+x+1) = 1)((x+1)×(x2+x)mod(x3+x+1)=1)，也就是说即使mod8\mod 8mod8不能构成一个有限域，但是由质多项式仍然可以为域中每个元素找到一个乘法逆元。
此时，多项式的乘法即可得到解决。例如，在GF(23)GF(2^3)GF(23)内：
x×(x2+1)=x×(x2+1)mod(x3+x+1)=(x3+x)mod(x3+x+1)=1x\times (x^2+1)=x\times (x^2+1)\mod (x^3+x+1)=(x^3+x)\mod (x^3+x+1)=1 x×(x2+1)=x×(x2+1)mod(x3+x+1)=(x3+x)mod(x3+x+1)=1

对于多项式除法，根据r(x)=q(x)t(x)+s(x)r(x) = q(x)t(x) + s(x)r(x)=q(x)t(x)+s(x)，其中r(x),q(x),t(x),s(x)r(x),q(x),t(x),s(x)r(x),q(x),t(x),s(x)均为GF(2w)GF(2^w)GF(2w)上的多项式，q(x)q(x)q(x)表示GF(2w)GF(2^w)GF(2w)上的质多项式。易知r(x)modq(x)=s(x)r(x) \mod q(x) = s(x)r(x)modq(x)=s(x)，且显然有dims(x)<dimq(x)dim\ s(x) < dim\ q(x)dim s(x)<dim q(x)。（dimf(x)dim\ f(x)dim f(x)表示的是f(x)f(x)f(x)的最高次的次数）。

根据上述理论，即可根据质多项式来得到GF(2w)GF(2^w)GF(2w)的所有元素。GF(2w)GF(2^w)GF(2w)中元素的生成过程如下：

1. 给定一个初始集合{0, 1, x};
2. 将这个集合的最后一个元素乘以x。如果得到的结果的度(dim) ≥ w，则将结果 mod q(x)，将结果加入集合中;
3. 重复2，直到集合中有2^w个元素。此时最后的一个元素乘以x再 mod q(x)的值必定为1.

举个栗子：
w=2w = 2w=2时，q(x)=x2+x+1q(x) = x^2 + x + 1q(x)=x2+x+1，初始集合为{0,1,x}\{0, 1, x\}{0,1,x}，

最后一个元素为xxx，乘以xxx后为x2x^2x2，dim(x2)=2dim\ (x^2) = 2dim (x2)=2，因此要对q(x)q(x)q(x)取模，得(x2)mod(x2+x+1)=x+1(x^2) \mod (x^2 + x + 1) = x + 1(x2)mod(x2+x+1)=x+1；此时集合更新为{0,1,x,x+1}\{0, 1, x, x+1\}{0,1,x,x+1}；

再将最后一个元素乘以xxx，为x2+xx^2 + xx2+x,仍需要对q(x)q(x)q(x)取模，得(x2+x)mod(x2+x+1)=1(x^2 + x) mod (x^2 + x + 1) = 1(x2+x)mod(x2+x+1)=1,结束上述过程。即GF(22)={0,1,x,x+1}GF(2^2) = \{0, 1, x, x+1\}GF(22)={0,1,x,x+1}。

也就是说，只需要找到对应GF(2w)GF(2^w)GF(2w)上的本原多项式q(x)q(x)q(x)即可通过xxx代换得到GF(2w)GF(2^w)GF(2w)，即：
GF(2w)=GF(2)[x]/q(x)GF(2^w) = GF(2)[x]/q(x) GF(2w)=GF(2)[x]/q(x)

常用的本原多项式有：

w = 4:   x^4 + x + 1
w = 8:   x^8 + x^4 + x^3 + x^2 + 1
w = 16:  x^16 + x^12 + x^3 + x + 1
w = 32:  x^32 + x^22 + x^2 + x + 1
w = 64:  x^64 + x^4 + x^3 + x + 1

例如，GF(24)GF(2^4)GF(24)的元素生成：(P(x)=x4+x+1P(x) = x^4 + x + 1P(x)=x4+x+1)

由上述计算结果，可以构建以下两个表（正表和反表）：

i	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
gflog[i]	-	0	1	4	2	8	5	10	3	14	9	7	6	13	11	12
gfilog[i]	1	2	4	8	3	6	12	11	5	10	7	14	15	13	9	-

其中，gflog[i]表（反表）是将二进制映射成多项式，例如将5映射成多项式，即gflog[5]=8gflog[5] = 8gflog[5]=8,即x8=5x^8 = 5x8=5；gfilog[i]表（正表）是将多项式映射成二进制(值),例如将5(即x5x^5x5)映射成值，即gfilog[5]=6gfilog[5] = 6gfilog[5]=6,即x5=6x^5 = 6x5=6.

对于乘法和除法，可以由如下方式进行：
xa×xb=xa+b;xa/xb=xa−b.x^a\times x^b=x^{a+b};\\ x^a/x^b=x^{a-b}. xa×xb=xa+b;xa/xb=xa−b.

因此，我们在计算可以a×ba×ba×b时，可以先得到a=xm,b=xna = x^m, b = x^na=xm,b=xn; 则a×b=xm+n,a/b=xm−na\times b = x^{m + n}, a/b = x^{m - n}a×b=xm+n,a/b=xm−n.

即：先将a，ba，ba，b查找反表（gflog[i]）得到m，nm，nm，n；计算m+nm + nm+n或者m−nm - nm−n，然后基于结果查找正表（gfilog[i]）得到乘法和除法的结果。需要注意的是，在计算m+nm + nm+n和m−nm - nm−n时，若结果大于等于15或小于0，则需要对15(2w−12^w - 12w−1)取模。因为0是单独存在的，在GF(2w)GF(2^w)GF(2w)中没有哪个元素的任意次方等于0.

再举个栗子：
7 × 9 = gfilog[gflog[7] + gflog[9]] = gfilog[10 + 14] = gfilog[9] = 10
(7×9=x10×x14=x10+14=x24mod15=x9=107 × 9 = x^{10} × x^{14} = x^{10 + 14} = x^{24 \mod 15} = x^9 = 107×9=x10×x14=x10+14=x24mod15=x9=10)
13 / 11 = gfilog[gflog[13] - gflog[11]] = gfilog[13 - 7] = gfilog[6] = 12
(13/11=x13/x7=x13−7=x6=1213 / 11 = x^{13} / x^7 = x^{13-7} = x^6 = 1213/11=x13/x7=x13−7=x6=12)

因此，构造正反表之后，可以大大简化乘法和除法的计算，可以直接通过加法和查表得到结果，而不需要再进行繁琐的取模等运算。

二、CRS码

（一）RS码简介

介绍CRS码之前先要弄明白RS码是什么。RS码就是Reed-Solomon码，是一种MDS码。编码时将k个数据块通过一定方式编码成m个校验块，我一般把它写成RS(k, m)，最多可以容忍m个数据块失效。解码时通过m个幸存块来恢复所有的数据块。那么RS码是通过什么方式来编码和解码的呢？
RS码通过构造生成矩阵来编码数据块。前k×k为单位阵，后m×k为生成矩阵。编码过程即为生成矩阵与数据块进行矩阵的乘法，即得到k个数据块和m个校验块。

当有不多于m个块失效时，为了恢复数据块，则取出k个幸存块对应的生成矩阵的行构成一个矩阵，对该矩阵的逆与对应的幸存块相乘即可恢复出数据块。
例如，数据块D2，D3D_2，D_3D2，D3和校验块P1P_1P1丢失，取出k个幸存块对应的生成矩阵的行，例如取出D1，D4，...，Dk，P2，P3D_1，D_4， ...，D_k，P_2，P_3D1，D4，...，Dk，P2，P3，一共kkk个数据块，有：

再将左边的子矩阵的逆乘到右侧，即有：

此时，所有的数据块就被恢复出来了。再要恢复校验块，很简单，再进行一次编码操作即可得到所有正确的校验块了。
然后问题来了，RS码的编解码步骤已经清楚了，现在就只有一个待解决的问题了，那就是得找到一个合适的生成矩阵，使得该生成矩阵的任意k行都线性无关，也就是说任意k行组成的矩阵的行列式不为0，此时即可保证任意k行组成的子矩阵必定可逆。因此，在RS码中，我们一般采用 范德蒙德矩阵（Vandermonde Matrix） 作为生成矩阵，此时该矩阵的任意k行组成的子矩阵一定可逆。
用范德蒙德矩阵进行编码的时间复杂度为O(mk)，解码的时间复杂度为O(k^3)。
在利用范德蒙德矩阵编码的时候，我们可以采用对数/反对数表的方法将乘法运算转换成加法运算，并且在迦罗华域中，加法运算转换成了XOR运算。
显然，这种方法要比直接复制副本存储的空间开销小得多，但是解码时花费的时间有点多了。所以，基于Reed-Solomon码，又有另一种纠删码被提出 – Cauchy Reed-Solomon Code。

（二）CRS Code简介

2.2.1 Cauchy Matrix

Cauchy矩阵是由两个交集为空的集合构成的矩阵。具体为：
令C = [c_ij]_m×\times×n，有集合X = {x1,x2,⋯,xmx_1, x_2, \cdots, x_mx1,x2,⋯,xm}，Y = {y1,y2,⋯,yny_1, y_2, \cdots, y_ny1,y2,⋯,yn}，且X∩Y = ∅。矩阵C中的元素cijc_{ij}cij满足
cij=1xi+yjc_{ij} = \frac{1}{x_i + y_j}cij=xi+yj1
则矩阵C为Cauchy矩阵。
Cauchy矩阵有一个很重要的性质，即Cauchy矩阵的任意一个子方阵必定可逆。且求逆的时间复杂度为O(n2n^2n2)，nnn为子方阵的维数。
例如，对于X={1,2}X = \{1, 2\}X={1,2}，Y={0,3,4,5,6}Y = \{0, 3, 4, 5, 6\}Y={0,3,4,5,6}，其编解码过程如下：

从编解码过程来看，柯西编解码最大的运算量是乘法和加法运算。柯西编解码为了降低乘法复杂度，采用了有限域上的元素都可以使用二进制矩阵表示的原理，将乘法运算和加法运算转换成了迦罗华域“与运算”和“XOR逻辑运算”，提高了编解码效率。
从数学的角度来看，在迦罗华有限域中，任何一个GF(2w2^w2w)域上的元素都可以映射到GF(2)二进制域，并且采用一个二进制矩阵的方式表示GF（2w2^w2w）中的元素。例如，GF（232^323）域中的元素可以表示成GF（2）域中的二进制矩阵：

其中，M(e)M(e)M(e)的第iii列为V(e∗2i−1)V(e * 2^{i-1})V(e∗2i−1).其中e∗2i−1e*2^{i-1}e∗2i−1为GF(23)GF(2^3)GF(23)上的乘法。通过以上规则，可以得到：
M(e1)∗V(e2)=V(e1e2)，M(e1)∗M(e2)=M(e1e2)M(e_1) * V(e_2) = V(e_1e_2)，\\ M(e_1) * M(e_2) = M(e_1e_2) \\ M(e1)∗V(e2)=V(e1e2)，M(e1)∗M(e2)=M(e1e2).
举个栗子：

根据上述理论，将Cauchy matrix的每个元素(eee)全部替换成M(e)M(e)M(e)，以此替换即可将计算从GF(2w)GF(2^w)GF(2w)上的计算变为GF(2)GF(2)GF(2)的计算(位运算)。为了可以使两矩阵能进行乘法计算，将每个数据块并没有划分为2w2^w2w大小，而是划分为www，即将每个数据块分为www块。例如，对前文所述的Cauchy matrix和GF(23)GF(2^3)GF(23)，有：

2.2.2 相关证明

上述结论可以证明。证明过程如下：
令P(x)P(x)P(x)为GF(2)[x]GF(2)[x]GF(2)[x]中的LLL次本原多项式。由前文介绍的Galois Field中的理论，可以得知GF(2L)GF(2^L)GF(2L)中的元素可以由f(x)=∑i=0L−1fixif(x) = \sum_{i = 0}^{L-1}f_ix^if(x)=∑i=0L−1fixi表出，其中fi∈GF(2)f_i\in GF(2)fi∈GF(2)。记列向量(f1,f2,⋯,fL−1)(f_1, f_2, \cdots, f_{L-1})(f1,f2,⋯,fL−1)为f(x)f(x)f(x)的系数向量（coefficient vector）。
对属于GF(2L)GF(2^L)GF(2L)的任意多项式f(x)f(x)f(x)，令τ(f)\tau(f)τ(f)为系数向量的矩阵，其中τ(f)\tau(f)τ(f)的第iii列为xi−1f(x)(modP(x))x^{i-1}f(x) \pmod {P(x)}xi−1f(x)(modP(x))的系数向量。
根据以上定义，有以下引理：
（1）τ(0)\tau(0)τ(0)是零矩阵；
（2）τ(1)\tau(1)τ(1)是单位矩阵；
（3）τ\tauτ是单射；
（4）τ(f)+τ(g)=τ(f+g)\tau(f) + \tau(g) = \tau(f+g)τ(f)+τ(g)=τ(f+g)；fff和ggg是GF(2L)GF(2^L)GF(2L)的多项式；
（5）τ(f)τ(g)=τ(fg)\tau(f)\tau(g) = \tau(fg)τ(f)τ(g)=τ(fg)；fff和ggg是GF(2L)GF(2^L)GF(2L)的多项式
其中(1)到(3)很容易证明，这里只证明(4)和(5):
Proof:Proof:Proof:
xi(f+g)=xif+xig(modP(x))x^i(f+g) = x^if + x^ig \pmod {P(x)}xi(f+g)=xif+xig(modP(x))
其中xi(f+g)x^i(f + g)xi(f+g)表示以τ(f+g)\tau(f + g)τ(f+g)的第iii列为系数的多项式，xifx^ifxif和xigx^igxig表示以τ(f)\tau(f)τ(f)和τ(g)\tau(g)τ(g)的第iii列为系数的多项式。显然，(4)得证。
令f(i)f^{(i)}f(i)表示τ(f)\tau(f)τ(f)的第iii列; 令{g0(j),g1(j),⋯,gL−1(j)}T\{g_0^{(j)}, g_1^{(j)}, \cdots, g_{L-1}^{(j)}\}^T{g0(j),g1(j),⋯,gL−1(j)}T为τ(g)\tau(g)τ(g)的第jjj列。则有：

∑i=0L−1gi(j)xi=xj−1gmodp(x)\sum_{i=0}^{L-1}g_i^{(j)}x^i = x^{j-1}g \mod p(x)i=0∑L−1gi(j)xi=xj−1gmodp(x)

(其中p(x)p(x)p(x)为GF(2L)GF(2^L)GF(2L)上的质多项式，g=g0(1)+g1(1)x+g2(1)x2+⋯+gL−1(1)xL−1g = g_0^{(1)}+g_1^{(1)}x + g_2^{(1)}x^2 + \cdots + g_{L-1}^{(1)}x^{L-1}g=g0(1)+g1(1)x+g2(1)x2+⋯+gL−1(1)xL−1)

因τ(f)τ(g)\tau(f)\tau(g)τ(f)τ(g)的第jjj列为：

[f(1)f(2)⋯f(L)][⋯g0(j)⋯⋯g1(j)⋯⋮⋮⋮⋯gL−1(j)⋯]\begin{bmatrix} f^{(1)} & f^{(2)} & \cdots & f^{(L)} \end{bmatrix} \begin{bmatrix} \cdots & g_0^{(j)} & \cdots \\ \cdots & g_1^{(j)} & \cdots \\ \vdots & \vdots & \vdots \\ \cdots & g_{L-1}^{(j)} & \cdots \\ \end{bmatrix} \\[10pt][f(1)f(2)⋯f(L)]⎣⎢⎢⎢⎢⎡⋯⋯⋮⋯g0(j)g1(j)⋮gL−1(j)⋯⋯⋮⋯⎦⎥⎥⎥⎥⎤

=g0(j)f(1)+g1(j)f(2)+⋯+gL−1(j)f(L)= g_0^{(j)}f^{(1)}+g_1^{(j)}f^{(2)}+ \cdots +g_{L-1}^{(j)}f^{(L)}=g0(j)f(1)+g1(j)f(2)+⋯+gL−1(j)f(L)

=∑i=0L−1gi(j)f(i+1)= \sum_{i=0}^{L-1}g_i^{(j)}f^{(i+1)}=∑i=0L−1gi(j)f(i+1)

这就是下面这个的多项式的系数向量：

∑i=0L−1gi(j)(xif)=f∑i=0L−1gi(j)xi=xj−1fgmodp(x)\sum_{i=0}^{L-1}g_i^{(j)}(x^if) = f\sum_{i=0}^{L-1}g_i^{(j)}x^i = x^{j-1}fg \mod p(x)i=0∑L−1gi(j)(xif)=fi=0∑L−1gi(j)xi=xj−1fgmodp(x)

即：τ(f)τ(g)\tau(f)\tau(g)τ(f)τ(g)的第j列与τ(fg)\tau(fg)τ(fg)的第j列相同
即：τ(f)τ(g)=τ(fg)\tau(f)\tau(g)=\tau(fg)τ(f)τ(g)=τ(fg)

也就是说，将GF(2w)GF(2^w)GF(2w)映射成τ(GF(2w))\tau(GF(2^w))τ(GF(2w))后，仍然有加法和乘法的不变性，即
M(e1)∗M(e2)=M(e1e2)M(e_1) * M(e_2) = M(e_1e_2)M(e1)∗M(e2)=M(e1e2)

若将M(e2)M(e_2)M(e2)的第一列提取出来，即为：
M(e1)∗V(e2)=V(e1e2)M(e_1) * V(e_2) = V(e_1e_2)M(e1)∗V(e2)=V(e1e2)

至此，前面的理论已完全证毕。

现在再来证明Cauchy matrix的求逆过程可以在O(n2)O(n^2)O(n2)的时间复杂度内完成。
Proof:Proof:Proof:
对于Cauchy矩阵的任意一个子方阵C=[cij]n×nC = [c_{ij}]_{n\times n}C=[cij]n×n:
[1x1+y11x1+y2⋯1x1+yn1x2+y11x2+y2⋯1x2+yn⋮⋮⋱⋮1xn+y11xn+y2⋯1xn+yn]\begin{bmatrix} \frac{1}{x_1+y_1} & \frac{1}{x_1+y_2} & \cdots & \frac{1}{x_1+y_n} \\[10pt] \frac{1}{x_2+y_1} & \frac{1}{x_2+y_2} & \cdots & \frac{1}{x_2+y_n} \\[10pt] \vdots & \vdots & \ddots & \vdots \\[10pt] \frac{1}{x_n+y_1} & \frac{1}{x_n+y_2} & \cdots & \frac{1}{x_n+y_n} \end{bmatrix}⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡x1+y11x2+y11⋮xn+y11x1+y21x2+y21⋮xn+y21⋯⋯⋱⋯x1+yn1x2+yn1⋮xn+yn1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤

其行列式det⁡C=∏i<j(xi−xj)∏i<j(yi−yj)∏i,j=1n(xi+yj)\det C = \frac{\prod_{i<j}(x_i-x_j)\prod_{i<j}(y_i-y_j)}{\prod_{i,j = 1}^n(x_i+y_j)}detC=∏i,j=1n(xi+yj)∏i<j(xi−xj)∏i<j(yi−yj)

令C−1=[dij]C^{-1} = [d_{ij}]C−1=[dij]，则有
dij=(−1)i+jdet⁡Cjidet⁡Cjd_{ij} = (-1)^{i+j}\frac{\det C_{ji}}{\det C_j}dij=(−1)i+jdetCjdetCji
其中cjic_{ji}cji表示删除了第jjj行和第iii列的矩阵CCC。令ak=∏i<k(xi−xk)∏k<j(xk−xj)a_k = \prod_{i<k}(x_i-x_k)\prod_{k<j}(x_k-x_j)ak=i<k∏(xi−xk)k<j∏(xk−xj)

bk=∏i<k(yi−yk)∏k<j(yj−yk)b_k = \prod_{i<k}(y_i-y_k)\prod_{k<j}(y_j-y_k)bk=i<k∏(yi−yk)k<j∏(yj−yk)

ck=∏i=1n(xk+yi)c_k = \prod_{i=1}^n(x_k+y_i)ck=i=1∏n(xk+yi)

fk=∏i=1n(yk+xi)f_k = \prod_{i=1}^n(y_k+x_i)fk=i=1∏n(yk+xi)

因此，有：
detC=∏k=1nakbk∏k=1nekfkdet C = \frac{\prod_{k=1}^na_kb_k}{\prod_{k=1}^ne_kf_k}detC=∏k=1nekfk∏k=1nakbk

det⁡Cji=det⁡Cejfiajbi(xj+yi)\det C_{ji} = \frac{\det C e_jf_i}{a_jb_i(x_j+y_i)}detCji=ajbi(xj+yi)detCejfi

因此：

dij=(−1)i+jejfiajbi(xj+yi)d_{ij} = (-1)^{i+j}\frac{e_jf_i}{a_jb_i(x_j+y_i)}dij=(−1)i+jajbi(xj+yi)ejfi

所以4n4n4n个未知量ak,bk,ck,fka_k, b_k, c_k, f_kak,bk,ck,fk可以在O(n2)O(n^2)O(n2)的时间复杂度内被解出，而dijd_{ij}dij则可由上述公式在O(1)O(1)O(1)的时间复杂度得到，因此，Cauchy矩阵的求逆过程可以在O(n2)O(n^2)O(n2)的时间复杂度内得到。

三、JErasure库源码分析

（一）Galois.c/Galois.h

这两个文件中定义的函数比较简单。首先：

static int prim_poly[33]; //对应的w的各个质多项式static int mult_type[33]; //对于不同的w采用不同的乘法计算static int nw[33];       //存放的是2^w的值，nwm1是2^w-1的值，即GF(2^w)的最大值static int *galois_log_tables[33];  //各个w(1~30)的反表
static int *galois_ilog_tables[33]; //各个w(1~30)的正表
static int *galois_mult_tables[33]; //各个w(1~30)的乘法表
static int *galois_div_tables[33];  //各个w(1~30)的除法表

以上是Galois.c中定义的一系列表，这些表将会在之后进行构造。

int galois_create_log_tables(int w)
//创建w的正反表，且将正表复制了两份，以便构建乘法表和除法表，然后将正表的初始位置(galois_log_tables[w][0]置为第二个副本的起始位置)，返回0表示构建成功，返回1表示失败int galois_logtable_multiply(int x, int y, int w)
//查正反表计算x*yint galois_logtable_divide(int x, int y, int w)
//查正反表计算x/yint galois_create_mult_tables(int w)
//构建大小为nw[w] * nw[w]的乘法表和除法表int galois_ilog(int value, int w)
//求value对应的二进制的值，即x^value对应的数值int galois_log(int value, int w)
//求value对应的多项式，即value对应的x^i的i值int galois_shift_multiply(int x, int y, int w)
//移位乘法，不需要任何表即可实现，但需要w次乘2，w次判断。
//对于x*y，先计算出y*2, y*2^2, y*2^3, ..., y*2^w-1的值，然后将x化为二进制，对x的第i位为1的，取出y*2^i进行异或得到结果。int galois_multtable_multiply(int x, int y, int w)
//查乘法表求x*yvoid galois_invert_binary_matrix(int *mat, int *inv, int rows)
//用初等变换法计算二进制矩阵mat的逆矩阵，mat的每一个元素表示一行，应为w，即w = row，若不可逆则停止运行程序int galois_inverse(int x, int w)
//若w = 23~32，用移位，否则直接用galois_single_multiply计算1/xint galois_shift_inverse(int x, int w)
//以移位的方式求y的倒数
//先由y得到mat2[w] = {y, y*2, y*2^2, ..., y*2^w-1}，通过galois_invert_binary_matrix函数得到mat2[w]的逆矩阵inv，inv[0]即为1/y的值int galois_single_multiply(int a, int b, int w)
//计算x*y，对于不同的w采用不同的乘法策略：
/********************************* w = 1~9: 构造乘法表，查乘法表* w = 10~22: 构建正反表，查正反表* w = 23~31: 用移位乘法* w = 32: 用裂乘********************************/int galois_single_divide(int a, int b, int w)
//计算a/b，对于不同的w采用不同的除法策略：
/********************************* w = 1~9: 构造除法表，查除法表* w = 10~22: 构建正反表，查正反表* w = 23~32: 先得到b的倒数(galois_inverse(b, w))，再用galois_single_multiply函数计算(分为w = 23~31和w = 32的两种策略)********************************/int galois_shift_divide(int x, int y, int w)
//利用galois_shift_inverse函数得到b的倒数，再利用移位乘法得到a/b的值int galois_multtable_divide(int x, int y, int w)
//查除法表得x/yint galois_create_split_w8_tables()
//创建裂乘表int galois_split_w8_multiply(int x, int y)
//裂乘

其中，裂乘的思路如下：
对于x×yx\times yx×y，以GF(232)GF(2^{32})GF(232)为例，xxx和yyy均为域内元素，则xxx和yyy可表示成：
x=x0+28×x1+22×8×x2+23×8×x3y=y0+28×y1+22×8×y2+23×8×y3x = x_0+2^8\times x_1+2^{2\times 8}\times x_2+2^{3\times 8}\times x_3 \\ y = y_0+2^8\times y_1+2^{2\times 8}\times y_2+2^{3\times 8}\times y_3 x=x0+28×x1+22×8×x2+23×8×x3y=y0+28×y1+22×8×y2+23×8×y3

其中xix_ixi和yiy_iyi均为unit8_t型数据，i∈{0,1,2,3}i \in \{0, 1, 2, 3\}i∈{0,1,2,3}.
则：
x×yx\times yx×y
=(x0+28×x1+22×8×x2+23×8×x3)×(y0+28×y1+22×8×y2+23×8×y3)= ( x_0+2^8\times x_1+2^{2\times 8}\times x_2+2^{3\times 8}\times x_3)\times (y_0+2^8\times y_1+2^{2\times 8}\times y_2+2^{3\times 8}\times y_3)=(x0+28×x1+22×8×x2+23×8×x3)×(y0+28×y1+22×8×y2+23×8×y3)
=x0y0=x_0y_0=x0y0
+28×(x0y1+x1y0)+\ 2^8\times (x_0y_1+x_1y_0)+ 28×(x0y1+x1y0)
+22×8×(x0y2+x1y1+x2y0)+\ 2^{2 \times 8}\times (x_0y_2+x_1y_1+x_2y_0)+ 22×8×(x0y2+x1y1+x2y0)
+23×8×(x0y3+x1y2+x2y1+x3y0)+\ 2^{3 \times 8}\times (x_0y_3+x_1y_2+x_2y_1+x_3y_0)+ 23×8×(x0y3+x1y2+x2y1+x3y0)
+24×8×(x1y3+x2y2+x3y1)+\ 2^{4 \times 8}\times (x_1y_3+x_2y_2+x_3y_1)+ 24×8×(x1y3+x2y2+x3y1)
+25×8×(x2y3+x3y2)+\ 2^{5 \times 8}\times (x_2y_3+x_3y_2)+ 25×8×(x2y3+x3y2)
+26×8×(x3y3)+\ 2^{6 \times 8}\times (x_3y_3)+ 26×8×(x3y3)

因此，设α\alphaα为unit32_t类型数中只有低8位的任意数字，那么根据上述x×yx\times yx×y的拆分结果可以归结为α×α,α×(α<<8),α×(α<<16),α×(α<<24),\alpha \times \alpha, \alpha \times (\alpha << 8), \alpha \times (\alpha<<16),\alpha \times (\alpha<<24),α×α,α×(α<<8),α×(α<<16),α×(α<<24), (α<<8)×(α<<24),(α<<16)×(α<<24),(α<<24)×(α<<24)(\alpha << 8)\times (\alpha<<24),(\alpha << 16)\times (\alpha<<24),(\alpha << 24)\times (\alpha<<24)(α<<8)×(α<<24),(α<<16)×(α<<24),(α<<24)×(α<<24)
一共7种运算。每一种运算可以用一个乘法表表示，每个乘法表大小为28×28×4Bytes2^8\times 2^8\times 4\ Bytes28×28×4 Bytes.

以上为Galois field上的一些基本运算。

void galois_region_xor(char *r1,         /* Region 1 */char *r2,         /* Region 2 */char *r3,         /* Sum region (r3 = r1 ^ r2) */int nbytes);      /* Number of bytes in region */void galois_w08_region_multiply(char *region,       /* Region to multiply */int multby,       /* Number to multiply by */int nbytes,       /* Number of bytes in region */char *r2,
/* If r2 != NULL, products go here.Otherwise region is overwritten */int add);
/* If (r2 != NULL && add) the produce is XOR'd with r2 */void galois_w16_region_multiply(char *region,       /* Region to multiply */int multby,       /* Number to multiply by */int nbytes,       /* Number of bytes in region */char *r2,
/* If r2 != NULL, products go here. Otherwise region is overwritten */int add);
/* If (r2 != NULL && add) the produce is XOR'd with r2 */void galois_w32_region_multiply(char *region,       /* Region to multiply */int multby,       /* Number to multiply by */int nbytes,       /* Number of bytes in region */char *r2,
/* If r2 != NULL, products go here. Otherwise region is overwritten */int add);
/* If (r2 != NULL && add) the produce is XOR'd with r2 */// 上述三个函数均为将region中各个元素乘以multby
// 1. r2 == NULL,则将结果写入region，结束
// 2. r2 != NULL 且 add is false,则将结果写入r2
// 3. r2 != NULL 且 add is true,将结果与r2异或再写入r2

以上四个函数为域内乘法，第一个表示将两个等大的数组r1r_1r1和r2r_2r2按元素进行异或，将结果置入r3r_3r3中。后面3个函数主要用来求点积，具体的实现步骤已经标注在上面的注释中了，点积的概念将在之后介绍到。

（二）JErasure.c/JErasure.h

3.2.1 函数的参数介绍

在分析函数功能之前，先弄明白这两个文件中定义的函数的参数。

KKK = 数据块数量

MMM = 校验块$(coding devices)数量

WWW = 字长

data_ptrs[][]data\_ptrs[][]data_ptrs[][] = 包含kkk个指向大小为size bytesbytesbytes的数据的指针的数组，即数据块指针，sizesizesize必须是sizeof(long)的整数倍

coding_ptrs[][]coding\_ptrs[][]coding_ptrs[][] = 指向编码数据（大小为字节）的m个指针数组，即校验块指针

packetsizepacketsizepacketsize = 使用bitmatrix编码的编码块(coding block)大小

Matrix[]Matrix[]Matrix[] = 包含k×mk\times mk×m个整数(integer)的数组，即编码矩阵，其中iii行jjj列元素为matrix[i×k+j]matrix[i\times k+j]matrix[i×k+j]

Bitmatrix[]Bitmatrix[]Bitmatrix[] = 包含kw×mwkw \times mwkw×mw个整数(integer)的数组，其中iii行jjj列元素为bitmatrix[i×kw+j]bitmatrix[i\times kw+j]bitmatrix[i×kw+j]

Erasures[]Erasures[]Erasures[] = 已删除的设备(devices)(块)的ididid的数组，用于记录哪些块丢失
id∈[0,k+m−1]id∈[0, k+m-1]id∈[0,k+m−1]
id∈[0,k−1]：数据块idid∈[0, k-1]：数据块idid∈[0,k−1]：数据块id
id∈[k,k+m−1]：编码块(校验块)idid∈[k, k+m-1]：编码块(校验块)idid∈[k,k+m−1]：编码块(校验块)id
例如：erasures[0]=0;例如：erasures[0] = 0;例如：erasures[0]=0; //第000个块丢失
erasures[1]=3;erasures[1] = 3;erasures[1]=3; //第333个块丢失
erasures[2]=−1;erasures[2] = -1;erasures[2]=−1; //结束标志

Operation[op][]Operation[op][]Operation[op][] = 包含5个整数的数组(五元组)
4:进行操作−0表示复制，1表示异或4 : 进行操作 - 0表示复制，1表示异或4:进行操作−0表示复制，1表示异或
0:源设备(device)[0,k+m−1](operation[op][0]=−1表示结束)0 : 源设备(device)[0,k+m-1](operation[op][0] = -1表示结束)0:源设备(device)[0,k+m−1](operation[op][0]=−1表示结束)
1:源数据包(packet)[0,w−1]1 : 源数据包(packet)[0,w-1]1:源数据包(packet)[0,w−1]
2:目标设备[0,k+m−1]2 : 目标设备[0,k+m-1]2:目标设备[0,k+m−1]
3:目标数据包[0,w−1]3 : 目标数据包[0,w-1]3:目标数据包[0,w−1]

我这里只介绍各个接口的实现和作用，各位最好还是对照源码来理解各个参数的含义。

3.2.2 位矩阵编解码的时空优化方案

因为位矩阵的大小变成了kw×mwkw\times mwkw×mw，比之前的矩阵大了w2−1w^2-1w2−1倍，空间利用率极低。而且在编码校验块的时候，中间可能会有大量的XOR操作是重复的，例如P1,1=D1,1⊗D1,2;P1,2=D1,1⊗D1,2⊗D2,2;P2,1=D1,1⊗D2,2⊗D3,1P_{1,1}=D_{1,1} \otimes D_{1,2};\\ P_{1,2}=D_{1,1} \otimes D_{1,2} \otimes D_{2,2};\\ P_{2,1}=D_{1,1} \otimes D_{2,2} \otimes D_{3,1}P1,1=D1,1⊗D1,2;P1,2=D1,1⊗D1,2⊗D2,2;P2,1=D1,1⊗D2,2⊗D3,1

因此，可以由P1,2=P1,1⊗D2,2P2,1=P1,2⊗D2,2⊗D3,1P_{1,2}=P_{1,1} \otimes D_{2,2} \\P_{2,1}=P_{1,2}\otimes D_{2,2}\otimes D_{3,1}P1,2=P1,1⊗D2,2P2,1=P1,2⊗D2,2⊗D3,1 来减少XOR次数，一般而言，kwkwkw和mwmwmw较大，所以以这种方式来减少XOR次数是可以减少编解码时间的。

刚刚那个栗子只减少了1次XOR操作，看起来8太行。现在来举个好一点的栗子：

由上面的图上所示，按照位矩阵的XOR操作得到校验块，可以得到：

但是很显然，中间有大量的XOR操作可以重复利用。如果将其重复的结果进行分类，则可以有以下结果：

很显然，这次是大大的减少所需的XOR操作。

以这种思路，JErasure库给出了两种方式来分别实现提高空间利用率和进一步减少XOR的操作数：

int **jerasure_dumb_bitmatrix_to_schedule(int k, int m, int w, int *bitmatrix);
//将bitmatrix转换为一系列操作。简单粗暴。
//返回一个operations[][]数组，具体实现方式为：
//若bitmatrix[index] == 1，则计算得到当前bitmatrix位置对应的数据块位置，如果在该行是第一个bitmatrix[index] == 1，则置operation为0(copy)，否则为1(XOR)int **jerasure_smart_bitmatrix_to_schedule(int k, int m, int w, int *bitmatrix);
// 采用与上一函数类似的方法，不过使用之前的结果来计算新的点积。

这两个函数都是将位矩阵(bitmatrix)转化为了调度(schedule)，既减少空间开销，也减少了XOR次数。
先来分析第一个函数：(dumb)
其实就是在模拟位矩阵与数据包的矩阵乘法过程。例如：

例如: C11=D11⊗D21⊗D22⊗D33⊗D41⊗D42⊗D43⊗D52C_{11}=D_{11} \otimes D_{21} \otimes D_{22} \otimes D_{33}\otimes D_{41} \otimes D_{42} \otimes D_{43} \otimes D_{52}C11=D11⊗D21⊗D22⊗D33⊗D41⊗D42⊗D43⊗D52
jerasure_dumb_bitmatrix_to_schedulejerasure\_dumb\_bitmatrix\_to\_schedulejerasure_dumb_bitmatrix_to_schedule即模拟D11,⋯,D52D_{11}, \cdots, D_{52}D11,⋯,D52的操作，且每个operation元素对应了一个源块和目标块的位置，以及要进行的操作。
例如，D11D_{11}D11对应的源块即为(1,1)(1,1)(1,1)，目标块C11C_{11}C11即为(1,1)(1,1)(1,1)，因为D11D_{11}D11在第一个位置上，所以操作为复制(copy)，即为(0)(0)(0)，因此，D11D_{11}D11对应的五元组为(1,1,1,1,0)(1,1,1,1,0)(1,1,1,1,0).
再例如，D21D_{21}D21对应的源块即为(2,1)(2,1)(2,1)，目标块C11C_{11}C11即为(1,1)(1,1)(1,1)，因为D21D_{21}D21不在第一个位置上，也就是说之前已经开始计算异或结果了，所以此时操作为异或(XOR)，即为(1)(1)(1)，代表XOR D21D_{21}D21，因此，D11D_{11}D11对应的五元组为(2,1,1,1,1)(2,1,1,1,1)(2,1,1,1,1).

现在来分析第二个函数：(smart)
前面那个函数只是开胃菜，这个函数才是对位运算的大大优化。基本思路与上一函数类似，但是这个会使用之前计算过的XOR结果(点积)来编码新的校验块。
下面是这种中间结果复用的实现方式：

先定义一些记号，以便后续说明：
from[]含有r个元素，初始置−1from[]含有r个元素，初始置-1from[]含有r个元素，初始置−1

diff[]含r个元素，diff[i]表示Mi中1的数量diff[]含r个元素，diff[i]表示Mi中1的数量diff[]含r个元素，diff[i]表示Mi中1的数量

flink[]和blink[]用于定位top，top用于找当前待计算的校验块flink[]和 blink[]用于定位top，top用于找当前待计算的校验块flink[]和blink[]用于定位top，top用于找当前待计算的校验块

bestdiff记录当前最小的diff[i]的值bestdiff记录当前最小的diff[i]的值bestdiff记录当前最小的diff[i]的值

bestrow记录当前bestdiff对应的diff[i]的i的值bestrow记录当前bestdiff对应的diff[i]的 i的值bestrow记录当前bestdiff对应的diff[i]的i的值

实现步骤：
假定MMM为bitmatrixbitmatrixbitmatrix，VVV为输入的数据块，UUU为生成的校验块，notdonenotdonenotdone代表待生成的校验块ididid的集合，初始为000 ~ mw−1mw-1mw−1，则步骤如下:

1. 选择未计算的校验块的packet且diff[i]最小的 i，即找到i∈notdone且diff[i]最小
2. 若from[i] = -1, 则令operation[op]的元素为M[i, j] = 1的所有对应的数据块;若from[i] ≠ -1, 则令operation[op]的第一个元素为校验块from[i]，之后的元素为第from[i]行与i行相异的bitmatrix中的元素对应的数据块
3. 从notdone中删除i
4. 对于notdone中的任一元素j，计算第j行与第i行中相异的元素的数量c，令x = c + 1，若x < diff[j]，则diff[j] = x，from[j] = i

也就是说，在某一次计算中，如果from[i]!=−1from[i] != -1from[i]!=−1，则说明第iii个校验块的计算中用前面的计算过的结果来计算第i个校验块比直接进行矩阵乘法用到的XOR数要少，如果from[i]==−1from[i] == -1from[i]==−1，则说明用前面的结果计算该校验块不如直接计算来的快，就直接从bitmatrixbitmatrixbitmatrix中进行异或运算。

效率分析：
例如，在解码过程中，要计算得到D0D_0D0和D1D_1D1，通过一般的矩阵乘法计算，一共需要124次XORs。但是，如果按照上述方案：

找到含有最少的1的一行(记为第iii行)，即计算D1,3D_{1,3}D1,3对应的一行，一共有8个1，因此共7次XORs。
计算其他9行(记为第jjj行)与D1,3D_{1,3}D1,3对应的bitmatrixbitmatrixbitmatrix的行中相异元素的数量，如果第jjj行中1的数量小于刚才计算出的数量+1，说明由第jjj行计算块得到结果需要的XOR数量多于由第jjj行计算块得到结果需要的XOR数量，此时将from[j]from[j]from[j]置为iii，即第jjj行的XOR数量由第iii行的结果来得到比当前方法更快。
重复之，直到计算得到了所有的块。

根据上述理论分析，可以得出以下计算步骤：

此时一共需要46次XORs操作。相比之前的124次XORs，这种方法得到了大幅优化。

3.2.3 无返回值(void)的函数

接下来是几个无返回值的函数，比较简单，所以拿出来放在一块儿分析了。

void jerasure_do_parity(int k, char **data_ptrs, char *parity_ptr, int size);
// 由data_ptrs直接进行XOR操作生成校验块，结果放入parity_ptr中void jerasure_matrix_encode(int k, int m, int w, int *matrix,char **data_ptrs, char **coding_ptrs, int size);
// 按常规方式编码，w必须是8，16，32void jerasure_bitmatrix_encode(int k, int m, int w, int *bitmatrix,char **data_ptrs, char **coding_ptrs, int size, int packetsize);
// 按常规方式对位矩阵编码，w ∈ [1, 32] and w ∈ N*void jerasure_schedule_encode(int k, int m, int w, int **schedule,char **data_ptrs, char **coding_ptrs, int size, int packetsize);
// 使用jerasure_smart_bitmatrix_to_schedule()或者jerasure_dumb_bitmatrix_to_schedule()进行编码，具体方案由schedule[][]确定

没有什么难懂的地方，我根据自己的理解加了两句注释，大家看看代码也很容易就看懂了。

3.2.4 返回整数(integer)的函数

这种返回整数的函数都是值返回000和−1-1−1的，返回000说明这个函数完整的进行了操作，返回成功的标志；返回−1-1−1说明参数有误，中途中断，返回失败的标志。

int jerasure_matrix_decode(int k, int m, int w, int *matrix, int row_k_ones, int *erasures,char **data_ptrs, char **coding_ptrs, int size);

kkk表示数据块，mmm表示校验块，matrixmatrixmatrix表示生成矩阵，row_k_onesrow\_k\_onesrow_k_ones用于记录编码的第一行是否为全1，用于优化，data_ptrsdata\_ptrsdata_ptrs表示指向数据块的指针，coding_ptrscoding\_ptrscoding_ptrs表示指向校验块的指针，erasureserasureserasures记录哪些块失效，sizesizesize表示块大小。
这个函数表示创建GF(2w)GF(2^w)GF(2w)上的解码矩阵进行解码，且最终会丢弃解码矩阵。

int jerasure_bitmatrix_decode(int k, int m, int w, int *bitmatrix, int row_k_ones, int *erasures,char **data_ptrs, char **coding_ptrs, int size, int packetsize);

表示创建GF(2)GF(2)GF(2)上的矩阵（位矩阵）进行解码，参数同上。packetsizepacketsizepacketsize是使用位矩阵进行编解码时一个packetpacketpacket的大小。

int jerasure_schedule_decode_lazy(int k, int m, int w, int *bitmatrix, int *erasures,char **data_ptrs, char **coding_ptrs, int size, int packetsize, int smart);

若smart==1smart == 1smart==1，则采用jerasure_smart_bitmatrix_to_schedule()jerasure\_smart\_bitmatrix\_to\_schedule()jerasure_smart_bitmatrix_to_schedule()来创建调度（schedule），否则采用jerasure_smart_bitmatrix_to_schedule()jerasure\_smart\_bitmatrix\_to\_schedule()jerasure_smart_bitmatrix_to_schedule()。

int jerasure_make_decoding_matrix(int k, int m, int w, int *matrix, int *erased, int *decoding_matrix, int *dm_ids);

创建解码矩阵，但是不解码，解码矩阵存放在decoding_matrixdecoding\_matrixdecoding_matrix中，dm_idsdm\_idsdm_ids含有kkk个整数，即用这kkk个存活块来构建解码矩阵，最终生成的decoding_matrixdecoding\_matrixdecoding_matrix含有k×kk\times kk×k个整数。

int jerasure_make_decoding_bitmatrix(int k, int m, int w, int *matrix, int *erased, int *decoding_matrix, int *dm_ids);

同上，创建解码位矩阵，但是不解码，解码矩阵存放在decoding_matrixdecoding\_matrixdecoding_matrix中，dm_idsdm\_idsdm_ids含有kwkwkw个整数，即用这kwkwkw个存活块来构建解码矩阵，最终生成的decoding_matrixdecoding\_matrixdecoding_matrix含有kw×kwkw\times kwkw×kw个整数。

int *jerasure_erasures_to_erased(int k, int m, int *erasures);

将erasures[]erasures[]erasures[]转换成erased[]erased[]erased[]返回，其中erased[]erased[]erased[]标识了哪些位置上的 packetpacketpacket不可用。例如erased[2]=1erased[2] = 1erased[2]=1说明第222个packetpacketpacket不可用（从第000个packetpacketpacket开始计数）。举个栗子：
若erasures[]={0,2,3}erasures[] = \{0,2,3\}erasures[]={0,2,3} 则erased[]={1,0,1,1}erased[] = \{1,0,1,1\} erased[]={1,0,1,1}

3.2.5 求点积（dot product）的函数

void jerasure_matrix_dotprod(int k, int w, int *matrix_row,int *src_ids, int dest_id,char **data_ptrs, char **coding_ptrs, int size);

简单来说，就是将生成矩阵（编码矩阵/解码矩阵）的一行元素与数据块/幸存块相乘（求点积）。

数据块/幸存块的id存放在src_idssrc\_idssrc_ids中，例如：

D1,1D_{1,1}D1,1的id为000，D4,3D_{4,3}D4,3的id为121212，C1,1C_{1,1}C1,1的id为161616等。

*注意：此处为了简单起见，给出的例子全部是以111开始计算id的，而一般情况下，应该从000开始。例如将上述D1→5D_{1\to5}D1→5和C1→2C_{1\to2}C1→2换成D0→4D_{0\to4}D0→4和C0→1C_{0\to1}C0→1，则D0,0D_{0,0}D0,0的id是000，D4,3D_{4,3}D4,3的id是111111等。

编码/解码生成的块的id存放在dest_iddest\_iddest_id中。若src_ids=nullsrc\_ids = nullsrc_ids=null，则说明是在编码过程中，直接利用data_ptrsdata\_ptrsdata_ptrs计算点积即可，若src_ids≠nullsrc\_ids\ne nullsrc_ids=null，则说明是在解码过程中，则需要从data_ptrsdata\_ptrsdata_ptrs和coding_ptrscoding\_ptrscoding_ptrs中得到对应的数据块/幸存块再计算点积。

该函数的对乘法计算进行了稍微的优化工作。因为点击计算需要GF(2w)GF(2^w)GF(2w)的乘法和加法，所以它会首先计算matrixmatrixmatrix中为000或111对应的块，然后再计算其他块，以简化计算的复杂度。例如：
[⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯13021037⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯][D1D2D3D4D5D6D7D8]=[⋯⋯⋯Ci⋯⋯⋯]\begin{bmatrix} \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 1 & 3 & 0 & 2 & 1 & 0 & 3 & 7 \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \end{bmatrix} \begin{bmatrix} D_1 \\[8pt] D_2 \\[8pt] D_3 \\[8pt] D_4 \\[8pt] D_5 \\[8pt] D_6 \\[8pt] D_7 \\[8pt] D_8 \end{bmatrix}= \begin{bmatrix} \cdots \\[8pt] \cdots \\[8pt] \cdots \\[8pt] C_i \\ \cdots \\[8pt] \cdots \\[8pt] \cdots \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡⋯⋯⋯1⋯⋯⋯⋯⋯⋯3⋯⋯⋯⋯⋯⋯0⋯⋯⋯⋯⋯⋯2⋯⋯⋯⋯⋯⋯1⋯⋯⋯⋯⋯⋯0⋯⋯⋯⋯⋯⋯3⋯⋯⋯⋯⋯⋯7⋯⋯⋯⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡D1D2D3D4D5D6D7D8⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡⋯⋯⋯Ci⋯⋯⋯⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤

一般情况下，计算得到Ci=1D1⊕3D2⊕0D3⊕2D4⊕1D5⊕0D6⊕3D7⊕7D8C_i=1D_1\oplus 3D_2\oplus 0D_3\oplus 2D_4\oplus 1D_5\oplus 0D_6\oplus 3D_7\oplus 7D_8Ci=1D1⊕3D2⊕0D3⊕2D4⊕1D5⊕0D6⊕3D7⊕7D8

该函数则先计算0和1对应的块，即D1,D3,D5,D6D_1,D_3,D_5,D_6D1,D3,D5,D6。D3,D6D_3,D_6D3,D6对应的矩阵元素为0，所以可以忽略(α⊕0=α\alpha \oplus 0=\alphaα⊕0=α)，即先得到D1⊕D5D_1\oplus D_5D1⊕D5。接下来再利用galois_w_region_multiplygalois\_w\_region\_multiplygalois_w_region_multiply计算剩余的非1的点积，然后将二者异或即得到CiC_iCi。

void jerasure_bitmatrix_dotprod(int k, int w, int *bitmatrix_row,int *src_ids, int dest_id,char **data_ptrs, char **coding_ptrs, int size, int packetsize);

含义同上一函数，方法为直接计算。从这里也可看出，GF(2)GF(2)GF(2)上的点积计算比GF(2w)GF(2^w)GF(2w)方便得多。

3.2.6 矩阵转置操作函数

int jerasure_invert_matrix(int *mat, int *inv, int rows, int w);

利用初等变换法求matmatmat的逆矩阵，将结果放入invinvinv中。大小为rows×rowsrows\times rowsrows×rows。

int jerasure_invert_bitmatrix(int *mat, int *inv, int rows);

利用初等变换法求matmatmat的逆矩阵，将结果放入invinvinv中。大小为rows×rowsrows\times rowsrows×rows。此时的w=1w=1w=1，因此也可以用jerasure_invert_matrixjerasure\_invert\_matrixjerasure_invert_matrix求逆，但是该函数的时间复杂度为O(n2)O(n^2)O(n2)，用上一函数令w=1w=1w=1的时间复杂度为O(n3)O(n^3)O(n3)。

int jerasure_invertible_matrix(int *mat, int rows, int w);

判断matmatmat是否可逆，不会求逆。

int jerasure_invertible_bitmatrix(int *mat, int rows);

同上，判断matmatmat是否可逆，不会求逆。

3.2.7 一些基本的矩阵操作函数

void jerasure_print_matrix(int *matrix, int rows, int cols, int w);

输出matrixmatrixmatrix，根据www的大小确定列间距。

void jerasure_print_bitmatrix(int *matrix, int rows, int cols, int w);

输出matrixmatrixmatrix，每www个字符后插入一个空格，每www行后插入一个空白行。

例如：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>void jerasure_print_matrix(int *m, int rows, int cols, int w)
{int i, j;int fw;char s[30];unsigned int w2;if (w == 32){fw = 10;}else{w2 = (1 << w);sprintf(s, "%u", w2 - 1);fw = strlen(s);}for (i = 0; i < rows; i++){for (j = 0; j < cols; j++){if (j != 0)printf(" ");printf("%*u", fw, m[i * cols + j]);}printf("\n");}
}void jerasure_print_bitmatrix(int *m, int rows, int cols, int w)
{int i, j;for (i = 0; i < rows; i++){if (i != 0 && i % w == 0)printf("\n");for (j = 0; j < cols; j++){if (j != 0 && j % w == 0)printf(" ");printf("%d", m[i * cols + j]);}printf("\n");}
}int main()
{int mat[] = {1, 5, 2, 7, 4, 5, 1, 3, 4, 7};int bitmat[6*15] = {   1, 0, 0,    1, 1, 0,    0, 0, 1,    1, 1, 1,    0, 1, 0,0, 1, 0,    0, 0, 1,    1, 0, 1,    1, 0, 0,    0, 1, 1,0, 0, 1,    1, 0, 0,    0, 1, 0,    1, 1, 0,    1, 0, 1,1, 1, 0,    1, 0, 0,    1, 0, 1,    0, 1, 0,    1, 1, 1,0, 0, 1,    0, 1, 0,    1, 1, 1,    0, 1, 1,    1, 0, 0,1, 0, 0,    0, 0, 1,    0, 1, 1,    1, 0, 1,    1, 1, 0};jerasure_print_matrix(mat, 2, 5, 3);putchar(10);jerasure_print_bitmatrix(bitmat, 2*3, 5*3, 3);system("pause");return 0;
}

输出结果为：

int *jerasure_matrix_multiply(int *m1, int *m2, int r1, int c1, int r2, int c2, int w);

计算矩阵乘法。m1m_1m1为r1×c1r_1\times c_1r1×c1，m2m_2m2为r2×c2r_2\times c_2r2×c2，计算m1×m2m_1\times m_2m1×m2，显然c1=r2c_1=r_2c1=r2。

3.2.8 监测函数

为评估编解码性能，JErasure库提供了一个监测函数。

void jerasure_get_stats(double *fill_in);

fill_infill\_infill_in是含有333个doubledoubledouble类型的向量，分别评价三个指标，其含义分别为：

使用galois_region_xor()galois\_region\_xor()galois_region_xor()进行了XOR运算的字节数；
使用galois_w08_region_multiply()galois\_w08\_region\_multiply()galois_w08_region_multiply(), galois_w16_region_multiply()galois\_w16\_region\_multiply()galois_w16_region_multiply() or galois_w32_region_multiply()galois\_w32\_region\_multiply()galois_w32_region_multiply()来乘以常数的字节数，即regionregionregion的长度；
使用memcpy()memcpy()memcpy()复制的字节数。
调用该函数之后，上述三个指标将会被清空。具体如下：

void jerasure_get_stats(double *fill_in)
{fill_in[0] = jerasure_total_xor_bytes;fill_in[1] = jerasure_total_gf_bytes;fill_in[2] = jerasure_total_memcpy_bytes;jerasure_total_xor_bytes = 0;jerasure_total_gf_bytes = 0;jerasure_total_memcpy_bytes = 0;
}

（三）cauchy.c/cauchy.h

JErasure库相关介绍相关推荐

axios队列 vue_(十三 )Vue 封装axios(四种请求)及相关介绍
Vue 封装axios(四种请求)及相关介绍首先axios是基于promise的http库 promise是什么? 1.主要用于异步计算 2.可以将异步操作队列化,按照期望的顺序执行,返回符合预期的 ...
python映射类型-python映射类型的相关介绍
映射类型是一类可迭代的键-值数据项的组合,提供了存取数据项及其键和值的方法,在python3中,支持两种无序的映射类型:内置的dict和标准库中的collections.defaultdict类型. ...
Android O 前期预研之二：HIDL相关介绍
在上一篇博客里,大致介绍了下Android O 中treble计划的一些背景与相关基本架构,这一篇中跟大家一起来探讨下HIDL相关的内容. Android HAL类型在此之前的ANDROID版本当 ...
Nginx工作原理及相关介绍
Nginx工作原理及相关介绍一.Nginx工作原理与模块介绍 1.Nginx基本工作原理 NGINX以高性能的负载均衡器,缓存,和web服务器闻名.Nginx由内核和模块组成,其中,内核的设计非常微 ...
数据分析与挖掘中常用Python库的介绍与实践案例
数据分析与挖掘中常用Python库的介绍与实践案例一.Python介绍现在python一词对我们来说并不陌生,尤其是在学术圈,它的影响力远超其它任何一种编程语言, 作为一门简单易学且功能强大的编程 ...
matlab图片导出无失真库export_fig介绍（半透明效果）
matlab图片导出无失真半透明等功能的库export_fig介绍首先,感谢export_fig的作者Yair Altman为相关方面做了很多介绍,本文主要结合新版本matlab,对作者的内容进行搬 ...
阿里MNN推理框架相关介绍
一.参考资料 MNN官网中文文档-语雀欢迎使用MNN文档 - MNN-Doc 2.1.1 documentation) 英文文档 MNN知识库 MNN 官方仓库二.相关介绍 1. MNN简介 M ...
昇腾Ascend处理器相关介绍
一.参考资料 modelzoo wiki 解密昇腾AI处理器–Ascend310简介 AI芯片:华为Ascend(昇腾)910结构分析解密昇腾AI处理器–DaVinci架构(计算单元) 二.相关介绍 ...

JErasure库相关介绍