从多项式除法来看Reed

多项式除法来看Reed"/>

从多项式除法来看Reed

前言

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​+m2​x+⋯+mk​xk−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;
}