快速相關(guān)攻擊算法A和B的實驗報告

第二屆全國密碼數(shù)學挑戰(zhàn)賽西南賽區(qū)

題目二

參賽單位:電子科技大學

參賽隊員:詹彤、王雨飛、方淳晟

指導(dǎo)老師：________________________

聯(lián)系電話：XXXXXXXXXXXXX

聯(lián)系由B箱：XXXXXXXXXXXXXXXXXXXXXX

摘要

快速相關(guān)攻擊算法是解決流密碼問題的一種重要方法，也是密碼學中重要研

究內(nèi)容之一。自Siegenthaler提出快速相關(guān)攻擊以來，不斷有更高效的算法被

提出。

本文基于Meier和Staffelbach提出的Fastcorrelationattacks中的兩

種算法，根據(jù)提供的被加密的序列密碼以及其和原序列的相關(guān)概率，利用概率等

相關(guān)知識，恢復(fù)線性移位寄存器的60位初始序列。

本文所用的算法，不僅可以解決快速相關(guān)攻擊中，密碼和輸出序列的相關(guān)概

率較高的流密碼問題，也可以解決對于線性移位寄存器的反饋多項式而言，抽頭

數(shù)不高的流密碼。因此有著廣泛的應(yīng)用前景。

2

目錄

1.序列密碼的相關(guān)攻擊介紹..........................4

2.快速相關(guān)攻擊研究現(xiàn)狀............................5

3.Meier-Staffelbach型算法的基本思想................5

4.快速相關(guān)攻擊算法I................................................................6

4.1算法I實現(xiàn)步驟..........................6

4.2算法I計算復(fù)雜度分析....................7

5.快速相關(guān)攻擊算法H................................................................8

5.1算法H實現(xiàn)步驟...........................8

5.2算法II計算復(fù)雜度分析....................9

6.結(jié)果部分.......................................10

7.相關(guān)代碼.......................................11

7.1校驗方程生成代碼.........................12

7.2算法I實現(xiàn)代碼...........................14

7.3算法I實現(xiàn)代碼...........................16

8.參考文獻.......................................20

3

1.序列密碼的相關(guān)攻擊介紹

本問題是基于線性反饋移位寄存器LFSR的相關(guān)攻擊問題。令LFSR的反饋多

項式為/*)，一個60位的初始密鑰a,a，…,a,經(jīng)過LFSR的線性反饋多項式

0I59

的運算，輸出序列u=(u,u,...u)，經(jīng)由二元無記憶對稱信道(BSC)加密，得到截

0IN-I

取序列z=(z,z,...z)o實現(xiàn)過程如圖1所示。

0IN-I

初始狀態(tài)為

(a。,Q1，...,。59)

圖1加密過程演示

如果線性函數(shù)/的某個輸入分量U與輸出Z之間存在相關(guān)性，即

nn

p=[P(u=z)]>0.5。那么可以窮搜索LFSR的初始序列，找出使LFSR的輸出與對

nn

應(yīng)密鑰流序列Z的距離最小的初始序列作為該LFSR的初始序列(種子密鑰)。

假設(shè)LFSR產(chǎn)生最大周期序列，設(shè)級數(shù)為n,則周期為2「1。該系統(tǒng)的密鑰量為:

K=2..-l(1)

如果窮搜索密鑰攻擊算法，n足夠大的時候，計算量是無法實現(xiàn)的，事實上,

隨著n的增大，算法的復(fù)雜度呈指數(shù)增長。如果相關(guān)攻擊不對種子密鑰進行枚舉

搜索，就稱為快速相關(guān)攻擊。

2.快速相關(guān)攻擊研究現(xiàn)狀

相關(guān)攻擊最早是由Siegenthaler提出，其原理是利用密鑰發(fā)生器的輸出序

列與其源LFSR的輸出序列之間具有的相關(guān)性還原LFSR的初始狀態(tài)。Meier和

Staffelbach在此基礎(chǔ)上提出利用概率迭代譯碼的快速相關(guān)攻擊算法⑴，該方法

4

接收到的密鑰流序列的每一比特位代入原序列的反饋多項式組成的方程組中，利

用方程組是否滿足原序列的反饋多項式的概率，通過一定次數(shù)的迭代，來最終找

到序列的初始狀態(tài)，實現(xiàn)解密。在LFSR的特征多項式重量較小(<10)的情況

下取得了較好的效果，之后有學者提出了一些改進算法。

LFSR的抽頭數(shù)一直是Meier-Staffelbach型相關(guān)攻擊算法的頸瓶。1999

年T.Johansson和F.Jonsson把快速相關(guān)攻擊應(yīng)用到LFSR的一般情形，他們的

攻擊的關(guān)鍵是把LFSR序列轉(zhuǎn)化成卷積碼，使用卷積碼的譯碼算法恢復(fù)LFSR的初

態(tài)，仿真結(jié)果表明，使用Viterbi譯碼算法，在抽頭數(shù)較大時，和以往算法比較，

攻擊成功時的相關(guān)系數(shù)減小；接著他們把卷積碼和迭代概率譯碼融合起來，把

LFSR序列轉(zhuǎn)變成turbo碼，繼而使用turbo譯碼算法恢復(fù)LFSR的初態(tài)。2000

年，V.Chepyzhov,T.Joansson和B.Smeets的部分思想⑵提出一個較為簡單

而有效的快速相關(guān)攻擊算法。T.Johansson和F.Jonsson提出不同于以往的BSC

攻擊模型即線性多項式重構(gòu)模型，吸取了多項式學習的已有成果進行攻擊。

近年來，快速相關(guān)攻擊技術(shù)有了很大的發(fā)展，許多文獻給出了新的快速相關(guān)

攻擊算法。大多數(shù)快速相關(guān)攻擊算法利用窮舉搜索猜測LFSR部分初始狀態(tài)，再

結(jié)合密鑰流和校驗方程集合逐步恢復(fù)其余的LFSR初始狀態(tài)?？焖傧嚓P(guān)攻擊新的

發(fā)展方向是：算法可由高速軟件或簡單硬件實現(xiàn)，并適于并行運算。

3.基于Meier-Staffelbach型的算法

本文提到的算法I和算法II,其基本思想是考慮到輸出序列z和原序列a之

間存在相關(guān)性〃，通常情況下p>0.53,因此可以將輸出序列z看作是原序列a的

一種擾亂，故可以把序列z的每一比特位代入LFSR反饋多項式組成的方程組中。

考慮在模二域GF(2)上有

f2i(X)=f(X2i)(2)

所以通過對/(x)的重復(fù)乘方可以得到多個LFSR的校驗多項式，對每一個街

區(qū)序列z上的比特位而言，如果某一比特位上，關(guān)于該比特位，這些校驗多項式

5

成立的個數(shù)越多，用該比特位代替LFSR序列的相應(yīng)比特是正確的概率就越大,

在求解概率問題中主要利用形如式⑶的概率迭代公式

S(p/)=pS(p,f-l)+(l-p)(l—S(p,f—1))(3)

以及求解條件概率公式

P(AI3)=上空2(4)

P(B)

貝葉斯公式主要用于求解算法n中的概率問題

-=咽)一⑷外⑸

iZ"P(B)P(AIB)

j=lJj

而上面提到的漢明距離，既是指表示兩個(相同長度)字對應(yīng)位不同的數(shù)量。

4.快速相關(guān)攻擊I

4.1算法I實現(xiàn)步驟

設(shè)反饋多項式級數(shù)為〃，抽頭數(shù)為f,相關(guān)概率p=P(〃=z),截取序列z的

nn

長度為N。

第1步：計算序列z每比特所需要的校驗方程平均數(shù)：

N

M=log(―)x(r+l)(6)

22n

第2步：計算截取序列名的比特的模二加(含有r個比特)和原LFSR序列的

相應(yīng)比特的模二加相同的概率S：

S(p,l)=p

s(pj)=pS(p,f-1)+(1-p)(l-S(p,f-l))(7)

第3步：計算在M個校驗方程中，至少有6個方程成立的概率Q(p,M,/O：

Q(p,",〃)=?。(q?S，x(l—q)?(l—S)，XSMT)

(8)

6

第4步：計算在M個校驗方程中，至少有〃個方程成立，且它是正確比特位

的概率，即z="概率為

nn

V(p,M,/?)=十C-L(q?S，x(l—S)MT

M

i=h(9)

因此，在給定M個方程中至少有八個方程成立的條件下，z=〃的概率為:

nn

T(P，M,/?)=E(10)

Q

第5步：求出一個最大的力，記為人,使。*N2〃。計算在〃個比特中錯誤

max

的平均數(shù)r：

r=(1-T)xn(11)

第6步：找出在M個方程中至少滿足〃個方程的比特位，用這些比特位作為

max

相應(yīng)下標位置上LFSR序列的“正確”比特，然后適當?shù)剡x取包含〃個比特的一

個集合記為I。，固定〃個連續(xù)比特位。

第7步：對I。中的每一比特位利用LFSR的反饋多項式組建成一個線性方程，

即把它表示成固定〃個比特位的一個線性組合，然后將這〃個線性方程組合在一

起組成一個線性方程組。

第8步：解線性方程，若是線性相關(guān)的，則用幾個附加比特位選擇一個線性無

關(guān)的子方程組，將解出來的"個連續(xù)比特擴展為整個LFSR序列，并與原來截取

到的序列z的N位進行比較，若相關(guān)概率在(p-e,p+e)之間(通常認為e取

0.05),則認為求解成功。

若相關(guān)概率不在區(qū)間之內(nèi)，則認為求解失敗，從而以漢明距1,2,...檢

驗/的所有變換，每進行以此漢明變換都轉(zhuǎn)入第7步。直到成功求解為止。

()

7

4.2算法I計算復(fù)雜度分析

算法I的主要求解時間耗費在進行漢明變換上。我們可以假設(shè)解方程和添加

附加位這一過程所用時間為，進行漢明變換的距離為"，則總時間為

T=A(n,d)xt=Cixt

i=0"(12)

則有TW2H<C,其中H(e)表示二進制的端函數(shù)，e=L整個算法復(fù)雜度為

n

0(2d),C取決于p、t、n,N,且0VCV1。

5.快速相關(guān)攻擊H

該算法的基本思想：根據(jù)截取序列的每一比特位所滿足的方程數(shù)。來計算截取

序列每一比特的新概率，當這些新概率小于給定值的數(shù)目超過另一給定值，或

者計算新概率的總的次數(shù)超過5,就對新概率小于給定值的截取序列的比特位進

行取補，接著以新序列代替原序列，并將各比特位的概率重新置為原相關(guān)概率，

重復(fù)上述過程，直到成功。

5.1算法n實現(xiàn)步驟：

第1步：計算序列z每比特所需要的校驗方程平均數(shù)

M=log2(N/2*n))*(t+l)(13)

第2步：計算截取序列z的比特的模二加(含有，個比特)和原LFSR序列的

相應(yīng)比特的模二加相同的概率S

s(p,D=p

S(p,r)=pxS(p,I)+(1-p)G-S(p1-D)(14)

第3步：計算在M個校驗方程中，至少有九個方程成立的概率。

8

Q(p,M,h)=Z。(qxS，x(l-S)MT+(1-q)x(l-S)ixS*)

M

i=h(15)

計算在"個校驗方程中，對不滿足〃個方程的比特位進行取補，取補后增

加的正確的比特位概率為

I(p,M,h)=ZC'(l-q)x(l-S)iC>?xS<X(1-S)M-,

MM

i=h\=h(16)

計算在〃個校驗方程中，有九個方程成立時每個比特位的后驗概率為

pn(p,M,h)=pSh(l-S)m-h/(pSh(1-S)m-h+(1-/?)(1-S)h(1-S)m-h)

(17)

第4步：找出使得I最大的人記為力，若I(p,M,h)<0,則表示該算法失

maxmax

敗，若>0,則計算P,用以表示各比特位的一個閾值，若后驗概率P<P,

thrnthr

就對該比特位進行取補

P=(p(p,A/,〃+l)+p{p,M,h))/2(18)

thrnmaxnmax

第5步：用N表示P<P的比特位的數(shù)目N,N>N則可結(jié)束此次

ihrnthrwwthr

迭代，也就是直接對尸<p的比特位進行取補，開始下一次迭代：

nihr

N=Q(p,M,h)xN(19)

thrmax

第6步：重置迭代次數(shù)i為。

第7步：對截取序列z的每一位比特位和原LFSR序列相應(yīng)比特位的和是相

同的概率S(p,p，…，p,t):

12t

S(pA)=p,

...(20)

S(p,p,...p,t)=pS(p,p,/-I)+(1-/?)(\-S(p,p

12tt]2Mt12f-1

第8步：由貝葉斯公式計算各比特位的新概率p*：

9

Y=pxSx...Sx(l-S4-l)x...x(l-S)

ijijhjhjM(21)

pi*=Y/(Y+(l-p)x(l-S)x...(l-S)x(l-s)x...xS))

ijijhj(h+l)jM

pi*：表示第i位比特的薪概率，而pi表示元概率

jl,j2…jh：在統(tǒng)計校驗方程成立的數(shù)目時，其表示使得方程成立的S

的下標。而剩下的人則表示不成立的S的下標。

第9步：確定P<P的總比數(shù)目N，若N>N或,Ya,則i++,轉(zhuǎn)到第7

nfhrwwthr

步。否則，對截取序列z的那些滿足p<p的比特位進行取補，并重置比特位

nthr

的概率位P。

第10步：若取補后的截取序列中出現(xiàn)了有不滿足反饋多項式的比特位時，

則轉(zhuǎn)到第6步。直到攻擊成功。

5.2算法n計算復(fù)雜度分析

其算法復(fù)雜度為oOvMf）。在p、f給定的情況下，復(fù)雜度隨LFSR的級數(shù)

增大而增大。在其他條件不變的情況下，復(fù)雜度隨，增大而增大，隨p增大而減

小。缺點也很明顯，在迭代過程中，可能出現(xiàn)分母為0的情況。

6.結(jié)果部分

下面給出各序列的從。至h的解。

059

（1）第一個序列Z

0,0,1,0,1,0,1,1,1,1,1,0,0,1,1,0,0,0,1,1,1,1,1,0,0,1,0,0,0,0,1,0,1,1,0,1,1,1,0,001,1,1,1,0

,1,0,1,1,0,0,0,1,0,0,0,1,1,1

（2）第二個序列Z

1,1,0,0,0,1,1,1,0,0,0,0,0,1,0,0,1,0,0,1,1,0,1,0,1,1,1,0,1,1,0,0,1,1,0,1,0,0,0,1,1,1,1,0,0,0

D

,0,1,0,0,1,1,0,1,1,0,0,0,1,1

（3）第三個序列z

1,1,0,1,1,1,1,1,0,0,0,1,1,1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,14,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1

,0,1,1,0,0,0,1,1,0,0,1,1,0,1

（4）第四個序列Z

1,0,1,1,1,0,0,1,1,0,1,0,1,1,0,1,1,1,0,1,0,1,1,0,0,100,1,1,0,1,1,1,0,0,1,1,1,1,1,0,1,1,1,0

,1,1,0,14,1,0,1,1,1,0,1,0,0,0

（5）第五個序列Z

1,0,1,0,1,0,1,1,0,1,1,0,0,1,1,0,0,0,1,1,1,1,1,0,0,1,0,0,0,0,1,0,1,1,0,1,1,1,0,0,0,1,1,1,1,0

（6）第六個序列Z

1,1,0,0,0,1,1,1,1,0,0,0,0,1,0,0,1,0,0,1,1,0,1,0,1,1,001,1,0,0,1,1,0,100,0,1,1,1,1,0,0,0

,0,1,0,0,1,1,0,140,0,0,1,1

（7）第七個序列Z

1,0,0,1,1,1,1,1,0,0,0,1,1,1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,1,14,0,0,1,0,0,0,1,1,0,0,0,1,1,1

,1,1,1,0,0,04,1,0,0,1,1,0,1

（8）第八個序列Z

0,o,1,1,0,1,LLL0,1,1,1,0,1,0,1,1,0,0,1,0,0,l,l,0,LLl,0,0,1,1,1,1,1,0,1,1,1,0,1

,1,0,1,1,1,0,141,0,1,0,0,0

（9）第九個序列Z

1,0,0,1,1,0,0,0,1,卬,1,0,0,1,0,0,0,0,1,0,1,1,0,1,11,0,0,0,1,1,1,1,0

30,1,1,1,0,0,1,0,0,0,1,1,1

（10）第十個序列z

1,1,0,0,0,1,1,1,0,0,0,1,0,1,001,0,0,1,1,0,1,0,1,1,1,0,1,1,0,0,1,1,0,1,0,0,0,1,1,1,1,0,0,0

,o,iAO,i,i,ojzi,o,o,o,i,i

（II）第H-一個序列Z

1,1,0,1,1,1,1,1,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,04,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1

,0,1,1,0,0,0,1,1,0,0,1,1,0,1

（12）第十二個序列Z

O,1,1,WU,1,0,1,0,1,1,0,1,1,1,0,1,0,1,1,0,0,1,1,0,1,1,0,1,1,1,0,0,1,1,1,0,1

n

,IQ,LI,141,1,1,0,1,0,0,0

7.相關(guān)代碼

7.1生成校驗方程

importscipy

importsys,time

fromscipy.specialimportcomb

importre

importmath

N=4000000#序列數(shù)

n=60#級數(shù)

p=0.75#相關(guān)概率

fx=[3,9,16,35,37,47,56,60]力反饋多項式的系數(shù)

#fx=[l,2,3,5]

t=len(fx)忻I?算抽頭數(shù)

importscipy

importsys,time

importre

importmath

#讀取序列z

defReadFile(filename):

datafile=open(filename)

try:

text=datafile.read()

returnevalCf+text[4：-1]+*Y)

finally:

datafile.close()

defwriteFile(filename,datawrite):

try:

file_object=open(filename,'w')

#file_object.write(10.75\n，+str(data).replace','')，')

file_object.writeC0.75\n+str(data_write).replace(***)[1：-l]+

finally:

file_object.closeO

N=4000000與序列數(shù)

艾

n=60#級數(shù)

fx=[3,9,16,35,37,47,56,60]#反饋多項式的系數(shù)

#fx=[l,2,3,5]

t=len(fx)#計算抽頭數(shù)

M=round(math,log(N/(2*n),2)*(t+1))#計算每位平均需要校驗方程數(shù)

imax=math.floor(math.log(N/n,2))

an=[[0foriinrange(t+1)]foriinrange(3*M)]#用「存儲最終校驗多項式

#print(imax)

defCreatePoly(w):#w為檢測第幾位

imax=math.floor(math.log(N/n,2))#乘方的最大值

count=0#記錄校驗方程的個數(shù)

al=[[0foriinrange(t)]forinrange(imax+1)]#用以保存校驗方程的序號

al[0]=fx

#生成倍式(第?行為原反饋多項式)

foriinranged,imax+1):

forjinrange(t):

al[i][j]=al[0][j]*pow(2,i)

#print(al)

性成校驗方程

a2=[[0foriinrange(t+1)]for:inrange(imax+1)]

#print(a2)

forcolumninrange(imax+1):#將校驗等式左邊寫通2的最后一列

a2[column][t]=al[column][t-1]

#print(a2)

foriinrange(imax+1):

forjinrange(t):

a2[i][t-j-l]=a2[i][t]-al[i][j]

#print(a2)

foriinrange(imax+1):

forjinrange(t+1):

ifw>=a2[i][j]:

offset=w-a2[i][j]

if(offset+a2[i][t])>N:

continue

else:

forkinrange(t+1):

an[count][k]=offset+a2[i][k]

count=count+l

returncount#返回校驗方程總數(shù)

B

deffindBitel():

rightnum=[]

data=ReadFileCdataT.txt")

forwinrange(4000000):

count=0

a=CreatePoly(w)

ifa>=92:

foriinrange(a):

try:

if

data[an[i][0]]^data[an[i][1]]data[an[i][2]]"data[an[i][3]]data[an[i][4]]\

data[an[i][5]]"data[an[i][6]]"data[an[i][7]]==data[an[i][8]]:

#f=open(zztestl,txt”,'a)

#f.write(str(an[i])+>,')

count=count+l

#f.close()

except:

print(an

print(count)

ifcount>92:

rightnum.append(w)

print(w)

deffindBite2():

f=open(*test2.txt*)

a=f.read()

data=eval(a)

f.close()

b=len(data)

an=[[0foriinrange(60)]]

foriinrange(60):

count=0

forkinrange(10000):

forjinrange(9):

ifdata[k][j]==12800+i:

anti]=12800+i

count=count+1

print("找到"，count,"個了")

break

■

print(an)

7.2算法I實現(xiàn)代碼

importscipy

importsys,time

fromscipy.specialimportcomb

importre

importmath

N=4000000#序列數(shù)

n=60#級數(shù)

P=0.75#相關(guān)概率

fx=[3,9,16,35,37,47,56,60]#反饋多項式的系數(shù)

#fx=[l,2,3,5]

t=len(fx)#計算抽頭數(shù)

#step1:

M=round(math.log(N/(2*n),2)*(t+l))#計算每位平均需要校驗方程數(shù)

#step2:

defS(p,t):

ift==l:

returnp

returnp*S(p,t-l)+(l-p)*(l-S(p,t-l))

#step3:

defQ(p,M,h):

q=0

foriinrange(h,M+l):

q=comb(M,i)*(p*pow(S(p,t),i)*pow((l-S(p,t)),M-i)+(l-p)*pow(l-S(p,t),i)*pow(S(p,t),M-i))+q

returnq

#print(Q(p,M,1))

ttstep4:

defV(p,M,h):

v=0

foriinrange(h,M+1):

v=comb(M,i)*(p*pow(S(p,i),i)*pow((l-S(p,i))?M-i))+v

E

returnv

defT(p,M,h):

t=V(p,M,h)/Q(p,M,h)

ttstep5:

deffindllmax():

count=0

forhinranged,M+l):

ifQ(p,M,h)*N>=n:

count=count+l

returncount

defcalR():

r=(1-T(p,M,h=findHmax()))*n

#step6

#從python下標為n的位的連續(xù)60個數(shù)組

defProspective(n,zn_temp2):

point=n

ttprint(zn_temp)

zn_temp=zntemp2.copy()

whilepoint:

temp=zn_temp[3]"zn_temp[12]-zn_temp[22]"zn_temp[24]"zn_temp[43]~zn_temp[50]zn_temp[56]zn_t

emp[59]

zn_temp.insert(0,temp)

zntemp.pop()

point=point-l

print(zntemp)

returnzn_temp

#step7

defLastCheck(zn_temp):

count_right=0

point=60

whilepoint<=4000000-1:

temp=zn_temp[0]zn_temp[4]zn_temp[13]zn_temp[23]zntemp[25]zn_temp[44]

zn_temp[51]"zn_temp[57]

zn_temp.append(temp)

zn_temp=zn_temp[1:]

if(zn_temp[59]==zn[point]):

count_right=count_right+l

point=point+l

36

returncount_right/(4000000-60)

#step8

defChangeHaming(zn_temp):

flag二1

count=0

zn_temp_copy=zn_temp.copy()

whileflagandcount<=59:

Pros_zn=Prospective(12800,zn_temp_copy)

result=LastCheck(Proszn)

print(result,count)

print('開頭序列'，Pros_zn)

print('選取序列'，zn_temp_copy)

ifresult>0.8orresult<0.70:

zn_temp_copy=zn_temp.copy()

zn_temp_copy[count]=zn_temp_copy[count]"1

count=count+l

else:

flag=0

zn=ReadFile(,data-1.txt')

ChangeHaming(zn[12800:12860])

7.3算法n實現(xiàn)代碼

importscipy

fromscipy.specialimportcomb

importsys,time

importrandom

an=[[0foriinrange(9)]foriinrange(180)]

data=[0foriinrange(4000000)]

list_p=[0foriinrange(160)]

list_S=[0foriinrange(160)]

pn=[0.75foriinrange(4000000)]的概率pn

#存儲校驗正確與錯誤的校驗方程系數(shù)

V

correctpos=[0foriinrange(160)]

errorpos=[0foriinrange(160)]

defReadFile(filename):

datafile=open(filename)

try:

text=datafile,read()

returneva1('1+text[4:T]+']')

finally:

datafile.closeO

defWriteFile(filename):

try:

file_object=open(filename,'w')

file_object.writeC0.75\n,+str(data).replaceC

finally:

file_object.close()

defCreatePoly(n):

imax=0

count=0

aO=[0,4,13,23,25,44,51,57,60]

ai=[[0,4,13,23,25,44,51,57,60]forinrange(20)]

whilepow(2,imax)*60<4000000:

forjinrange(9):

ai[imax][j]=pow(2,imax)*aO[j]

imax=imax+l

foriinrange(imax):

forjinrange(9):

offset=n-ai[i][j]

ifoffset>=0:

ifai[i][8]+offset>=4000000:

continue

forzinrange(9):

an[count][z]=ai[i][z]+offset

count=count+l

#print("成功生成T,count,〃個校驗方程')

#foriinrange(count):

#print(an[i])

returncount

#s

defCalusp_t(p,t):

ift==0:

E

returnpow(l-p,t)

return(pow(2*p-1,t-1)*(p-0.5)+0.5)

#將pi的值打表

defMakeListOfValue_p(p):

foriinrange(160):

list_p[i]=Calusp_t(p,i)

#第i個校驗式的t位相等的概率

defCalcuS(i,t):

ift==l:

returnpn[an[i-l][0]]

else:

return(l-pn[an[i-l][t-l]])*(l-CalcuS(i-l,t-l))+pn[an[i-l][t-l]]*CalcuS(i-l,t-1)

defMakeListOfValue_S():

foriinrange(160):

list_S[i]=CalcuS(i,8)

#計算第i個位置的條件概;神4

defCalcuP(h,i,M):

Y=1

foriinranged,h+1):

Y=CalcuS(correctpos[i],8)*Y

foriinrange(h+1,M+l):

Y=Y*(1-CalcuS(errorposti],8))

Y=Y*pn[i]

T=1

foriinranged,h+1):

T=(l-CalcuS(correctpos[i],8))*T

foriinrange(h+1,M+l):

T=T*CalcuS(errorpos[i],8)

T=(l-pn[i])*T

returnY/(Y+T)

#Q

defCalcuU(p,M,h):

Q=0

s=CalcuS(O,8)

foriinrange(h,M+l):

Q=comb(M,i)*(p*pow(s,i)*pow((1-s),M-i)+(1-p)*pow(1-s,i)*pow(s,M-i))+Q

returnQ

#R

defCalcuV(p,M,h):

T=0

s=CalcuS(0,8)

foriinranged,h):

T=comb(M,i)*(p*pow(s,i)*pow((l-s),M-i))+T

returnT

defCalcuT(p,M,h):

returnCalcuV(p,M,h)/CalcuU(p,M,h)

defCalcuW(p,M,h):

T=0

s=CalcuS(0,8)

foriinrange(1,h+1):

T=comb(M,i)*((1-p)*pow(1-s,i)*pow(s,M-i))+T

returnT

deffindllmax(p,M):

Hmax=0

Imax=0

foriinrange(M):

Sprint(CalcuW(p,M,i),CalcuV(p,M,i))

ifCalcuW(p,M,i)-CalcuV(p,M,i)>Imax:

Hmax=i

#print(Hmax)

return0.5*(CalcuP(Hmax,p,M)+CalcuP(Hmax+l,p,M)),Umax

defGetEvenVerify(n):

M=CreatePoly(