密码体制定义

一个密码体制是满足以下条件的五元组(P,C,K,E,D):

  1. P表示所有可能的明文组成的有限集

  2. C表示所有可能的密文组成的有限集

  3. K代表密钥空间，由所有可能的密钥组成的有限集

  4. 对每一个k∈K，都存在一个加密规则e(k)∈E和相应的解密规则d(k)∈D。并且对每对e(k):P→C, d(k):C→P，满足条件：对每一个明文x∈P，均有 $d_k(e_k(x)) = x$

密码分析的方法

  • 唯密文攻击

  • 已知明文攻击

  • 选择明文攻击

  • 选择密文攻击

  • 自适应选择密文攻击

Kerchhoff假设：攻击者Oscar知道所使用的密码体制

目标是设计在Kerchhoff假设下安全的密码体制

古典密码

代换密码（substitution cipher）是明文中的每一个字符被替换成密文中的另一个字符

置换密码（permutation cipher）又称换位密码（transposition cipher），加密过程中的明文的字母保持相同 ，但顺序被打乱了

代换密码

明文字母表A：
\(Z_q=\{0,1,...,q-1\}\\\)

明文消息是一个长度为L的字母串
\(m=(m_0,m_1,...,m_{L-1}) ,m_i \in Z_q,m \in Z_q^L\\\)

密文字母表A’
\(Z_{q'} = \{0,1,...,q'-1\}\\\)

密文组长度为L’
\(c = (c_0,c_1,...,c_{L'-1}),c_i \in Z_{q'},c \in Z_{q'}^{L'}\\\)

加密变换：明文空间到密文空间的映射
\(f:m \to c,m \in M,c \in C\\\)

在一到一的映射下，存在有逆映射
\(f^{-1}(c)=f^{-1}f(m) = m,m \in M,c\in C\\\)

加密变换通常是在密钥k的控制下变化的，在f和k的作用下由
\(Z_q^L \to Z_{q'}^{L'}\\\)

的映射，称这种密码为代换密码

L = 1时称作单字母或单码代换，也称流密码。L>1时称多字母代换或多码代换，也称分组密码

一般选择q=q’，即明文和密文字母表相同，此时，f可以构造成一到一的映射，密码没有数据扩展，L=L’

L < L’，则有数据扩展，函数f为一到多的映射，成为多名或同音代换密码

L > L’，则明文数据将被压缩，函数f不可逆，常用在数据认证系统中

单表代换：在A=A’, q=q’ 和L=1时，对所有明文字母，都用一个固定的代换进行加密

多表代换：在A=A’, q=q’ 和L=1时，使用一个以上的代换表进行加密

单表代换

移位代换密码

加密变换
\(E_k(i)=(i+k) \equiv j \pmod {q},i \geq 0, j < q \\ K =\{k \mid 0 \leq k \leq q\}\\\)

k=0时为恒等变换

解密变换
\(D(j)=E_{q-k}(j) \equiv j+q-k \equiv i+k-k \equiv i \pmod {q}\)

仿射密码

令 $P=E=Z_{26}$且$K={(a,b)\in Z_{26} \times Z_{26}:gcd(a,26)=1}$

对任意的$K=(a,b) \in K,x、y \in Z_{26}$定义加密变换为
\(e_k(x) = (ax+b) \bmod 26\\\)

解密变换为
\(d_k(y)=a^{-1}(y-b) \bmod 26\\\)

a = 1时就得到移位密码，q=26时可能的密钥数为（26*12）-1=311个

单字母代换的密码都无法抗击字母频率攻击

对抗频率分析：多名或同音代替密码、多表代替密码、多字母代替密码

多表代换

明文字母序列$m=m_1m_2…$

代换序列$\pi =(\pi_1,\pi_2,…)$

密文字母序列$c=E_k(m)=\pi(m)=\pi_1(m_1)\pi_2(m_2)…$

非周期多表代换密码，$\pi$为非周期的无限序列。对每个明文字母都采用不同的代换表进行加密，称作一次一密钥

周期多表代换密码，$\pi$为周期序列，重复地使用

当代换序列为1的时候就退化为单表代换

维吉尼亚密码

设m是一个正整数。定义P=E=K=$(Z_{26})^m$，对任意的密钥K=$(k_1,k_2,…k_m)$定义
\(e_k(x_1,x_2,...,x_m)=(x_1+k_2,x_2+k_2,...,x_m+k_m)\\ d_k(y_1,y_2,...,y_m)=(y_1-k_1,y_2-k_2,...,y_m-k_m)\\\)

以上所有的运算都是在$Z_{26}$上进行

维吉尼亚密码分析：

  • 首先确定密钥的长度m，使用kakasi测试法和重合指数法
  • 然后再使用统计特性来进行分析
kakasi测试法

密文中的重复性可以暗示出密钥的长度，如果两个相同明文序列之间的距离是密钥长度的整数倍，那么产生的密文序列也是相同的

搜索长度至少为3的相同的密文段，记下它们之间的距离，$d_1,d_2,…$密文长度m则是它们最大公因数的因子

重合指数法

设x为英语文本串,字母a,b,c,…z出现的概率分别为$p_0,p_1,…,p_{25}$，那么我们期望
\(I_c(x) \approx \sum_{i=0}^{25}p_i^2 =0.065\\\)

如果密文的重合指数较低，那么就可能是多表代换密码。维吉尼亚密码将密文分行，每行是单表代换密码。猜测密钥长度时，计算每一行的重合指数，若均趋近0.065，则很可能就确定了密钥的长度。

在单表替代时，明文的字母被其他字母代替，但不影响文本的统计属性，即加密后密文的重合指数仍不变

确定具体的密钥

令$1 \leq i \leq m，f_0,f_1,…f_{25}$分别表示串$y_i$中字母A，B，…，Z出现的频率。再令$n’ = n/m$表示串$y_i$的长度，则26个字母在$y_i$中的概率分布为：
\(\frac {f_0}{n'},\frac{f_1}{n'},...,\frac{f_{25}}{n'}\\\)

考虑到子串$y_i$是由对应的待加密的明文子集中的字母移动$k_i$个位置所得的。因此，移位后的概率分布
\(\frac{f_{k_i}}{n'},...,\frac{f_{25+k_i}}{n'}\\\)

应该近似于统计得到的各字母的概率分布$p_0,p_1,…p_{25}$

假设$0 \leq g \leq 25$，定义数值
\(M_g = \sum_{i=0}^{25}\frac{p_if_{i+g}}{n'}\\\)

如果$g=k_i$，类似前面重合指数的讨论，应有
\(M_g \approx \sum_{i=0}^{25}p_i^2=0.065\)

希尔密码

基本思想：将m个明文字母通过线性变换，将它们转换为m个密文字母。解密只需做一次逆变换即可

密钥K={$Z_{26}$上的$m\times m$可逆矩阵}，明文x=$(x_1,…x_m)$与密文y=$(y_1,…y_m)$均为m维向量

加密
\((y_1,y_2,...y_m)=(x_1,x_2,...x_m) \begin{pmatrix} k_{1,1}&k_{1,2}&\cdots&k_{1,m}\\ k_{2,1}&k_{2,2}&\cdots&k_{2,m}\\ \vdots&\vdots&\ddots&\vdots\\ k_{m,1}&k_{m,2}&\cdots&k_{m,m} \end{pmatrix}\\\)

解密$x = yK^{-1}$

希尔密码分析

线性变换的安全性很弱，易被已知明文攻击击破

对于一个$m\times m$的希尔密码，假定有m个明密文对，明文密文长度都为m，可以得到$m\times m$的方阵X和Y，有
\(Y = XK\\ K=X^{-1}Y\)

置换密码

令m为一正整数，P=C=$Z_{26}^m$，K由所有定义在集合{1,2,…,m}上的置换组成，对任意的密钥（即置换）$\pi$，定义
\(e_\pi(x_1,x_2,...x_m)=(x_{\pi(1)},x_{\pi(2)},...,x_{\pi(m)})\\ d_\pi(y_1,y_2,...,y_m)=(y_{\pi^{-1}(1)},y_{\pi^{-1}(2)},...,y_{\pi^{-1}(m)}\\\) 其中$\pi^{-1}$为置换$\pi$的逆置换

置换密码实质上是希尔密码的特例，置换可以使用一个矩阵来表示，置换矩阵为
\(K_\pi =(k_{ij})_{m\times m}\\ k_{ij}=\begin {cases} 1&\text{if } j=\pi(i)\\ 0&\text{otherwise} \end {cases}\\\)

如置换$\pi=(3\quad5\quad 1\quad 6\quad 4\quad 2)$，对应的置换矩阵为
\(K_\pi= \begin{bmatrix} 0&0&1&0&0&0\\ 0&0&0&0&1&0\\ 1&0&0&0&0&0\\ 0&0&0&0&0&1\\ 0&0&0&1&0&0\\ 0&1&0&0&0&0 \end{bmatrix}\\\)

置换密码不能抗击已知明密文攻击

密码体制的分类

依据对信息元素的处理方式可以分为：序列（流）密码、分组（块）密码

实验

维吉尼亚密码分析

首先使用kakasi分析法找出密钥长度，接下来根据密钥长度对密文进行分组，每一组密文对应同一个单个字符的密钥，对每一组进行重合指数法的分析，找出这一组对应的单个字符。由此，我们可以得到整串密钥。最后，再使用这一串密钥对密文进行解密。

主要函数及功能

void kakasi(string cipher, int length =3)

使用kakasi分析法找出长度为length的相同的字符串之间的距离的最大公因数。

int find_gcd(vector sequence)

寻找一系列整数的最可能的最大公因数

int gcd(int a, int b)

使用辗转相除法计算两个数的最大公因数

string find_key(string cipher, int key_len)

将密文cipher根据密钥长度key_len进行分组，每一组密文调用find_single_key找出本组对应的单独的密钥（单个字符），最终返回整个密钥（字符串）

string find_single_key(string sub_cipher)

寻找本组密文对应的密钥，首先统计每个英文字母在本组密文中出现的频率，使用cipher_frequency数组存储，同时列出已知的各个英文字母频率，使用standard_frequency数组存储。设本组密码移位量为g，则明文a对应的密文就为g，根据重合指数法，我们需要计算M_g的值。遍历g从0到25，取所计算得各M_g值中与0.065相差最小的那个移位量g作为本组的密钥，并作为结果返回

void decrypt(string cipher, string key)

在已知密文和密钥的情况下进行解密。因为密文是明文通过密钥移位得到的，因此，只要使用减法移回去，注意模26以及ascii值的转换即可。

代码

#include<iostream>
#include<fstream>
#include<map>
#include<vector>
#include<string>
#define debug 0
using namespace std;

float standard_frequency[26] = {0.082,0.015,0.028,0.043,0.127,0.022,0.02,0.061,0.07,0.002,0.008,0.04,0.024,0.067,0.075,0.019,0.001,0.06,0.063,0.091,0.028,0.01,0.023,0.001,0.02,0.001};
char find_single_key(string sub_cipher){
    double diff = 10;
    int possible_g;
    float M_g = 0;
    float cipher_frequency[26]={0};
    for(int i = 0; i < sub_cipher.size();i++) 
        cipher_frequency[sub_cipher[i]-'A'] += 1;
    for(int i = 0; i < 26 ;i++){
        cipher_frequency[i] = cipher_frequency[i]/sub_cipher.size();
    }
    for(int g = 0; g < 26; g++){
        M_g = 0;
        for(int j = 0; j < 26; j++){
            M_g += (standard_frequency[j]*cipher_frequency[(j+g)%26]);
        }
        if( abs(M_g - 0.065) < diff){
            diff = abs(M_g - 0.065);
            possible_g = g;
        }
    }
    return possible_g+'A';
}
string find_key(string cipher,int key_len){
    string key = "";
    string sub_cipher = "";
    for(int i = 0; i < key_len; i++){
        sub_cipher = "";
        for(int j = 0; j < cipher.size(); j++){
            if(j % key_len == i){
                sub_cipher += cipher[j];
            }
        }
        key += find_single_key(sub_cipher);
        //cout<<"group"<<i<<" subcipher:"<<endl<<sub_cipher<<endl;
    }
    return key;
}

//suppose a > b, calculate gcd(a,b)
int gcd(int a, int b){
    if(b == 0)
        return a;
    return gcd(b,a%b);
}
int find_gcd(vector<int> sequence){
    map<int,int> gcd_times;
    int temp_gcd;
    int temp_times;
    for(int i = 0; i < sequence.size(); i++){
        for(int j = 0; j < sequence.size();j++){
            if(i == j) continue;
            if(sequence[i]>sequence[j])
                temp_gcd = gcd(sequence[i],sequence[j]);
            else
                temp_gcd = gcd(sequence[j],sequence[i]);
            auto iter = gcd_times.find(temp_gcd);
            if(iter != gcd_times.end()){
                temp_times = gcd_times[temp_gcd];
                gcd_times.erase(iter);
                gcd_times.insert(make_pair(temp_gcd,temp_times+1));
            }
            else{
                gcd_times.insert(make_pair(temp_gcd,1));
            }
        }
    }
    int res, max_appear_times = -1;
    for(auto it = gcd_times.begin(); it != gcd_times.end(); it++){
        if(it->second > max_appear_times){
            res = it->first;
            max_appear_times = it->second;
        }
    }
    return res;
}
void decrypt(string cipher, string key){
    fstream output("plain.txt");
    string plain = "";
    if(debug){
        cout<<"cipher to decrypt is:"<<endl<<cipher<<endl;
        cout<<"key is"<<key<<endl;
    }
    int group = key.size();
    for(int i = 0; i < cipher.size(); i++){
        plain += ((cipher[i]-'A'-(key[i%group]-'A')+26)%26)+'A';
    }
    output << plain;
    output.close();
    cout<<"Plaintext:"<<endl;
    cout<<plain<<endl;
}
//use kakasi test to find out the same part in the cipher
int kakasi(string cipher, int length = 3){
    map<string,int> substr_fre;
    map<string,vector<int>> substr_pos;
    string check_substr;
    int times;
    for(int i = 0; i < cipher.size(); i++){
        check_substr = cipher.substr(i,length);
        auto iter = substr_fre.find(check_substr);
        if(iter == substr_fre.end()){
            vector<int> temp;
            temp.push_back(i);
            substr_fre.insert(make_pair(check_substr,1));
            substr_pos.insert(make_pair(check_substr,temp));
        }
        else{
            times = iter->second;
            substr_fre.erase(iter);
            substr_fre.insert(make_pair(check_substr,times+1));
            substr_pos[check_substr].push_back(i);
        }
    }
    string most = "";
    int appear = -1;
    for(auto it = substr_fre.begin();it != substr_fre.end();it++){
        if(it->second > appear){
            appear = it->second;
            most = it->first;
        }
    }
    vector<int> most_appear = substr_pos[most];
    vector<int> gap;
    for(int i = 0; i < most_appear.size()-1;i++){
        gap.push_back(most_appear[i+1]-most_appear[i]);
    }
    return find_gcd(gap);
}
int main(){
    fstream input;
    input.open("ciphertext.txt");
    string cipher;
    input >> cipher;
    input.close();
    cout<<"Read cipher finished."<<endl;
    int key_length = kakasi(cipher,3);
    cout<<"Key length is "<<key_length<<endl;
    string key = find_key(cipher,key_length);
    cout<<"Key is "<<key<<endl;
    decrypt(cipher,key);
}

仿射希尔密码分析

求解仿射希尔密码的密钥在知道足够长度明密文对的情况下即是解一个由两个方程组成的线性方程组。在本题中，矩阵的运算是在模26下进行的。除法可以转化为乘法。除以某个数即是乘以这个数的逆。因此，还需要能够对某个数求模逆。

具体的求解过程在下面的函数里。

关键函数及功能

int gcdEx(int a, int b, int* x, int* y)

扩展欧几里得算法，找出a，b的最大公因数，并计算出ax+by = gcd(a,b) 中的x和y

int modinv(int ele, int m)

调用上述扩展欧几里得算法求出ele模m的逆

void known_plain_attack_hill(stirng plain, string cipher, int m = M)

已知明文攻击破解仿射希尔密码。

要求解一个m=M的仿射希尔密码，我们需要两个线性方程，且需要知道长度为2 *M *M的明文。开始时判断一下输入的明文密文对是否足够长。

接下来使用2个M阶方阵矩阵表示明文，2 个M阶方阵表示密文，记为p1,p2,c1,c2，密钥L和B同样用M阶方阵表示，记为keyL，keyB。我们有
\(c1 = p1 \times keyL + keyB \\ c2 = p2 \times keyL + keyB\)

用下式减上式，有
\((c2-c1) = (p2-p1) \times keyL\)

记c2-c1的结果为cSub，p2-p1的结果为pSub，我们有
\(cSub = pSub \times keyL\\ keyL = cSub \times pSub^{-1}\)

那么，我们只需要计算出pSub的逆，就能够很快得到keyL的值了

这里使用初等行变换计算矩阵pSub在模26运算下的逆矩阵。

新创建一个M*2M的矩阵，由pSub和M阶单位阵拼接而成。在此例中M=3，即
\(\left[ \begin{matrix} a_1 & a_2 & a_3 & 1 & 0 & 0\\ b_1 & b_2 & b_3 & 0 & 1 & 0\\ c_1 & c_2 & c_3 & 0 & 0 & 1 \end{matrix} \right ]\)

将左边的Psub变成一个上三角矩阵,使用双重for循环第一重循环遍历每一行，将每一行的对应元素变为1，比如说在第i行，调用modinv函数求出第i行第i列元素在模26运算下的逆，不妨记为ele_inv，我们的目标是将第i行第i个元素变成1，因此需要使用初等行变换的标量乘运算，那么进入第二重循环，遍历当行的每一个元素，每一个元素都乘以ele_inv，并模26。循环完毕后，我们得到下面这样一个矩阵
\(\left[ \begin{matrix} 1 & a'_2 & a'_3 & a'_4 & a'_5 & a'_6\\ 0 & 1 & b'_3 & b'_4 & b'_5 & b'_6\\ 0 & 0 & 1 & c'_4 & c'_5 & c'_6 \end{matrix} \right ]\)

接下来再将左边变成一个单位阵。解决方法是使用三重for循环，第一重循环自下而上遍历每一行，比如第i行；第二重循环从当行的第M-1个元素开始，向左遍历到第i-1个元素，目的是依次将这些元素变为0，因此找出需要减去的值，记为iden_coe，假设现在的位置为j，进入第三重循环，对于当行的每一个元素，使用初等行变换，用矩阵的第j行乘以（-iden_coe）再加到这一行来，因为第j行第j列的值正好为1，且左边方阵其余位置均为0，那么右边得到的便是pSub的逆矩阵了。
\(\left[ \begin{matrix} 1 & 0 & 0 & a''_4 & a''_5 & a''_6\\ 0 & 1 & 0 & b''_4 & b''_5 & b''_6\\ 0 & 0 & 1 & c''_4 & c''_5 & c''_6 \end{matrix} \right ]\)

我们使用cSub乘以pSub的逆矩阵，就得到了KeyL 的值。

这样keyB的值也非常好计算了，只要将keyL的值带入下式其中之一即可
\(keyB = c1 - p1 \times keyL\\keyB = c2 - p2 \times keyL\)

代码

#include<iostream>
#define M 3
#define debug 1
using namespace std;
int gcdEx(int a, int b, int *x, int *y) {
    if(b==0){
        *x = 1,*y = 0;
        return a;
    }
    else{
        int r = gcdEx(b, a%b, x, y);
        int t = *x;
        *x = *y;
        *y = t - a/b * *y;
        return r;
    }
}
int modinv(int ele, int m){
    int* x = new int[1];
    int* y = new int[1];
    gcdEx(m,ele,x,y);
    int res = *y;
    while(res < 0) res += 26;
    delete[] x;
    delete[] y;
    return res;
}
void print_M_col_matrix(int matrix[][M], int row){
    for(int i = 0; i < row; i++){
        for(int j = 0; j < M; j++){
            cout<<matrix[i][j]<<" ";
        }
        cout<<endl;
    }
}
void print_2M_col_matrix(int matrix[][2*M], int row){
    for(int i = 0; i < row; i++){
        for(int j = 0; j < 2*M; j++){
            cout<<matrix[i][j]<<" ";
        }
        cout<<endl;
    }
}
//plain and cipher should be the same length and can divide m
void known_plain_attack_hill(string plain, string cipher, int m = M){
    if((plain.size() != cipher.size())||plain.size()<2*M*M){
        cout<<"Length of plaintext and ciphertext not correspond, or plain length not long enough"<<endl;
        return;
    }
    int keyL[M][M];
    int p1[M][M], p2[M][M], c1[M][M],c2[M][M];
    for(int i = 0; i < M; i++){
        for(int j = 0; j < M; j++){
            p1[i][j] = plain[i*M+j]-'a';
            c1[i][j] = cipher[i*M+j]-'a';
            p2[i][j] = plain[M*M+i*M+j]-'a';
            c2[i][j] = cipher[M*M+i*M+j]-'a';
        }
    }
    //now we have two equations:
    //c1 = p1*keyL+b
    //c2 = p2*keyL+b
    //now we calculate c2-c1 to solve the keyL
    int c_sub[M][M];
    int p_sub[M][M];
    for(int i = 0; i < M; i++){
        for(int j = 0; j < M; j++){
            c_sub[i][j] = c2[i][j] - c1[i][j];
            p_sub[i][j] = p2[i][j] - p1[i][j];
        }
    }
    //we need to calculate the inverse matrix of p_sub
    //c_sub = p_sub * keyL
    //use elementary row operation
    int row_operation[M][2*M];
    for(int i = 0; i < M; i++){
        for(int j = 0; j < 2*M; j++){
            if(j < M){
                row_operation[i][j] = p_sub[i][j];
            }
            else{
                if(i == j-M)
                    row_operation[i][j] = 1;
                else
                    row_operation[i][j] = 0;
            }
        }
    }
    if(debug){
        cout<<"orignial matrix"<<endl;
        print_2M_col_matrix(row_operation,M);
    }
    //turn it into a uppper triangular matrix
    //for every row
    int ele_inv;
    int sub_coe;
    for(int i = 0; i < M; i++){
        //turn its diagonal entry to 1
        ele_inv = modinv(row_operation[i][i],26);
        //div_coe = row_operation[i][i];
        for(int k = 0; k < 2*M; k++)
            row_operation[i][k] = (row_operation[i][k]*ele_inv)%26;
            //row_operation[i][k] /= div_coe;

        //then turn the entry under diagonal entry to 0,
        //by subtracting the ith row by a coefficient
        for(int j = i+1; j < M; j++){
            sub_coe = row_operation[j][i];
            for(int r = 0; r < 2*M; r++){
                row_operation[j][r] -= sub_coe*row_operation[i][r];
                while(row_operation[j][r] < 0){
                    row_operation[j][r] += 26;
                }
                row_operation[j][r] %= 26;      
            }
                
        }
    }
    if(debug){
        cout<<"triangular matrix"<<endl;
        print_2M_col_matrix(row_operation,M);
    }
    //then change other non-diagonal element to 0
    //for every row, begin from the last row, 
    int iden_coe;
    for(int i = M-2; i > -1; i--){
        //check its position that needs to turn to 0
        for(int j = M-1; j > i; j--){
            iden_coe = row_operation[i][j];
            for(int k = 0; k < 2*M; k++ ){
                row_operation[i][k] -= row_operation[j][k]*iden_coe;
                while(row_operation[i][k] < 0) row_operation[i][k] += 26;
                row_operation[i][k] %= 26;
            }
        }
    }
    if(debug){
        cout<<"Identity matrix"<<endl;
        print_2M_col_matrix(row_operation,M);
    }
    int p_sub_inv[M][M];
    for(int i = 0; i < M; i++){
        for(int j = 0; j < M; j++){
            p_sub_inv[i][j] = row_operation[i][j+M];
        }
    }
    cout<<"p_sub_inv:"<<endl;
    print_M_col_matrix(p_sub_inv,3);
    cout<<"c_sub"<<endl;
    for(int i = 0; i < 3;i++){
        for(int j = 0; j < 3; j++)
            cout<<c_sub[i][j]<<" ";
        cout<<endl;
    }
    //now calculate p_sub_inv*c_sub to get keyL
    for(int i = 0; i < M; i++){
        for(int j = 0; j < M; j++){
            keyL[i][j] = 0;
            for(int k = 0; k < M; k++){
                keyL[i][j] += p_sub_inv[i][k]*c_sub[k][j];
                while(keyL[i][j] < 0) keyL[i][j] += 26;
                keyL[i][j] %= 26;
            }
        }
    }
    cout<<"keyL is:"<<endl; 
    print_M_col_matrix(keyL,M);
    //now calculate keyB, c1 = p1*keyL + keyB
    //keyB = c1 - p1*keyL
    int keyB[M][M]={0};
    int temp_sum;
    for(int i = 0; i < M; i++){
        for(int j = 0; j < M; j++){
            temp_sum = 0;
            for(int k = 0; k < M; k++){
                temp_sum += p1[i][k]*keyL[k][j];
            }
            keyB[i][j] = c1[i][j] - temp_sum;
            while(keyB[i][j] < 0) keyB[i][j] += 26;
            keyB[i][j] %= 26;
        }
    }
    cout<<"keyB is:"<<endl;
    print_M_col_matrix(keyB,M);
}
int main(){
    string plain = "adisplayedequation";
    string cipher = "dsrmsioplxljbzullm";
    known_plain_attack_hill(plain,cipher);
}