随机数检测（四）

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

3.12 累加和检测方法

累加和检测方法如下图。

以下实现代码供参考。

def cusum(epsilon:bitarray, direction):

    n = int(len(epsilon))
    X = [0]*n
    S = [0]*n
    epsilonList = epsilon.tolist()
    for i in range(n):
        X[i]=2*epsilonList[i]-1
    if direction == 0:#forward cumum
        
        S[0] = X[0]
        SMax = abs(S[0])
        for i in range(1,n):
            S[i] = S[i-1]+X[i]
            if abs(S[i]) > SMax:
                SMax = abs(S[i])
    else:#backward cusunm
        S[0] = X[n-1]
        SMax = abs(S[n-1])
        for i in range(1,n):
            S[i] = S[i-1]+X[n-1-i]
            if abs(S[i]) > SMax:
                SMax = abs(S[i])
    sum1 = 0.0
    for i in range(int((1-(n/SMax))/4), int(((n/SMax)-1)/4)):
        sum1 += norm.cdf((4*i + 1)*SMax/sqrt(n)) - norm.cdf((4*i - 1)*SMax/sqrt(n))
    sum2 = 0.0
    for i in range(int((-3-(n/SMax))/4), int(((n/SMax)-1)/4)):
        sum2 += norm.cdf((4*i + 3)*SMax/sqrt(n)) - norm.cdf((4*i + 1)*SMax/sqrt(n))
    pvalue =  1.0 - sum1 + sum2
    qvalue = pvalue
    return {'p':pvalue, 'q':qvalue}

3.13 近似熵检测方法

近似熵检测方法如下图。


以下实现代码供参考。

def getselfCorrelationM(len):
    if len < 100000000:
        m = [2,5]
    else:
        m = [5,7]
    return m    

def alterOffsetBit(m, offset, listSub:list):
    ret = []
    listSub[offset] = 0
    if offset == (m-1):
        ret.append(listSub.copy())
    else: 
        ret += (alterOffsetBit(m, offset+1, listSub))

    listSub[offset] = 1
    if offset == (m-1):
        ret.append(listSub.copy())
    else: 
        ret += (alterOffsetBit(m, offset+1, listSub))
    return ret
    

def seedSequence(m):
    sequencelist = [0]*m
    return alterOffsetBit(m, 0, sequencelist)

def getPhyN(epsilonList:list, mlen):
    n = len(epsilonList)
    prefixM = epsilonList[:mlen-1:1]
    epsilonList += prefixM

    Vm = getV(epsilonList, mlen)
    C = [Vm[i]/n for i in range(2**mlen)]
    
    phy_mlist = [C[i]*math.log(C[i],math.e) for i in range(2**mlen)]
    return sum(phy_mlist)
    
def getV(epsilonList:list, mlen):
    seedList = seedSequence(mlen)
    V = [0] * (2**mlen)
    for i in range(len(epsilonList)-mlen+1):
        seedindex = seedList.index(epsilonList[i:i+mlen:1])
        V[seedindex] += 1
    return V

def approximateEntropy(epsilon:bitarray,m):
    n = int(len(epsilon))
    epsilonList = epsilon.tolist()
    
    apen_m =  getPhyN(epsilonList, m) - getPhyN(epsilonList, m+1) 
    V = 2.0*n*(math.log(2,math.e) - apen_m)
    pvalue = scipy.special.gammaincc(2**(m-1), V/2)
    qvalue = pvalue
    return {'p':pvalue, 'q':qvalue}

3.14 线性复杂度检测

线性复杂度检测方法如下图。


以下实现代码供参考。

def linearComplexity(len):
    if len < 100000000:
        m = [500,1000]
    else:
        m = [5000]
    return m    

def Berlekamp_Massey_algorithm(sequence):
    N = len(sequence)
    s = sequence[:]

    for k in range(N):
        if s[k] == 1:
            break
    f = set([k + 1, 0])
    l = k + 1

    g = set([0])
    a = k
    b = 0

    for n in range(k + 1, N):
        d = 0
        for ele in f:
            d ^= s[ele + n - l]

        if d == 0:
            b += 1
        else:
            if 2 * l > n:
                f ^= set([a - b + ele for ele in g])
                b += 1
            else:
                temp = f.copy()
                f = set([b - a + ele for ele in f]) ^ g
                l = n + 1 - l
                g = temp
                a = b
                b = n - l + 1
   
    return l#(print_poly(f), l)

def getV(epsilonList:list, N:int, M:int):
    V = [0]*7
    mean = M/2.0 + (9.0+1)/36.0 - 1.0/(2**M) * (M/3.0 + 2.0/9.0)
    temp = [0]*M
    for ii in range(N):
        B_ = [0]*M
        C = [0]*M
        T = [0]*M
        P = [0]*M
        # print('loop N', ii)
        L = 0
        m = -1
        d = 0
        C[0] = 1
        B_[0] = 1
        
        for N_ in range(M):
            d = epsilonList[ii*M+N_]
            for i in range(1,L+1):
                d += C[i] * epsilonList[ii*M+N_-i]
            d = d%2
            if d == 1:
                T = C.copy()
                P = temp.copy()
                for j in range(M):
                    if B_[j] == 1:
                        P[j+N_-m] = 1
                for i in range(M):
                    C[i] = (C[i] + P[i])%2
                if L <= int(N_/2) :
                    L = N_ + 1 - L
                    m = N_
                    B_ = T.copy()
        
        T_ = (L - mean) + 2.0/9.0

        if T_ <= -2.5:
            V[0]+=1
        elif  T_ > -2.5 and T_ <= -1.5:
            V[1]+=1
        elif  T_ > -1.5 and T_ <= -0.5:
            V[2]+=1
        elif  T_ > -0.5 and T_ <= 0.5:
            V[3]+=1
        elif  T_ > 0.5 and T_ <= 1.5:
            V[4]+=1
        elif  T_ > 1.5 and T_ <= 2.5:
            V[5]+=1
        else:
            V[6]+=1
    return V
    
def linearComplexity(epsilon:bitarray,m):

    n = len(epsilon)
    N = int(len(epsilon)/m)
    epsilonList = epsilon.tolist()
    
    L = [ Berlekamp_Massey_algorithm(epsilonList[i*m:(i+1)*m]) for i in range(N)]
    mu = (m/2) + (9 + (-1)**(m+1))/36 - (m/3 + 2/9)/(2**m)
    T = [((-1)**m)*(L[i] - mu) + 2/9 for i in range(N)]
    V = [0]*7
    for i in range(N):
        if T[i] <= -2.5:
            V[0] += 1
        elif T[i] <= -1.5 and T[i] > -2.5:
            V[1] += 1
        elif T[i] <= -0.5 and T[i] > -1.5:
            V[2] += 1
        elif T[i] <= 0.5 and T[i] > -0.5:
            V[3] += 1
        elif T[i] <= 1.5 and T[i] > 0.5:
            V[4] += 1 
        elif T[i] <= 2.5 and T[i] > 1.5:
            V[5] += 1
        elif T[i] > 2.5:
            V[6] += 1
    pi = [0.010417, 0.03125, 0.12500, 0.5000, 0.25000, 0.06250, 0.020833]
    resultV = [((V[i] - N*pi[i])**2)/(N*pi[i]) for i in range(7)]
    pvalue = scipy.special.gammaincc(3, sum(resultV)/2)
    qvalue = pvalue
    return {'p':pvalue, 'q':qvalue}

3.15 Maurer通用统计检测

Maurer通用统计检测方法如下图。


以下实现代码供参考。

expected_value = [0, 0, 0, 0, 0, 0, 5.2177052, 6.1962507, 7.1836656,
                8.1764248, 9.1723243, 10.170032, 11.168765,
                12.168070, 13.167693, 14.167488, 15.167379]
variance = [ 0, 0, 0, 0, 0, 0, 2.954, 3.125, 3.238, 3.311, 3.356, 3.384,
                3.401, 3.410, 3.416, 3.419, 3.421 ]

def getLQ(n:int):
    
    if  n >= 1000000:
        return (7,1280)
    else:
        return (7,1280)
    

def universal(epsilon:bitarray, L):

    n = len(epsilon)
    # (L,Q) = getLQ(n)
    Q = 10*(2**L)
    K = int(n/L)-Q
    p = 2**L
    c = 0.7 - 0.8/L + (4 + 32/L)*(K**(-3/L))/15
    sigma = c * sqrt(variance[L]/K)
    sqrt2 = sqrt(2)
    sum = 0.0
    T = [0]*p
    for i in range(1,Q+1):    
        decRep = 0
        for j in range(L):
            decRep = (decRep << 1) + epsilon[(i-1)*L + j]

        T[decRep] = i
    for i in range(Q+1, Q+K+1):
        decRep = 0
        for j in range(L):
            decRep = (decRep << 1) + epsilon[(i-1)*L + j]
        sum += np.log2(i - T[decRep])
        T[decRep] = i
    
    phi =(sum/K)
    V = (phi-expected_value[L])/sigma
    pvalue = math.erfc(abs(V)/sqrt2)
    qvalue = math.erfc(V/sqrt2)/2
    return {'p':pvalue, 'q':qvalue}

3.16 离散傅里叶检测方法

离散傅里叶检测方法如下图。


以下实现代码供参考。

def discreteFourierTransform(epsilon:bitarray):

    n = len(epsilon)
    X = [2*epsilon[i]-1 for i in range(n)]
    f = np.fft.fft(X)
    m = [abs(i) for i in f[0:int(n/2)]]
    count = 0
    upperBound = sqrt(2.995732274*n)
    for i in range(int(n/2)):
        if m[i] < upperBound:
            count +=1
    N_l = count
    N_o = int(0.95*n/2.0)
    d = (N_l - N_o)/sqrt(n/3.8*0.95*0.05)

    pvalue = math.erfc(abs(d)/sqrt(2))
    qvalue = math.erfc(d/sqrt(2))/2