对分组密码算法进行结构性测试，采取明密文独立性测试的方法验证

#include <iostream>
#include <random>
#include <bitset>
#include <ctime>
#include <cassert>
#include <cmath>

const int maxm = 10; //最大明文组数

using std::bitset;
using std::cout;
using std::endl;
using std::mt19937_64;

// des部分
class des
{
private:
    bool Init = false;     //初始化
    bitset<64> key;        //64位密钥
    bitset<48> subKey[16]; // 16 个子密钥
    // 初始 IP 置换
    int IP[64] = {58, 50, 42, 34, 26, 18, 10, 2,
                  60, 52, 44, 36, 28, 20, 12, 4,
                  62, 54, 46, 38, 30, 22, 14, 6,
                  64, 56, 48, 40, 32, 24, 16, 8,
                  57, 49, 41, 33, 25, 17, 9, 1,
                  59, 51, 43, 35, 27, 19, 11, 3,
                  61, 53, 45, 37, 29, 21, 13, 5,
                  63, 55, 47, 39, 31, 23, 15, 7};
    // IP 逆置换
    int IP_[64] = {40, 8, 48, 16, 56, 24, 64, 32,
                   39, 7, 47, 15, 55, 23, 63, 31,
                   38, 6, 46, 14, 54, 22, 62, 30,
                   37, 5, 45, 13, 53, 21, 61, 29,
                   36, 4, 44, 12, 52, 20, 60, 28,
                   35, 3, 43, 11, 51, 19, 59, 27,
                   34, 2, 42, 10, 50, 18, 58, 26,
                   33, 1, 41, 9, 49, 17, 57, 25};
    // E - 扩展规则
    int E[48] = {32, 1, 2, 3, 4, 5,
                 4, 5, 6, 7, 8, 9,
                 8, 9, 10, 11, 12, 13,
                 12, 13, 14, 15, 16, 17,
                 16, 17, 18, 19, 20, 21,
                 20, 21, 22, 23, 24, 25,
                 24, 25, 26, 27, 28, 29,
                 28, 29, 30, 31, 32, 1};

    // 8 个 S - 盒
    int S[8][4][16] = {
        {{14, 4, 13, 1, 2, 15, 11, 8, 3, 10, 6, 12, 5, 9, 0, 7},
         {0, 15, 7, 4, 15, 2, 13, 1, 10, 6, 12, 11, 9, 5, 3, 8},
         {4, 1, 14, 8, 13, 6, 2, 11, 15, 12, 9, 7, 3, 10, 5, 0},
         {15, 12, 8, 2, 4, 9, 1, 7, 5, 11, 3, 14, 10, 0, 6, 13}},

        {{15, 1, 8, 14, 6, 11, 3, 4, 9, 7, 2, 13, 12, 0, 5, 10},
         {3, 13, 4, 7, 15, 2, 8, 14, 12, 0, 1, 10, 6, 9, 11, 5},
         {0, 14, 7, 11, 10, 4, 13, 1, 5, 8, 12, 6, 9, 3, 2, 15},
         {13, 8, 10, 1, 3, 15, 4, 2, 11, 6, 7, 12, 0, 5, 14, 9}},

        {{10, 0, 9, 14, 6, 3, 15, 5, 1, 13, 12, 7, 11, 4, 2, 8},
         {13, 7, 0, 9, 3, 4, 6, 10, 2, 8, 5, 14, 12, 11, 15, 1},
         {13, 6, 4, 9, 8, 15, 3, 0, 11, 1, 2, 12, 5, 10, 14, 7},
         {1, 10, 13, 0, 6, 9, 8, 7, 4, 15, 14, 3, 11, 5, 2, 12}},

        {{7, 13, 14, 3, 0, 6, 9, 10, 1, 2, 8, 5, 11, 12, 4, 15},
         {12, 8, 11, 5, 6, 15, 0, 3, 4, 7, 2, 12, 1, 10, 14, 9},
         {10, 6, 9, 0, 12, 11, 7, 13, 15, 1, 3, 14, 5, 2, 8, 4},
         {3, 15, 0, 6, 10, 1, 13, 8, 9, 4, 5, 11, 12, 7, 2, 14}},

        {{2, 12, 4, 1, 7, 10, 11, 6, 8, 5, 3, 15, 13, 0, 14, 9},
         {14, 11, 2, 12, 4, 7, 13, 1, 5, 0, 15, 10, 3, 9, 8, 6},
         {4, 2, 1, 11, 10, 13, 7, 8, 15, 9, 12, 5, 6, 3, 0, 14},
         {11, 8, 12, 7, 1, 14, 2, 13, 6, 15, 0, 9, 10, 4, 5, 3}},

        {{12, 1, 10, 15, 9, 2, 6, 8, 0, 13, 3, 4, 14, 7, 5, 11},
         {10, 15, 4, 2, 7, 12, 9, 5, 6, 1, 13, 14, 0, 11, 3, 8},
         {9, 14, 15, 5, 2, 8, 12, 3, 7, 0, 4, 10, 1, 13, 11, 6},
         {4, 3, 2, 12, 9, 5, 15, 10, 11, 14, 1, 7, 6, 0, 8, 13}},

        {{4, 11, 2, 14, 15, 0, 8, 13, 3, 12, 9, 7, 5, 10, 6, 1},
         {13, 0, 11, 7, 4, 9, 1, 10, 14, 3, 5, 12, 2, 15, 8, 6},
         {1, 4, 11, 13, 12, 3, 7, 14, 10, 15, 6, 8, 0, 5, 9, 2},
         {6, 11, 13, 8, 1, 4, 10, 7, 9, 5, 0, 15, 14, 2, 3, 12}},

        {{13, 2, 8, 4, 6, 15, 11, 1, 10, 9, 3, 14, 5, 0, 12, 7},
         {1, 15, 13, 8, 10, 3, 7, 4, 12, 5, 6, 11, 0, 14, 9, 2},
         {7, 11, 4, 1, 9, 12, 14, 2, 0, 6, 10, 13, 15, 3, 5, 8},
         {2, 1, 14, 7, 4, 10, 8, 13, 15, 12, 9, 0, 3, 5, 6, 11}}};

    // P - 置换
    int P[32] = {16, 7, 20, 21,
                 29, 12, 28, 17,
                 1, 15, 23, 26,
                 5, 18, 31, 10,
                 2, 8, 24, 14,
                 32, 27, 3, 9,
                 19, 13, 30, 6,
                 22, 11, 4, 25};

    // PC-1 置换
    int PC1[56] = {57, 49, 41, 33, 25, 17, 9,
                   1, 58, 50, 42, 34, 26, 18,
                   10, 2, 59, 51, 43, 35, 27,
                   19, 11, 3, 60, 52, 44, 36,
                   63, 55, 47, 39, 31, 23, 15,
                   7, 62, 54, 46, 38, 30, 22,
                   14, 6, 61, 53, 45, 37, 29,
                   21, 13, 5, 28, 20, 12, 4};

    // PC-2 压缩置换
    int PC2[48] = {14, 17, 11, 24, 1, 5,
                   3, 28, 15, 6, 21, 10,
                   23, 19, 12, 4, 26, 8,
                   16, 7, 27, 20, 13, 2,
                   41, 52, 31, 37, 47, 55,
                   30, 40, 51, 45, 33, 48,
                   44, 49, 39, 56, 34, 53,
                   46, 42, 50, 36, 29, 32};

    // 子密钥生成时，每轮左移位置数
    int loop[16] = {1, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1};

    // 左移操作
    bitset<28> LS(bitset<28> A, int shift);

    // 子密钥生成
    void generateSubKey();

    // 轮函数
    bitset<32> f(bitset<32> R, bitset<48> K);

public:
    void init(bitset<64> k); //初始化(就是使用加密密钥生成轮密钥)

    bitset<64> encrypt(bitset<64> m);

    bitset<64> decrypt(bitset<64> c);
};

bitset<28> des::LS(bitset<28> A, int shift)
{
    bitset<28> output = A;
    if (shift == 1)
    {
        for (int i = 0; i < 27; i++)
        {
            if (i - shift < 0)
                output[i - shift + 28] = A[i];
            else
                output[i] = A[i + shift];
        }
    }
    if (shift == 2)
    {
        for (int i = 0; i < 26; i++)
        {
            if (i - shift < 0)
                output[i - shift + 28] = A[i];
            else
                output[i] = A[i + shift];
        }
    }
    return output;
}
void des::generateSubKey()
{
    bitset<56> reducedKey;
    bitset<28> Ci;
    bitset<28> Di;
    bitset<48> Ki;

    // (1) 对 K 的 56 个非校验位实行置换 PC-1，得到 C0D0，其中 C0 和 D0 分别由 PC-1 置换后的前 28 位和后 28 位组成。
    for (int i = 0; i < 56; i++)
    {
        reducedKey[i] = key[PC1[i] - 1];
    }

    // (2) 计算 Ci = LSi (Ci-1) 和 Di = LSi (Di-1)，即左移
    // (3) 对 56 位的 CiDi 实行 PC-2 压缩置换，得到 48 位的 Ki。
    // (4) 如果已经得到 K16，密钥调度过程结束。
    for (int times = 0; times < 16; times++)
    {
        // 上下分开
        for (int i = 0; i < 28; i++)
        {
            Ci[i] = reducedKey[i];
        }
        for (int i = 28; i < 56; i++)
        {
            Di[i - 28] = reducedKey[i];
        }
        Ci = LS(Ci, loop[times]);
        Di = LS(Di, loop[times]);

        // 上下合并
        for (int i = 0; i < 56; i++)
        {
            if (i < 28)
            {
                reducedKey[i] = Ci[i];
            }
            else
            {
                reducedKey[i] = Di[i - 28];
            }
        }

        // PC-2 置换，得到 48 位的 subKey
        for (int i = 0; i < 48; i++)
        {
            Ki[i] = reducedKey[PC2[i] - 1];
        }

        // 添加到子密钥集合中
        subKey[times] = Ki;
    }
}

bitset<32> des::f(bitset<32> R, bitset<48> K)
{
    bitset<48> Ei;
    bitset<32> result;

    // (1) 将长度为 32 位的 R 作 E - 扩展，成为 48 为的串
    for (int i = 0; i < 48; i++)
    {
        Ei[i] = R[E[i] - 1];
    }

    // (2) 将 Ei 和长度为 48 位的子密钥 K 作 48 位二进制串按位异或运算。
    for (int i = 0; i < 48; i++)
    {
        Ei[i] = Ei[i] ^ K[i];
    }

    // (3) 将 (2) 得到的结果平均分成 8 个分组（每个分组长度为 6），各个分组分别经过 8 个不同的 S - 盒进行 6-4 转换，得到 8 个长度分别为 4 位的分组。
    int j = 0; //result 的下标
    for (int i = 0; i < 48; i += 6)
    {
        //input: b1b2b3b4b5b6
        //n = (b1b6) 行号
        //m = (b2b3b4b5) 列号
        int n = Ei[i] * 2 + Ei[i + 5];
        int m = Ei[i + 1] * 8 + Ei[i + 2] * 4 + Ei[i + 3] * 2 + Ei[i + 4];
        int value = S[i / 6][n][m];

        // (4) 将 (3) 得到的分组结果顺序连接得到长度为 32 位的串
        bitset<4> newValue(value); //int to binary
        result[j] = newValue[3];
        result[j + 1] = newValue[2];
        result[j + 2] = newValue[1];
        result[j + 3] = newValue[0];
        j += 4;
    }

    bitset<32> output = result;
    // (5) 将 (4) 的结果经过 P - 置换，得到的结果作为论函数 f 的最终 32 位输出。
    for (int i = 0; i < 32; i++)
    {
        output[i] = result[P[i] - 1];
    }

    return output;
}

void des::init(bitset<64> k)
{
    key = k;
    generateSubKey();
    Init = true;
}

bitset<64> des::encrypt(bitset<64> m)
{
    bitset<64> M0;
    bitset<32> Li;
    bitset<32> Ri;
    bitset<32> temp;
    bitset<64> c;

    assert(Init); //断言是否初始化

    // 初始置换 IP
    for (int i = 0; i < 64; i++)
    {
        M0[i] = m[IP[i] - 1];
    }

    for (int i = 0; i < 32; i++)
    {
        Li[i] = M0[i];
    }
    for (int i = 32; i < 64; i++)
    {
        Ri[i - 32] = M0[i];
    }

    // 16 轮迭代 T
    for (int times = 0; times < 16; times++)
    {
        temp = Ri;
        Ri = Li ^ f(Ri, subKey[times]);
        Li = temp;
    }

    // 逆置换 IP^-1
    for (int i = 0; i < 32; i++)
    {
        c[i] = Ri[i];
    }
    for (int i = 32; i < 64; i++)
    {
        c[i] = Li[i - 32];
    }

    M0 = c; // 使用中间变量存储 c
    for (int i = 0; i < 64; i++)
    {
        c[i] = M0[IP_[i] - 1];
    }

    return c;
}

bitset<64> des::decrypt(bitset<64> c)
{
    bitset<64> M0;
    bitset<32> Li;
    bitset<32> Ri;
    bitset<32> temp;
    bitset<64> m;

    assert(Init); //断言是否初始化

    // 初始置换 IP
    for (int i = 0; i < 64; i++)
    {
        M0[i] = c[IP[i] - 1];
    }

    for (int i = 0; i < 32; i++)
    {
        Li[i] = M0[i];
    }
    for (int i = 32; i < 64; i++)
    {
        Ri[i - 32] = M0[i];
    }

    // 16 轮迭代 T (逆序)
    for (int times = 15; times >= 0; times--)
    {
        temp = Ri;
        Ri = Li ^ f(Ri, subKey[times]);
        Li = temp;
    }

    // 逆置换 IP^-1
    for (int i = 0; i < 32; i++)
    {
        m[i] = Ri[i];
    }
    for (int i = 32; i < 64; i++)
    {
        m[i] = Li[i - 32];
    }

    M0 = m;
    for (int i = 0; i < 64; i++)
    {
        m[i] = M0[IP_[i] - 1];
    }

    return m;
}

int main()
{
    bitset<64> m[maxm], c[maxm], d[maxm];
    uint64_t seed = time(NULL);
    mt19937_64 rnd(seed); //初始化PRNG
    bitset<64> key = rnd();
    cout << "1)随机密钥key:" << key << endl;
    cout << "2)随机生成"<<maxm<<"个明文"<<endl;
    for (int i = 0; i < maxm; i++)
    {
        m[i] = rnd();
        cout << "随机明文m" << i+1 << ":" << m[i] << endl;
    }
    cout << "3)开始加密(des)" << endl;
    des des1;
    des1.init(key);
    for (int i = 0; i < maxm; i++)
    {
        c[i] = des1.encrypt(m[i]);
        cout << "加密得到c" << i+1 << ":" << c[i] << endl;
    }
    cout << "4)计算明密文距离:" << endl;
    for (int i = 0; i < maxm; i++)
    {
        d[i] = m[i] ^ c[i];
        cout << "明密文距离d" << i+1 << ":" << d[i] << endl;
    }
    cout << "5)计算di汉明重量" << endl;
    for (int i = 0; i < maxm; i++)
    {
        cout << "汉明重量w" << i+1 << "汉明重量:" << d[i].count() << endl;
    }
    cout << "6)统计汉明重量分组vi"<<endl;
    int wset[6]={0,28,30,32,34,64};
    int v[5]={0};
    for (int i = 0; i < maxm; i++)
    {
        int cnt=d[i].count();
        if(cnt<=wset[1]) v[0]++;
        else if(cnt<=wset[2]) v[1]++;
        else if(cnt<=wset[3]) v[2]++;
        else if(cnt<=wset[4]) v[3]++;
        else v[4]++;
    }
    for(int i=0 ;i<5;i++){
        cout << "分组v"<<i+1<<"从"<<((i==0)?wset[i]:wset[i]+1)<<"到"<<wset[i+1]<<":	"<<v[i]<<endl;
    }
    double obs=0;
    double pi[5]={0.190866,0.163124,0.292020,0.163124,0.190866};//pi_i的值
    for(int i=0;i<5;i++){
        obs+=pow((double)v[i]-maxm*pi[i],2)/(maxm*pi[i]);
    }
    cout <<"obs:"<<obs<<endl;
    return 0;
}