BCH码的基本原理参考了此处
【举例子详细分析】BCH码(BCH code)
模二多项式环
模二多项式环中的一个多项式，和二进制数 是有一个由多项式的系数决定的，一一对应的关系的。
比如，有一个多项式
将其系数为0的项补全后，可以得到
，
按顺序取出系数为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码
原理大概是，根据编码纠错个数， 有如下形式
要发送的信息，对其进行编码操作，得到的发送的信息S为
若发送过程不出错，则直接除以即可得到原文
但若发送信息的过程中收到干扰，某几位比特出现错误，即接受到的信息R为
由于之前和
的性质，
也满足
所以
这三个值都是已知的(可以求出的)，根据这三个值，通过某些方法解出 的值，即可找出错误位置并纠正，得到
，除以
即可得到原文
对信息编码
假设要发送的信息 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
原文中
此处 和
在
下是相同的
全部代码
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}')
文章出处登录后可见！
已经登录？立即刷新
