C++实现伽罗华域生成及四则运算（三）

上一篇文章介绍了用C++实现伽罗华域$GF(2^n)$数的次幂、矩阵数乘、方阵次幂、求行列式、求伴随矩阵、逆矩阵的函数，本篇继续，对求行列式的函数进行改进，最后结合幻方4阶矩阵探索一下求逆矩阵的其他方法。

原来的求行列式代码

/*方阵行列式,r是阶数
代数余子式法
*/
template<uint8_t n>
uint16_t GFM<n>::delta(uint16_t* A, uint16_t r) {if (r == 0) throw TypeException("矩阵阶数出错");if (r == 1) return A[0];if (r == 2) return this->multi(A[0],A[3]) ^ this->multi(A[1],A[2]); //对角线法// 按第0行展开代数余子式uint16_t* B = (uint16_t*)malloc(sizeof(uint16_t) * (r-1) * (r-1));uint16_t res = 0;for (uint16_t omit = 0; omit < r; omit++) { //指定列号，求指定列的代数余子式for (uint16_t i = 1, k = 0; i < r; i++) { //从第1行开始读取，写入B[k]for (uint16_t j = 0; j < r; j++) { //逐列读取if (j == omit) continue; //跳过指定的列B[k++] = A[i*r+j];}}res ^= this->multi(A[omit],delta(B, r-1)); //按A的第0行代数余子式展开 }free(B); //释放内存return res;
}

是利用嵌套展开代数余子式的方法，uint16_t* B = (uint16_t*)malloc(sizeof(uint16_t) * (r-1) * (r-1));申请内存，用于存放比方阵A低1阶的子方阵，而对于该子方阵，只用再次调用delta函数就可以了，理解上应该没有困难，但仔细一想，太费内存了。
假如A是16阶，则B就是15阶，B的子方阵就是14阶，依次类推，直到最后到达2阶if (r == 2) return this->multi(A[0],A[3]) ^ this->multi(A[1],A[2]);才结束，这前后一共申请内存$(16^2+15^2+\cdots +2^2)\ast sizeof(uint16\_t)=2982Bytes$，足足约3K。

改进后的求行列式代码

/*方阵行列式,r是阶数
上三角对角线法
*/
template<uint8_t n>
uint16_t GFM<n>::delta(uint16_t* A, uint16_t r) {if (r == 0) throw TypeException("矩阵阶数出错");if (r == 1) return A[0];//复制Auint16_t* B = (uint16_t*)malloc(sizeof(uint16_t) * r * r);memcpy(B, A, sizeof(uint16_t) * r * r); uint16_t res = 1;//将方阵转化成上三角形for (uint16_t i = 0; i < r-1; i++) { //指定与第i行异或if (B[i*r+i] == 0) { //第i行i列若已经是0，则要与不为0的下一行交换一整行，行列式变号，但GF(2^n)域内特殊，不用变号for (uint16_t j = i+1; j < r; j++) {if (B[j*r+i] == 0) continue;uint16_t temp[r];memcpy(temp, B+i*r, sizeof(uint16_t)*r);memcpy(B+i*r, B+j*r, sizeof(uint16_t)*r);memcpy(B+j*r, temp, sizeof(uint16_t)*r);}if (B[i*r+i] == 0) return 0; //全部查询一遍，结果仍为0，则结果必然为0};for (uint16_t j = i+1; j < r; j++) { //指定当前参与变换的行if (B[j*r+i] == 0) continue; //已经是0，刚好不用乘uint16_t row_times = this->multi(this->adver(B[i*r+i]),B[j*r+i]); //计算商=B[j][i]/B[i][i]for (uint16_t k = i; k < r; k++) { //当前行加上第i行乘row_times，行列式不变，省略前i个0B[j*r+k] ^= this->multi(B[i*r+k],row_times);}}}//计算对角线乘积for (uint16_t i = 0; i < r; i++) {res = this->multi(B[i*r+i], res);}
cout << res << endl;free(B); //释放内存return res;
}

看着有些长，但它不嵌套，采用的是变换上三角方法，最后直接把对角线作乘积就得出了其行列式值。这里用到2个定理：

• 一行的倍数加到另一行，行列式不变
  $$\left|\begin{array}{ccc}a& b& c\\ d& e& f\\ g& h& i\end{array}\right|=\left|\begin{array}{ccc}a& b& c\\ d+a\ast k& e+b\ast k& f+c\ast k\\ g& h& i\end{array}\right|$$
• 交换两行，行列式变号
  $$\left|\begin{array}{ccc}a& b& c\\ d& e& f\\ g& h& i\end{array}\right|=-\left|\begin{array}{ccc}d& e& f\\ a& b& c\\ g& h& i\end{array}\right|$$

GF(256)域幻方4阶矩阵求逆

幻方4阶矩阵形如下：
$$\left[\begin{array}{cccc}a& b& c& d\\ d& a& b& c\\ c& d& a& b\\ b& c& d& a\end{array}\right]或\left[\begin{array}{cccc}a& b& c& d\\ b& c& d& a\\ c& d& a& b\\ d& a& b& c\end{array}\right]$$
其特点是每行依次左移或右移1位，在很多加密算法里都比较常见。
那么如何求它的逆矩阵呢？用前2点的方法肯定可以，现在我介绍第3种方法，先看推导：
$${\left[\begin{array}{cccc}a& b& c& d\\ d& a& b& c\\ c& d& a& b\\ b& c& d& a\end{array}\right]}^4={\left(\left[\begin{array}{cccc}a& b& c& d\\ d& a& b& c\\ c& d& a& b\\ b& c& d& a\end{array}\right]\times \left[\begin{array}{cccc}a& b& c& d\\ d& a& b& c\\ c& d& a& b\\ b& c& d& a\end{array}\right]\right)}^2$$

$$={\left[\begin{array}{cccc}a\ast a+b\ast d+c\ast c+d\ast b& a\ast b+b\ast a+c\ast d+d\ast c& a\ast c+b\ast b+c\ast a+d\ast d& a\ast d+b\ast c+c\ast b+d\ast a\\ d\ast a+a\ast d+b\ast c+c\ast b& d\ast b+a\ast a+b\ast d+c\ast c& d\ast c+a\ast b+b\ast a+c\ast d& d\ast d+a\ast c+b\ast b+c\ast a\\ c\ast a+d\ast d+a\ast c+b\ast b& c\ast b+d\ast a+a\ast d+b\ast c& c\ast c+d\ast b+a\ast a+b\ast d& c\ast d+d\ast c+a\ast b+b\ast a\\ b\ast a+c\ast d+d\ast c+a\ast b& b\ast b+c\ast a+d\ast d+a\ast c& b\ast c+c\ast b+d\ast a+a\ast d& b\ast d+c\ast c+d\ast b+a\ast a\end{array}\right]}^2$$

$$={\left[\begin{array}{cccc}a\ast a+c\ast c& 0& b\ast b+d\ast d& 0\\ 0& a\ast a+c\ast c& 0& b\ast b+d\ast d\\ b\ast b+d\ast d& 0& a\ast a+c\ast c& 0\\ 0& b\ast b+d\ast d& 0& a\ast a+c\ast c\end{array}\right]}^2$$

$$=\left[\begin{array}{cccc}(a\ast a+c\ast c{)}^2+(b\ast b+d\ast d{)}^2& 0& (a\ast a+c\ast c)(b\ast b+d\ast d)+(b\ast b+d\ast d)(a\ast a+c\ast c)& 0\\ 0& (a\ast a+c\ast c{)}^2+(b\ast b+d\ast d{)}^2& 0& (a\ast a+c\ast c)(b\ast b+d\ast d)+(b\ast b+d\ast d)(a\ast a+c\ast c)\\ (b\ast b+d\ast d)(a\ast a+c\ast c)+(a\ast a+c\ast c)(b\ast b+d\ast d)& 0& (b\ast b+d\ast d{)}^2+(a\ast a+c\ast c{)}^2& 0\\ 0& (b\ast b+d\ast d)(a\ast a+c\ast c)+(a\ast a+c\ast c)(b\ast b+d\ast d)& 0& (b\ast b+d\ast d{)}^2+(a\ast a+c\ast c{)}^2\end{array}\right]$$

$$=\left[\begin{array}{cccc}(a\ast a+c\ast c{)}^2+(b\ast b+d\ast d{)}^2& 0& 0& 0\\ 0& (a\ast a+c\ast c{)}^2+(b\ast b+d\ast d{)}^2& 0& 0\\ 0& 0& (a\ast a+c\ast c{)}^2+(b\ast b+d\ast d{)}^2& 0\\ 0& 0& 0& (a\ast a+c\ast c{)}^2+(b\ast b+d\ast d{)}^2\end{array}\right]$$

$$=\left[\begin{array}{cccc}{a}^4+c^4+b^4+d^4& 0& 0& 0\\ 0& a^4+c^4+b^4+d^4& 0& 0\\ 0& 0& a^4+c^4+b^4+d^4& 0\\ 0& 0& 0& a^4+c^4+b^4+d^4\end{array}\right]=A$$

令$x=a^4+b^4+c^4+d^4$因此
$$A÷x=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

由此可以推断：
$${\left[\begin{array}{cccc}a& b& c& d\\ d& a& b& c\\ c& d& a& b\\ b& c& d& a\end{array}\right]}^3÷x与\left[\begin{array}{cccc}a& b& c& d\\ d& a& b& c\\ c& d& a& b\\ b& c& d& a\end{array}\right]互为逆矩阵$$

证：
$${\left[\begin{array}{cccc}a& b& c& d\\ d& a& b& c\\ c& d& a& b\\ b& c& d& a\end{array}\right]}^3÷x=\left[\begin{array}{cccc}a\ast a+c\ast c& 0& b\ast b+d\ast d& 0\\ 0& a\ast a+c\ast c& 0& b\ast b+d\ast d\\ b\ast b+d\ast d& 0& a\ast a+c\ast c& 0\\ 0& b\ast b+d\ast d& 0& a\ast a+c\ast c\end{array}\right]\times \left[\begin{array}{cccc}a& b& c& d\\ d& a& b& c\\ c& d& a& b\\ b& c& d& a\end{array}\right]÷x$$
$$=\left[\begin{array}{cccc}{a}^3+a{c}^2+b^{2}c+c{d}^2& a^{2}b+b{c}^2+b^{2}d+d^3& a^{2}c+c^3+a{b}^2+a{d}^2& a^{2}d+c^{2}d+b^3+b{d}^2\\ a^{2}d+c^{2}d+b^3+b{d}^2& a^3+a{c}^2+b^{2}c+c{d}^2& a^{2}b+b{c}^2+b^{2}d+d^3& a^{2}c+c^3+a{b}^2+a{d}^2\\ a^{2}c+c^3+a{b}^2+a{d}^2& a^{2}d+c^{2}d+b^3+b{d}^2& a^3+a{c}^2+b^{2}c+c{d}^2& a^{2}b+b{c}^2+b^{2}d+d^3\\ a^{2}b+b{c}^2+b^{2}d+d^3& a^{2}c+c^3+a{b}^2+a{d}^2& a^{2}d+c^{2}d+b^3+b{d}^2& a^3+a{c}^2+b^{2}c+c{d}^2\end{array}\right]÷x$$
$$=\left[\begin{array}{cccc}\frac{a^3+a{c}^2+b^{2}c+c{d}^2}{x}& \frac{a^{2}b+b{c}^2+b^{2}d+d^3}{x}& \frac{a^{2}c+c^3+a{b}^2+a{d}^2}{x}& \frac{a^{2}d+c^{2}d+b^3+b{d}^2}{x}\\ \frac{a^{2}d+c^{2}d+b^3+b{d}^2}{x}& \frac{a^3+a{c}^2+b^{2}c+c{d}^2}{x}& \frac{a^{2}b+b{c}^2+b^{2}d+d^3}{x}& \frac{a^{2}c+c^3+a{b}^2+a{d}^2}{x}\\ \frac{a^{2}c+c^3+a{b}^2+a{d}^2}{x}& \frac{a^{2}d+c^{2}d+b^3+b{d}^2}{x}& \frac{a^3+a{c}^2+b^{2}c+c{d}^2}{x}& \frac{a^{2}b+b{c}^2+b^{2}d+d^3}{x}\\ \frac{a^{2}b+b{c}^2+b^{2}d+d^3}{x}& \frac{a^{2}c+c^3+a{b}^2+a{d}^2}{x}& \frac{a^{2}d+c^{2}d+b^3+b{d}^2}{x}& \frac{a^3+a{c}^2+b^{2}c+c{d}^2}{x}\end{array}\right]$$

$${\left[\begin{array}{cccc}a& b& c& d\\ d& a& b& c\\ c& d& a& b\\ b& c& d& a\end{array}\right]}^3÷x\times \left[\begin{array}{cccc}a& b& c& d\\ d& a& b& c\\ c& d& a& b\\ b& c& d& a\end{array}\right]$$

$$=\left[\begin{array}{cccc}\frac{a^4+d^4+c^4+b^4}{x}& \frac{0}{x}& \frac{0}{x}& \frac{0}{x}\\ \frac{0}{x}& \frac{a^4+d^4+c^4+b^4}{x}& \frac{0}{x}& \frac{0}{x}\\ \frac{0}{x}& \frac{0}{x}& \frac{a^4+d^4+c^4+b^4}{x}& \frac{0}{x}\\ \frac{0}{x}& \frac{0}{x}& \frac{0}{x}& \frac{a^4+d^4+c^4+b^4}{x}\end{array}\right]=E$$
证毕。

测试代码如下：

//测试4阶幻方矩阵GFM<8> gf8;// uint16_t a[16] = {0xe,0x9,0x8,0x7, 0x7,0xe,0x9,0x8, 0x8,0x7,0xe,0x9, 0x9,0x8,0x7,0xe};uint16_t a[16] = {0xe,0x9,0x8,0x7, 0x9,0x8,0x7,0xe, 0x8,0x7,0xe,0x9, 0x7,0xe,0x9,0x8};uint16_t res[16];uint16_t q[16];uint16_t x = gf8.adver(gf8.calc("0xe^4+0x9^4+0x8^4+0x7^4"));gf8.matrixPow(q, a, 4, 3);gf8.matrixNumMulti(q, q, 4, 4, x);gf8.matrixMulti(res, a, 4, 4, q, 4, 4);cout << "验证:" << endl;for (int i = 0; i < 4; i++) {for (int j = 0; j < 4; j++){cout << hex << setw(2) << res[i*4+j] << " ";}cout << endl;}

致敬

向伽罗瓦、莱布尼茨等老一辈数学家致敬！