如果商用密码产品认证中遇到问题，欢迎加微信symmrz或13720098215沟通。

3.8 块内最大游程检测方法

块内最大游程检测方法如下图。

以下实现代码供参考。

#get blockFrequency param, follow GM/T 0005-2021 appendix A
def getlongestRunOfOnesM(len):if len < 6272:m = 8K = 3pi = [0.2148, 0.3672,0.2305,0.1875]Vparam = [1,2,3,4]elif len < 750000:m =  128K = 5pi = [0.1174,0.2430,0.2494,0.1752,0.1027,0.1124]Vparam = [4,5,6,7,8,9]else:m = 10000K = 6pi = [0.086632, 0.208201, 0.248419,0.193913,0.121458,0.068011,0.073366]Vparam = [10,11,12,13,14,15,16]return {'m':m, 'K':K, 'pi':pi, 'vparam':Vparam}  #run'1' and run'0' is decided by input
def longestRunOfOnes(epsilon:bitarray, runtype:int):n = int(len(epsilon))param = getlongestRunOfOnesM(n)K = param['K']m = param['m']pi = param['pi']vparam = param['vparam']N = n/mepsilonList = epsilon.tolist()nu = [0]*(K+1)for i in range(int(N)):v_n_obs = 0run = 0for j in range(m):if epsilonList[i*m+j] == runtype:run += 1if  run > v_n_obs:v_n_obs = runelse:run = 0if v_n_obs < vparam[0]:nu[0] += 1elif v_n_obs > vparam[K]:nu[K] += 1else:for j in range(K+1):if v_n_obs == vparam[j]:nu[j] += 1breakV_value = 0.0for i in range(K+1):V_value += ((nu[i] - N*pi[i])**2)/(N * pi[i])pvalue = scipy.special.gammaincc(K/2, V_value/2)qvalue = pvaluereturn {'p':pvalue, 'q':qvalue}

3.9 二元推导检测方法

二元推导检测方法如下图。

以下实现代码供参考。

def getK(n):if n < 100000000:k = [3,7]else:k = [3,7,15]return k
def binaryDerivate(epsilon:bitarray, k):epsilonList = epsilon.tolist()n = len(epsilonList)if n < 100:print('error n len!')for loop in range(k):for i in range(n-1-loop):epsilonList[i] = epsilonList[i]^epsilonList[i+1]del epsilonList[n - k: n]sn = epsilonList.count(1) - epsilonList.count(0)V = sn/sqrt(n-k)# print('V:', V)pvalue = math.erfc(abs(V)/sqrt(2))qvalue = math.erfc(V/sqrt(2))/2return {'p':pvalue, 'q':qvalue}

3.10 自相关检测

自相关检测方法如下图。

以下实现代码供参考。

def getselfCorrelationM(len):if len < 1000000:d = [2,8,16]elif len < 100000000:d = [1,2,8,16]else:d = [1,2,8,16,32]return d    def selfCorrelation(epsilon:bitarray, d):n = int(len(epsilon))epsilonList = epsilon.tolist()ad = 0for i in range(n-d):ad += epsilonList[i]^epsilonList[i+d]V = 2*(ad-((n-d)/2))/sqrt(n-d)pvalue = math.erfc(abs(V)/sqrt(2))qvalue = math.erfc(V/sqrt(2))/2return {'p':pvalue, 'q':qvalue}

3.11 矩阵秩检测

矩阵秩检测方法如下图。

以下实现代码供参考。

MATRIX_FORWARD_ELIMINATION  = 0
MATRIX_BACKWARD_ELIMINATION = 1def perform_elementary_row_operations(flag:int, i:int, M:int, Q:int, A:np.matrix):if flag == MATRIX_FORWARD_ELIMINATION:for j in range(i+1,M):#( j=i+1; j<M;  j++ )if A[j,i] == 1: for k in range(i,Q):#( k=i; k<Q; k++ ) A[j,k] = (A[j,k] + A[i,k]) % 2else:for j in range(i-1,-1,-1):#( j=i-1; j>=0;  j-- )if A[j,i] == 1:for k in range(Q):#k=0; k<Q; k++ )A[j,k] = (A[j,k] + A[i,k]) % 2def find_unit_element_and_swap(flag:int, i:int, M:int, Q:int, A:np.matrix):row_op=0if flag == MATRIX_FORWARD_ELIMINATION:index = i+1while (index < M) and (A[index,i] == 0): index += 1if index < M:row_op = swap_rows(i, index, Q, A)else:index = i-1while (index >= 0) and (A[index,i] == 0): index -= 1if index >= 0:row_op = swap_rows(i, index, Q, A)return row_opdef swap_rows(i:int, index:int, Q:int, A:np.matrix):A[[i,index],:] = A[[index,i],:]return 1def determine_rank(m:int, M :int, Q:int, A:np.matrix):# DETERMINE RANK, THAT IS, COUNT THE NUMBER OF NONZERO ROWS rank = mfor i in range(M):if np.sum(A[i] == 1) == 0:rank -= 1return rankdef computeRank(M:int, Q:int, matrix:np.matrix, countflag):m=min(M,Q);# FORWARD APPLICATION OF ELEMENTARY ROW OPERATIONSfor i in range(m-1):if matrix[i,i] == 1:perform_elementary_row_operations(MATRIX_FORWARD_ELIMINATION, i, M, Q, matrix)else:if find_unit_element_and_swap(MATRIX_FORWARD_ELIMINATION, i, M, Q, matrix) == 1:perform_elementary_row_operations(MATRIX_FORWARD_ELIMINATION, i, M, Q, matrix)for i in range(m-1, 0,-1):#( i=m-1; i>0; i-- ) {if matrix[i,i] == 1:perform_elementary_row_operations(MATRIX_BACKWARD_ELIMINATION, i, M, Q, matrix)else:# {    /* matrix[i][i] = 0 */if find_unit_element_and_swap(MATRIX_BACKWARD_ELIMINATION, i, M, Q, matrix) == 1:perform_elementary_row_operations(MATRIX_BACKWARD_ELIMINATION, i, M, Q, matrix)rank = determine_rank(m, M, Q, matrix)return rankdef rank(epsilon:bitarray):M = 32Q = 32MMutipleQ = M*Qn = int(len(epsilon))N = int(n/MMutipleQ)epsilonList = epsilon.tolist()R = []if N==0:print('rank N error!')return 1for i in range(N):A = np.mat([epsilonList[i*MMutipleQ + j*Q:i*MMutipleQ + j*Q + Q:1] for j in range(M)],int)R.append(computeRank(M, Q, A, i==3)) # print('R:', R)FM = R.count(M)FM_1 = R.count(M-1)V = (FM-0.2888*N)**2/(0.2888*N) + (FM_1-0.5776*N)**2/(0.5776*N) + (N-FM-FM_1-0.1336*N)**2/(0.1336*N) pvalue = scipy.special.gammaincc(1, V/2)qvalue = pvaluereturn {'p':pvalue, 'q':qvalue}