随机数检测（三）

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

3.8 块内最大游程检测方法

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

以下实现代码供参考。

#get blockFrequency param, follow GM/T 0005-2021 appendix A
def getlongestRunOfOnesM(len):
    if len < 6272:
        m = 8
        K = 3
        pi = [0.2148, 0.3672,0.2305,0.1875]
        Vparam = [1,2,3,4]
    elif len < 750000:
        m =  128
        K = 5
        pi = [0.1174,0.2430,0.2494,0.1752,0.1027,0.1124]
        Vparam = [4,5,6,7,8,9]
    else:
        m = 10000
        K = 6
        pi = [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/m

    epsilonList = epsilon.tolist()

    nu = [0]*(K+1)
    for i in range(int(N)):
        v_n_obs = 0
        run = 0
        for j in range(m):
            if epsilonList[i*m+j] == runtype:
                run += 1
                if  run > v_n_obs:
                    v_n_obs = run
            else:
                run = 0
        
        if v_n_obs < vparam[0]:
            nu[0] += 1
        elif v_n_obs > vparam[K]:
            nu[K] += 1
        else:
            for j in range(K+1):
                if v_n_obs == vparam[j]:
                    nu[j] += 1
                    break

    V_value = 0.0
    for 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 = pvalue
    return {'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))/2
    return {'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 = 0
    for 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))/2
    return {'p':pvalue, 'q':qvalue}

3.11 矩阵秩检测

矩阵秩检测方法如下图。

以下实现代码供参考。

MATRIX_FORWARD_ELIMINATION  = 0
MATRIX_BACKWARD_ELIMINATION = 1

def 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]) % 2
    
    else:
        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]) % 2
    

def find_unit_element_and_swap(flag:int, i:int, M:int, Q:int, A:np.matrix):
    row_op=0
    if flag == MATRIX_FORWARD_ELIMINATION:
        index = i+1
        while (index < M) and (A[index,i] == 0): 
            index += 1
        if index < M:
            row_op = swap_rows(i, index, Q, A)
    
    else:
        index = i-1
        while (index >= 0) and (A[index,i] == 0): 
            index -= 1
        if index >= 0:
            row_op = swap_rows(i, index, Q, A)
    return row_op

def swap_rows(i:int, index:int, Q:int, A:np.matrix):
    A[[i,index],:] = A[[index,i],:]
    return 1

def determine_rank(m:int, M :int, Q:int, A:np.matrix):
    # DETERMINE RANK, THAT IS, COUNT THE NUMBER OF NONZERO ROWS 
    rank = m
    for i in range(M):
        if np.sum(A[i] == 1) == 0:
            rank -= 1
    return rank

def computeRank(M:int, Q:int, matrix:np.matrix, countflag):
    m=min(M,Q);
    # FORWARD APPLICATION OF ELEMENTARY ROW OPERATIONS
    for 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 rank

def rank(epsilon:bitarray):
    M = 32
    Q = 32
    MMutipleQ = M*Q
    
    n = int(len(epsilon))
    N = int(n/MMutipleQ)
    epsilonList = epsilon.tolist()
    R = []
    if N==0:
        print('rank N error!')
        return 1

    for 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 = pvalue
    return {'p':pvalue, 'q':qvalue}