模二多项式环

模二多项式环中的一个多项式，和二进制数是有一个由多项式的系数决定的，一一对应的关系的。
比如，有一个多项式
F(x)=x4+x+1F(x) = x^4 + x + 1F(x)=x4+x+1
将其系数为0的项补全后，可以得到
F(x)=1∗x4+0∗x3+0∗x2+1∗x+1F(x) = 1*x^4 + 0*x^3 + 0*x^2 + 1*x + 1F(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:passdef 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,afor i in range(len(a)):try:a[-1-i] = (b[-1-i] + a[-1-i]) % 2except:breakreturn 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]*ielse:sa = asa = [0] * (len(r)-len(sa)) + sa#r += np.array(sa)#r %= 2r = [(r[t] + sa[t])%2 for t in range(len(r))]return Polynomial2(r)def __sub__(self,oo):# 在模二多项式环下，加减相同return self + oodef div(self,other):r,b = self.param[::-1],other.param[::-1]if len(r) < len(b):return Polynomial2([0]),selfq=[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)):     #除去列表最左端无意义的0if r[0]==0:r.remove(0)else:breakelse:breakreturn 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 *= selfreturn tdef 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.onesdef Poly(ones):if len(ones) == 0:return Polynomial2('0')ll = [0 for i in range(max(ones)+1)]for i in ones:ll[i] = 1return 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 + 1sage: p(x^3)
x^12 + x^3 + 1sage: factor(p) # 分解多项式
x^4 + x + 1sage: 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)=x4+x+1p(x)=x^4 + x + 1p(x)=x4+x+1
p3(x)p_3(x)p3(x) 满足 p3(x3)≡0modpp_3(x^3)\equiv 0\mod pp3(x3)≡0modp
选定 p3(x)=x4+x3+x2+x+1p_3(x)=x^4 + x^3 + x^2 + x + 1p3(x)=x4+x3+x2+x+1
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-7Tg8CO4k-1659454509053)(Untitled%202.assets/9669eadfe0fc4983a64329cccd80c540.png)]
p5(x)p_5(x)p5(x) 满足 p5(x5)≡0modpp_5(x^5)\equiv 0\mod pp5(x5)≡0modp
选定 p5(x)=x2+x+1p_5(x)=x^2 + x + 1p5(x)=x2+x+1
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Nm5xHjrq-1659454509054)(Untitled%202.assets/d3bdf9ab65e84ffebd58f2e5a4d03f36.png)]
假设需要纠错3位，则
Q(x)=p(x)∗p3(x)∗p5(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)+xe1+xe2+xe3R(x) = Q(x) * m(x) + x^{e_1}+x^{e_2}+x^{e_3} R(x)=Q(x)∗m(x)+xe1+xe2+xe3

由于之前p3(x)p_3(x)p3(x)和p5(x)p_5(x)p5(x)的性质，Q(x)Q(x)Q(x)也满足
Q(x3)≡0modpQ(x5)≡0modpQ(x^3)\equiv0\mod p\\ Q(x^5)\equiv0\mod p Q(x3)≡0modpQ(x5)≡0modp
所以
R(x)≡xe1+xe2+xe3modpR(x3)≡x3⋅e1+x3⋅e2+x3⋅e3modpR(x5)≡x5⋅e1+x5⋅e2+x5⋅e3modpR(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
这三个值都是已知的(可以求出的)，根据这三个值，通过某些方法解出 e1,e2,e3e_1,e_2,e_3e1,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 * p5print('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%pprint('rx1 =',rx1)
print('rx3 =',rx3)
print('rx5 =',rx5)

rx1 = x
rx3 = x^3 + x^2 + 1
rx5 = 0

原文中 R(x3)modp=x13R(x^3) \mod p=x^{13}R(x3)modp=x13
此处 x3+x2+1x^3 + x^2 + 1x3+x2+1 和 x13x^{13}x13 在modp\mod pmodp下是相同的

全部代码

from random import randintDEBUG = 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:passdef 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,afor i in range(len(a)):try:a[-1-i] = (b[-1-i] + a[-1-i]) % 2except:breakreturn 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]*ielse:sa = asa = [0] * (len(r)-len(sa)) + sa#r += np.array(sa)#r %= 2r = [(r[t] + sa[t])%2 for t in range(len(r))]return Polynomial2(r)def __sub__(self,oo):# 在模二多项式环下，加减相同return self + oodef div(self,other):r,b = self.param[::-1],other.param[::-1]if len(r) < len(b):return Polynomial2([0]),selfq=[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)):     #除去列表最左端无意义的0if r[0]==0:r.remove(0)else:breakelse:breakreturn 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 *= selfreturn tdef 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.onesdef Poly(ones):if len(ones) == 0:return Polynomial2('0')ll = [0 for i in range(max(ones)+1)]for i in ones:ll[i] = 1return Polynomial2(ll[::-1])from tqdm import trange
from number import *PP = Polynomial2
P = Poly
# 简化名称,按长度区分 P 和 PPp = Polynomial2('10011')
p3 = Polynomial2('11111')
Q = p*p3def getMorder(r,e1e2):return ((r+Poly(e1e2))//Q).orderdef 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)%pt3 = (Poly([3*i])+Poly([3*j]))%pif 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 = ansif DEBUG:print('过滤后ans:',res,'出错位置:',(ie1,ie2))input()return resif __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}')