最近信息安全作业要写hill密码,涉及到求矩阵的逆。不想用一个函数就用那些乱七八糟的库,大概找了找又没现成的,只能自己撸了。
写的时候才发现线代都忘光了QAQ
1 //Gauss-Jordan Elimination method to get inverse matrix 2 template <typename T=double> 3 std::vector<std::vector<double>> matrixInversion(const std::vector<std::vector<T>>& a) 4 { 5 const static double eps = 1e-6; 6 std::vector<std::vector<double>> augma; //augmented matrix 7 for (auto& vec : a) 8 { 9 std::vector<double> tmp; 10 std::transform(vec.begin(), vec.end(), std::back_inserter(tmp), [](const T interger) 11 { 12 return static_cast<double>(interger); 13 }); 14 augma.push_back(std::move(tmp)); 15 } 16 int rank = a.size(); 17 //init 18 for (int i = 0; i < rank; ++i) 19 { 20 augma[i].resize(rank * 2); 21 augma[i][rank + i] = 1; 22 } 23 for (int i = 0; i < rank; ++i) 24 { 25 if (::abs(augma[i][i]) < eps) 26 { 27 int j; 28 for (j = i + 1; j < rank; ++j) 29 { 30 if (::abs(augma[j][i]) > eps) break; 31 } 32 if (j == rank) return {}; 33 for (int r = i; r < 2 * rank; ++r) 34 { 35 augma[i][r] += augma[j][r]; 36 } 37 } 38 double ep = augma[i][i]; 39 for (int r = i; r < 2 * rank; ++r) 40 { 41 augma[i][r] /= ep; 42 } 43 44 for (int j = i + 1; j < rank; ++j) 45 { 46 double e = -1 * (augma[j][i] / augma[i][i]); 47 for (int r = i; r < 2 * rank; ++r) 48 { 49 augma[j][r] += e * augma[i][r]; 50 } 51 } 52 } 53 for (int i = rank - 1; i >= 0; --i) 54 { 55 for (int j = i - 1; j >= 0; --j) 56 { 57 double e = -1 * (augma[j][i] / augma[i][i]); 58 for (int r = i; r < 2 * rank; ++r) 59 { 60 augma[j][r] += e * augma[i][r]; 61 } 62 } 63 } 64 65 for (auto& vec : augma) 66 { 67 vec.erase(vec.begin(), vec.begin() + rank); 68 } 69 70 return augma; 71 }