BCH码的基本原理参考了此处
【举例子详细分析】BCH码(BCH code)
模二多项式环
模二多项式环中的一个多项式,和二进制数 是有一个由多项式的系数决定的,一一对应的关系的。
比如,有一个多项式
F
(
x
)
=
x
4
+
x
+
1
F(x) = x^4 + x + 1
F(x)=x4+x+1
将其系数为0的项补全后,可以得到
F
(
x
)
=
1
∗
x
4
+
0
∗
x
3
+
0
∗
x
2
+
1
∗
x
+
1
F(x) = 1*x^4 + 0*x^3 + 0*x^2 + 1*x + 1
F(x)=1∗x4+0∗x3+0∗x2+1∗x+1,
按顺序取出系数为10011
模二多项式环中的运算有关特性不再赘述
代码实现
class Polynomial2():
'''
模二多项式环,定义方式有三种
一是
>>> Polynomial2([1,1,0,1])
x^3 + x^2 + 1
二是
>>> Polynomial2('1101')
x^3 + x^2 + 1
三是
>>> Poly([3,1,4]) # 直接给出系数为1的项的阶
x^4 + x^3 + x
>>> Poly([])
0
'''
def __init__(self,ll):
if type(ll) == str:
ll = list(map(int,ll))
self.param = ll[::-1]
self.ones = [i for i in range(len(self.param)) if self.param[i] == 1] # 系数为1的项的阶数列表
self.Latex = self.latex()
self.b = ''.join([str(i) for i in ll]) # 二进制形式打印系数
self.order = 0 # 最高阶
try:self.order = max(self.ones)
except:pass
def format(self,reverse = True):
'''
格式化打印字符串
reverse = False时,可以低位在左
但是注意定义多项式时只能高位在右
'''
r = ''
if len(self.ones) == 0:
return '0'
if reverse:
return ((' + '.join(f'x^{i}' for i in self.ones[::-1])+' ').replace('x^0','1').replace('x^1 ','x ')).strip()
return ((' + '.join(f'x^{i}' for i in self.ones)+' ').replace('x^0','1').replace('x^1 ','x ')).strip()
def __call__(self,x):# 懒得写
print(f'call({x})')
def __add__(self,other):
a,b = self.param[::-1],other.param[::-1]
if len(a) < len(b):a,b = b,a
for i in range(len(a)):
try:a[-1-i] = (b[-1-i] + a[-1-i]) % 2
except:break
return Polynomial2(a)
def __mul__(self,other):
a,b = self.param[::-1],other.param[::-1]
r = [0 for i in range(len(a) + len(b) - 1)]
for i in range(len(b)):
if b[-i-1] == 1:
if i != 0:sa = a+[0]*i
else:sa = a
sa = [0] * (len(r)-len(sa)) + sa
#r += np.array(sa)
#r %= 2
r = [(r[t] + sa[t])%2 for t in range(len(r))]
return Polynomial2(r)
def __sub__(self,oo):
# 在模二多项式环下,加减相同
return self + oo
def div(self,other):
r,b = self.param[::-1],other.param[::-1]
if len(r) < len(b):
return Polynomial2([0]),self
q=[0] * (len(r) - len(b) + 1)
for i in range(len(q)):
if len(r)>=len(b):
index = len(r) - len(b) + 1 # 确定所得商是商式的第index位
q[-index] = int(r[0] / b[0])
# 更新被除多项式
b_=b.copy()
b_.extend([0] * (len(r) - len(b)))
b_ = [t*q[i] for t in b_]
r = [(r[t] - b_[t])%2 for t in range(len(r))]
for j in range(len(r)): #除去列表最左端无意义的0
if r[0]==0:
r.remove(0)
else:
break
else:
break
return Polynomial2(q),Polynomial2(r)
def __floordiv__(self,other): # 只重载了整除,即//
return self.div(other)[0]
def __mod__(self,other):
return self.div(other)[1]
def __repr__(self) -> str:
return self.format()
def __str__(self) -> str:
return self.format()
def __pow__(self,a):
# 没有大数阶乘的需求,就没写快速幂
t = Polynomial2([1])
for i in range(a):
t *= self
return t
def latex(self,reverse=True):
# Latex格式打印...其实就是给大于一位长度的数字加个括号{}
def latex_pow(x):
if len(str(x)) <= 1:
return str(x)
return '{'+str(x)+'}'
r = ''
if len(self.ones) == 0:
return '0'
if reverse:
return (' + '.join(f'x^{latex_pow(i)}' for i in self.ones[::-1])+' ').replace('x^0','1').replace(' x^1 ',' x ').strip()
return (' + '.join(f'x^{latex_pow(i)}' for i in self.ones)+' ').replace('x^0','1').replace(' x^1 ',' x ').strip()
def __eq__(self,other):
return self.ones == other.ones
def Poly(ones):
if len(ones) == 0:
return Polynomial2('0')
ll = [0 for i in range(max(ones)+1)]
for i in ones:
ll[i] = 1
return Polynomial2(ll[::-1])
from tqdm import trange
from number import *
PP = Polynomial2
P = Poly
# 简化名称,按长度区分 P 和 PP
'''
定义方式可改为
PP('11010')
或者
P([5,6,7])
'''
另外,在sagemath中有更强大的更完备的多项式整数环
sage: P.<x> = PolynomialRing(Zmod(2))
sage: p = x^4 + x + 1
sage: p(x^3)
x^12 + x^3 + 1
sage: factor(p) # 分解多项式
x^4 + x + 1
sage: factor(x^4+x^2+1) # 分解多项式
(x^2 + x + 1)^2
对这个不太熟悉,但功能没问题,主要用来检验结果
BCH码
原理大概是,根据编码纠错个数,
Q
(
x
)
Q(x)
Q(x) 有如下形式
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-FT2gSQOi-1659454509052)(Untitled%202.assets/9c2418c5823d42d89f917cb668638fa4.png)]
p
(
x
)
p(x)
p(x) 是本原多项式primitive polynomial
选定
p
(
x
)
=
x
4
+
x
+
1
p(x)=x^4 + x + 1
p(x)=x4+x+1
p
3
(
x
)
p_3(x)
p3(x) 满足
p
3
(
x
3
)
≡
0
m
o
d
p
p_3(x^3)\equiv 0\mod p
p3(x3)≡0modp
选定
p
3
(
x
)
=
x
4
+
x
3
+
x
2
+
x
+
1
p_3(x)=x^4 + x^3 + x^2 + x + 1
p3(x)=x4+x3+x2+x+1
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-7Tg8CO4k-1659454509053)(Untitled%202.assets/9669eadfe0fc4983a64329cccd80c540.png)]
p
5
(
x
)
p_5(x)
p5(x) 满足
p
5
(
x
5
)
≡
0
m
o
d
p
p_5(x^5)\equiv 0\mod p
p5(x5)≡0modp
选定
p
5
(
x
)
=
x
2
+
x
+
1
p_5(x)=x^2 + x + 1
p5(x)=x2+x+1
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Nm5xHjrq-1659454509054)(Untitled%202.assets/d3bdf9ab65e84ffebd58f2e5a4d03f36.png)]
假设需要纠错3位,则
Q
(
x
)
=
p
(
x
)
∗
p
3
(
x
)
∗
p
5
(
x
)
Q(x)=p(x)*p_3(x)*p_5(x)
Q(x)=p(x)∗p3(x)∗p5(x)
要发送的信息
m
(
x
)
m(x)
m(x),对其进行编码操作,得到的发送的信息S为
S
(
x
)
=
Q
(
x
)
∗
m
(
x
)
S(x) = Q(x) * m(x)
S(x)=Q(x)∗m(x)
若发送过程不出错,则直接除以
Q
(
x
)
Q(x)
Q(x)即可得到原文
m
(
x
)
m(x)
m(x)
但若发送信息的过程中收到干扰,某几位比特出现错误,即接受到的信息R为
R
(
x
)
=
Q
(
x
)
∗
m
(
x
)
+
x
e
1
+
x
e
2
+
x
e
3
R(x) = Q(x) * m(x) + x^{e_1}+x^{e_2}+x^{e_3}
R(x)=Q(x)∗m(x)+xe1+xe2+xe3
由于之前
p
3
(
x
)
p_3(x)
p3(x)和
p
5
(
x
)
p_5(x)
p5(x)的性质,
Q
(
x
)
Q(x)
Q(x)也满足
Q
(
x
3
)
≡
0
m
o
d
p
Q
(
x
5
)
≡
0
m
o
d
p
Q(x^3)\equiv0\mod p\\ Q(x^5)\equiv0\mod p
Q(x3)≡0modpQ(x5)≡0modp
所以
R
(
x
)
≡
x
e
1
+
x
e
2
+
x
e
3
m
o
d
p
R
(
x
3
)
≡
x
3
⋅
e
1
+
x
3
⋅
e
2
+
x
3
⋅
e
3
m
o
d
p
R
(
x
5
)
≡
x
5
⋅
e
1
+
x
5
⋅
e
2
+
x
5
⋅
e
3
m
o
d
p
R(x)\equiv x^{e_1}+x^{e_2}+x^{e_3}\mod p\\ R(x^3)\equiv x^{3·e_1}+x^{3·e_2}+x^{3·e_3}\mod p\\ R(x^5)\equiv x^{5·e_1}+x^{5·e_2}+x^{5·e_3}\mod p
R(x)≡xe1+xe2+xe3modpR(x3)≡x3⋅e1+x3⋅e2+x3⋅e3modpR(x5)≡x5⋅e1+x5⋅e2+x5⋅e3modp
这三个值都是已知的(可以求出的),根据这三个值,通过某些方法解出
e
1
,
e
2
,
e
3
e_1,e_2,e_3
e1,e2,e3 的值,即可找出错误位置并纠正,得到
S
(
x
)
S(x)
S(x),除以
Q
(
x
)
Q(x)
Q(x)即可得到原文
m
(
x
)
m(x)
m(x)
对信息编码
假设要发送的信息 m = 11010
这里原paper里面是低位在左边
这里需要反序生成多项式
p = Polynomial2('10011')
p3 = Polynomial2('11111')
p5 = Polynomial2('111')
Q = p * p3 * p5
print('p =',p.b)
print('p =',p)
print('p3 =',p3.b)
print('p5 =',p5.b)
m = PP('01011')
send = Q*m
print('send =',send.b)
p = 10011
p = x^4 + x + 1
p3 = 11111
p5 = 111
send = 010011011100001
ei1 = 2
ei2 = 5
e1 = Poly([ei1])
e2 = Poly([ei2])
r = send + e1 + e2
# r.ones是r中系数为1的项的阶数
# [3*i for i in r.ones]就是把所有阶数乘了3
r3 = Poly([3*i for i in r.ones]) # R(x^3)
r5 = Poly([5*i for i in r.ones]) # R(x^5)
rx1 = r%p
rx3 = r3%p
rx5 = r5%p
print('rx1 =',rx1)
print('rx3 =',rx3)
print('rx5 =',rx5)
rx1 = x
rx3 = x^3 + x^2 + 1
rx5 = 0
原文中
R
(
x
3
)
m
o
d
p
=
x
13
R(x^3) \mod p=x^{13}
R(x3)modp=x13
此处
x
3
+
x
2
+
1
x^3 + x^2 + 1
x3+x2+1 和
x
13
x^{13}
x13 在
m
o
d
p
\mod p
modp下是相同的
全部代码
from random import randint
DEBUG = 0
# 如果开启DEBUG模式,则会打印中间值信息,并且按enter才会进入下一次恢复
# 不开启DEBUG,则while 1不断运行,捕捉到Ctrl+C中断时打印当前统计的成功率
class Polynomial2():
'''
模二多项式环,定义方式有三种
一是
>>> Polynomial2([1,1,0,1])
x^3 + x^2 + 1
二是
>>> Polynomial2('1101')
x^3 + x^2 + 1
三是
>>> Poly([3,1,4]) # 直接给出系数为1的项的阶
x^4 + x^3 + x
>>> Poly([])
0
'''
def __init__(self,ll):
if type(ll) == str:
ll = list(map(int,ll))
self.param = ll[::-1]
self.ones = [i for i in range(len(self.param)) if self.param[i] == 1] # 系数为1的项的阶数列表
self.Latex = self.latex()
self.b = ''.join([str(i) for i in ll]) # 二进制形式打印系数
self.order = 0 # 最高阶
try:self.order = max(self.ones)
except:pass
def format(self,reverse = True):
'''
格式化打印字符串
reverse = False时,可以低位在左
但是注意定义多项式时只能高位在右
'''
r = ''
if len(self.ones) == 0:
return '0'
if reverse:
return ((' + '.join(f'x^{i}' for i in self.ones[::-1])+' ').replace('x^0','1').replace('x^1 ','x ')).strip()
return ((' + '.join(f'x^{i}' for i in self.ones)+' ').replace('x^0','1').replace('x^1 ','x ')).strip()
def __call__(self,x):# 懒得写
print(f'call({x})')
def __add__(self,other):
a,b = self.param[::-1],other.param[::-1]
if len(a) < len(b):a,b = b,a
for i in range(len(a)):
try:a[-1-i] = (b[-1-i] + a[-1-i]) % 2
except:break
return Polynomial2(a)
def __mul__(self,other):
a,b = self.param[::-1],other.param[::-1]
r = [0 for i in range(len(a) + len(b) - 1)]
for i in range(len(b)):
if b[-i-1] == 1:
if i != 0:sa = a+[0]*i
else:sa = a
sa = [0] * (len(r)-len(sa)) + sa
#r += np.array(sa)
#r %= 2
r = [(r[t] + sa[t])%2 for t in range(len(r))]
return Polynomial2(r)
def __sub__(self,oo):
# 在模二多项式环下,加减相同
return self + oo
def div(self,other):
r,b = self.param[::-1],other.param[::-1]
if len(r) < len(b):
return Polynomial2([0]),self
q=[0] * (len(r) - len(b) + 1)
for i in range(len(q)):
if len(r)>=len(b):
index = len(r) - len(b) + 1 # 确定所得商是商式的第index位
q[-index] = int(r[0] / b[0])
# 更新被除多项式
b_=b.copy()
b_.extend([0] * (len(r) - len(b)))
b_ = [t*q[i] for t in b_]
r = [(r[t] - b_[t])%2 for t in range(len(r))]
for j in range(len(r)): #除去列表最左端无意义的0
if r[0]==0:
r.remove(0)
else:
break
else:
break
return Polynomial2(q),Polynomial2(r)
def __floordiv__(self,other): # 只重载了整除,即//
return self.div(other)[0]
def __mod__(self,other):
return self.div(other)[1]
def __repr__(self) -> str:
return self.format()
def __str__(self) -> str:
return self.format()
def __pow__(self,a):
# 没有大数阶乘的需求,就没写快速幂
t = Polynomial2([1])
for i in range(a):
t *= self
return t
def latex(self,reverse=True):
# Latex格式打印...其实就是给大于一位长度的数字加个括号{}
def latex_pow(x):
if len(str(x)) <= 1:
return str(x)
return '{'+str(x)+'}'
r = ''
if len(self.ones) == 0:
return '0'
if reverse:
return (' + '.join(f'x^{latex_pow(i)}' for i in self.ones[::-1])+' ').replace('x^0','1').replace(' x^1 ',' x ').strip()
return (' + '.join(f'x^{latex_pow(i)}' for i in self.ones)+' ').replace('x^0','1').replace(' x^1 ',' x ').strip()
def __eq__(self,other):
return self.ones == other.ones
def Poly(ones):
if len(ones) == 0:
return Polynomial2('0')
ll = [0 for i in range(max(ones)+1)]
for i in ones:
ll[i] = 1
return Polynomial2(ll[::-1])
from tqdm import trange
from number import *
PP = Polynomial2
P = Poly
# 简化名称,按长度区分 P 和 PP
p = Polynomial2('10011')
p3 = Polynomial2('11111')
Q = p*p3
def getMorder(r,e1e2):
return ((r+Poly(e1e2))//Q).order
def main(mm):
m = Polynomial2(bin(mm)[2:])
send = Q*m # 发送的Send多项式
ie1 = randint(0,send.order) # 随机选两个不相同的出错位置
ie2 = randint(0,send.order)
while ie1 == ie2:
ie2 = randint(0,send.order)
e1 = Poly([ie1])
e2 = Poly([ie2])
recv = send + e1 + e2 # 接收到的出错多项式
rx = Poly([i for i in recv.ones])
rx3 = Poly([3*i for i in recv.ones])
if DEBUG:
ff = max(len(recv.b),len(send.b))
print('mm =',bin(mm)[2:])
print('m =',m)
print('p =',p)
print('p3 =',p3)
print()
print('send:\t',send.b.zfill(ff))
print('recv:\t',recv.b.zfill(ff))
print('e:\t',(e1+e2).b.zfill(ff))
print('e1,e2 =',ie1,ie2)
print()
print('R(x) =',rx)
print('R(x) =',rx%p)
print('R(x^3) =',rx3)
print('R(x^3) =',rx3%p)
ans = []
for i in range(17):
_e1 = Poly([i])
for j in range(i,17):
if i == j:continue
_e2 = Poly([j])
t1 = (_e1+_e2)%p
t3 = (Poly([3*i])+Poly([3*j]))%p
if DEBUG and i == ie1 and j == ie2:
print('t1: ',t1)
print('t3: ',t3)
print(t1 == rx%p,t3 == rx3%p)
if t1 == rx%p and t3 == rx3%p:
#print('find',int(_e1.b,2),int(_e2.b,2),i,j)
ans.append((i,j))
if DEBUG:print('过滤前ans:',ans)
if len(ans) > 1:
# 若找到多组 e1e2,则判断(R-e1-e2)//p 的阶数是否小于8
# 因为 m 的阶数不会超过8
res = [i for i in ans if getMorder(rx,i)<=7]
else:
res = ans
if DEBUG:
print('过滤后ans:',res,'出错位置:',(ie1,ie2))
input()
return res
if __name__ == '__main__':
cnt = [0,0,0,0,0]
while 1:
try:
# 随机生成一个8bit明文
cnt[len(main(randint(0,255)))] += 1
# 返回值是结果列表 如果列表里只有一个元素则恢复成功
# cnt 记录的是返回的列表的长度
except:
print('\n',cnt)
if sum(cnt) != 0:
sucP = (cnt[1]+cnt[2]/2+cnt[4]/4)/sum(cnt)
print(f'p : {sucP},p^25:{sucP**25}')