前言
Reed-Solomon是一种应用广泛的纠错码。 ( n , k ) (n,k) (n,k)的RS码可以容忍 n − k + 1 2 \frac{n-k+1}{2} 2n−k+1个错误,编码 k k k个信息,需要 n n n个码字。也就是说,当错误的码字小于 n − k + 1 2 \frac{n-k+1}{2} 2n−k+1时,可以完全恢复 k k k个信息。
编码
假设
q
q
q是一个素数,
F
q
\mathbb{F}_q
Fq是有
q
q
q个元素的有限域,
n
,
k
,
d
n,k,d
n,k,d是整数,满足
1
≤
k
<
n
≤
q
,
d
=
n
−
k
+
1
1\leq k<n\leq q, d=n-k+1
1≤k<n≤q,d=n−k+1。
有
n
n
n个不同的
F
q
\mathbb{F}_q
Fq中的元素
a
1
,
a
2
,
⋯
,
a
n
a_1,a_2,\cdots,a_n
a1,a2,⋯,an,我们要加密
k
k
k个信息
m
1
,
m
2
,
⋯
,
m
k
m_1,m_2,\cdots,m_k
m1,m2,⋯,mk,其中
m
i
m_i
mi是
F
q
\mathbb{F}_q
Fq中的元素。
那么
k
k
k个信息作为系数,可以得到一个
k
−
1
k-1
k−1次的多项式
f
(
x
)
=
m
1
+
m
2
x
+
⋯
+
m
k
x
k
−
1
.
f(x)=m_1+m_2x+\cdots+m_kx^{k-1}.
f(x)=m1+m2x+⋯+mkxk−1.
计算码字
c
i
=
f
(
a
i
)
∈
F
q
.
c_i=f(a_i) \in \mathbb{F}_q.
ci=f(ai)∈Fq.
则我们得到的编码为
c
=
(
c
1
,
c
2
,
⋯
,
c
n
)
.
\mathbf{c}=(c_1,c_2,\cdots,c_n).
c=(c1,c2,⋯,cn).
多项式除法解码
首先计算
g
0
(
x
)
=
∏
i
=
1
n
(
x
−
a
i
)
∈
F
q
[
x
]
.
g_0(x)=\prod_{i=1}^{n}(x-a_i) \in \mathbb{F}_q[x].
g0(x)=∏i=1n(x−ai)∈Fq[x].
第一步:输入一个编码
(
c
1
,
c
2
,
⋯
,
c
n
)
(c_1,c_2,\cdots,c_n)
(c1,c2,⋯,cn),执行多项式插值,得到
g
1
(
x
)
∈
F
q
[
x
]
g_1(x) \in \mathbb{F}_q[x]
g1(x)∈Fq[x],使得:
c
i
=
g
1
(
a
i
)
,
1
≤
i
≤
n
.
c_i=g_1(a_i) , 1\leq i\leq n.
ci=g1(ai),1≤i≤n.
第二步:对
g
0
(
x
)
,
g
1
(
x
)
g_0(x),g_1(x)
g0(x),g1(x)使用扩展欧几里得除法,直到余数
g
(
x
)
g(x)
g(x)的次数小于
n
+
k
2
\frac{n+k}{2}
2n+k.这时有:
g
(
x
)
=
u
(
x
)
g
0
(
x
)
+
v
(
x
)
g
1
(
x
)
.
g(x)=u(x)g_0(x)+v(x)g_1(x).
g(x)=u(x)g0(x)+v(x)g1(x).
第三步:
g
(
x
)
g(x)
g(x)整除
v
(
x
)
v(x)
v(x),即
g
(
x
)
=
f
1
(
x
)
v
(
x
)
+
r
(
x
)
g(x)=f_1(x)v(x)+r(x)
g(x)=f1(x)v(x)+r(x)满足
v
(
x
)
v(x)
v(x)的次数大于
r
(
x
)
r(x)
r(x)的次数。
如果
f
1
(
x
)
f_1(x)
f1(x)的次数小于
k
k
k并且
r
(
x
)
=
=
0
r(x)==0
r(x)==0输出
f
1
(
x
)
f_1(x)
f1(x),其系数就是m,否则输出解码失败。
c++代码
该代码使用了NTL库
#include<random>
#include<NTL/ZZ_pEX.h>
#include<NTL/ZZ_pEXFactoring.h>
#include<NTL/ZZ_pXFactoring.h>
#include<NTL/ZZ_pX.h>
using namespace std;
using namespace NTL;
int n=15;
int k=7;
int prime=17;
Vec<ZZ_pX> partialXgcd(ZZ_pX a,ZZ_pX b,ZZ_pX u0=(ZZ_pX)1,ZZ_pX u1=(ZZ_pX)0,ZZ_pX v0=(ZZ_pX)0,ZZ_pX v1=(ZZ_pX)1){
ZZ_pX r;
ZZ_pX q;
DivRem(q,r,a,b);
if(deg(b)<(n+k)/2){
Vec<ZZ_pX> res(INIT_SIZE,3);
res[0]=b;
res[1]=u1;
res[2]=v1;
return res;
}
ZZ_pX u2=u0-q*u1;
ZZ_pX v2=v0-v1*q;
return partialXgcd(b,r,u1,u2,v1,v2);
}
ZZ_pX reedSolomon(vector<vector<ZZ_p>> points){
//ZZ_p::init(ZZ(prime));
ZZ_pX one(1);
ZZ_pX g0=one;
vec_ZZ_p as(INIT_SIZE,points.size()),bs(INIT_SIZE,points.size());
for(int i=0;i<points.size();i++){
ZZ_pX tt;
tt.SetLength(2);
SetCoeff(tt,1);
SetCoeff(tt,0,-points[i][0]);
g0=g0*tt;
as[i]=points[i][0];
bs[i]=points[i][1];
}
ZZ_pX g1=interpolate(as,bs);
auto res=partialXgcd(g0,g1);
res[0]=res[0]/res[2];
res[0].normalize();
if(deg(res[0])<k)
return res[0];
cout<<"Decoding error"<<endl;
return res[0];
}
int main(){
ZZ_p::init(ZZ(prime));
vector<vector<ZZ_p>> points(n,vector<ZZ_p>(2));
cout<<"Please input a degree "<<k<<" polynomial"<<endl; //the input should be like [1 2 3 4 5 6 7]
ZZ_pX poly;
cin>>poly;
for(int i=0;i<n;i++){
points[i][0]=i+1;
points[i][1]=eval(poly,(ZZ_p)i+1); //function eval computes the value of polynomial at a point.
}
ZZ_pX result=reedSolomon(points);
for(int i=0;i<n;i++){
ZZ_p tt(points[i][0]);
cout<<eval(result,tt)<<' '<<points[i][1]<<endl;
}
cout<<endl;
}