高斯消元模板

xiaoxiao2021-02-28  99

入门去看:英雄哪里出来和诚叙(两位大佬)。

理解后的部分修改的两个板子:整型和浮点型。

//整型版 #include <string.h> #include <cstdio> #define ll long long using namespace std; struct GAUSS { const static int mod = 1e9+7; const static int maxn = 1005; int r, c, free_x, cnt; //cnt表示x取值的值域 ll matrix[maxn][maxn]; ll x[maxn]; ll ans; void init(int _r, int _c, int _cnt) { r = _r, c = _c, cnt = _cnt; memset(matrix, 0, sizeof matrix); } //构造增广矩阵(系数, b) void Link(int r, int c, ll val) { matrix[r][c] = val; } ll _abs(ll x) { return x < 0? -x: x; } ll qpow(int base, int n) { ll ans = 1; for(int i = 1; i <= n; ++i) ans *= base; return ans; } void ex_gcd(ll a, ll b, ll &d, ll &x, ll &y) { if(b == 0) { x = 1, y = 0, d = a; return; } ex_gcd(b, a%b, d, y, x); y -= a/b*x; } void swap_row(int a, int b) { for(int i = 0; i <= c; ++i) { ll tmp = matrix[a][i]; matrix[a][i] = matrix[b][i]; matrix[b][i] = tmp; } } void swap_col(int a, int b) { for(int i = 0; i < r; ++i) { ll tmp = matrix[i][a]; matrix[i][a] = matrix[i][b]; matrix[i][b] = tmp; } } bool gauss() { int row, col, maxrow, i, j; for(row = 0, col = 0; row < r && col < c; ++row) { maxrow = row; for(i = row+1; i < r; ++i) { if(_abs(matrix[i][col]) > _abs(matrix[maxrow][col])) maxrow = i; } if(matrix[maxrow][col] == 0) { ++col; --row; continue; } if(maxrow != row) swap_row(row, maxrow); for(i = row+1; i < r; ++i) { if(matrix[i][col]) { ll tmp = matrix[i][col]; //matrix[i][col]必须暂时保存,在首次计算时就会被改变 for(j = col; j <= c; ++j) { matrix[i][j] = (matrix[i][j]*matrix[row][col]-matrix[row][j]*tmp)%mod; if(matrix[i][j] < 0) matrix[i][j] += mod; } } } ++col; } for(i = row; i < r; ++i) if(matrix[i][c]) return false; //无解 free_x = c-row; //转化为上三角 for(i = 0; i < r && i < c; ++i) { if(matrix[i][i] == 0) { for(j = i+1; j < c; ++j)//虽然此处条件是判断是否<c,但实际上要么一定能找到系数矩阵中一列,要么连带b向量中全找不到 if(matrix[i][j]) { swap_col(i, j); break; } } } ans = qpow(cnt, free_x); return true; } //如果存在唯一解,下面计算每个未知数的值 void calc_x() { for(int i = c-1; i >= 0; --i) { ll tmp = matrix[i][c]; for(int j = i+1; j < c; ++j) tmp = (tmp-matrix[i][j]*x[j]+mod)%mod; ll D, X, Y; ex_gcd(matrix[i][i], mod, D, X, Y); X = (X%mod+mod)%mod; x[i] = X*tmp%mod; } } } AC; void solve() { int r, c, cnt; scanf("%d %d %d", &r, &c, &cnt); AC.init(r, c, cnt); for(int i = 0; i < r; ++i) for(int j = 0; j <= c; ++j) //输入增广矩阵(系数, b) { int x; scanf("%lld", &x); AC.Link(i, j, x); } if(!AC.gauss()) printf("No Answer.\n"); { printf("%lld\n", AC.ans); if(AC.ans == 1) { AC.calc_x(); for(int i = 0; i < c; ++i) printf("%lld ", AC.x[i]); printf("\n"); } } } int main() { solve(); return 0; }

//浮点数版 #include <string.h> #include <cstdio> #include <cmath> #define LL long long using namespace std; struct GAUSS { const static int maxn = 105; const static double eps = 1e-12; int r, c, free_x; double matrix[maxn][maxn]; double x[maxn]; void init(int _r, int _c) { r = _r, c = _c; //浮点数也可以这么清零 memset(matrix, 0, sizeof matrix); } //构造增广矩阵(系数, b) void Link(int r, int c, double val) { matrix[r][c] = val; } inline bool sgn(double x) //判断绝对值是否小于极小数, = 0代表为极小数 { return (x > eps) - (x < -eps); } void swap_row(int a, int b) { for(int i = 0; i <= c; ++i) { double tmp = matrix[a][i]; matrix[a][i] = matrix[b][i]; matrix[b][i] = tmp; } } void swap_col(int a, int b) { for(int i = 0; i < r; ++i) { double tmp = matrix[i][a]; matrix[i][a] = matrix[i][b]; matrix[i][b] = tmp; } } bool gauss() { int row, col, maxrow, i, j; for(row = 0, col = 0; row < r && col < c; ++row) { maxrow = row; for(i = row+1; i < r; ++i) //先进行此处判断,因为绝对值小于极小数的数都被认为是0 { if(fabs(matrix[i][col]) > fabs(matrix[maxrow][col])) maxrow = i; } if(sgn(matrix[maxrow][col]) == 0) { ++col; --row; continue; } if(maxrow != row) swap_row(row, maxrow); for(i = row+1; i < r; ++i) { if(sgn(matrix[i][col]) != 0) { double tmp = matrix[i][col]; for(j = col; j <= c; ++j) { matrix[i][j] -= tmp/matrix[row][col]*matrix[row][j]; } } } ++col; } for(i = row; i < r; ++i) if(sgn(matrix[row][c]) != 0) return false; //无解 free_x = c-row; for(i = 0; i < r && i < c; ++i) { if(sgn(matrix[i][i]) == 0) { for(j = i+1; j < c; ++j) if(sgn(matrix[i][j]) != 0) { swap_col(i, j); break; } } } return true; } void calc_x() { for(int i = c-1; i >= 0; --i) { double tmp = matrix[i][c]; for(int j = i+1; j < c; ++j) tmp = tmp-matrix[i][j]*x[j]; x[i] = tmp/matrix[i][i]; } } } AC; void solve() { int r, c; scanf("%d %d", &r, &c); AC.init(r, c); for(int i = 0; i < r; ++i) for(int j = 0; j < c; ++j) { double x; scanf("%lf", &x); AC.Link(i, j, x); } AC.gauss(); AC.calc_x(); for(int i = 0; i < c-1; ++i) printf("%.2f ", AC.x[i]+AC.eps); //避免-0.0的输出 printf("%.2f\n", AC.x[c-1]+AC.eps); } int main() { solve(); return 0; }

附:-0.0的产生是因为被用来表示一个可以四舍五入为零的负数,或者是一个从负方向上趋近于零的数。所以浮点数运算时最终结果需要加一个eps避免-0.0。

Go on.
转载请注明原文地址: https://www.6miu.com/read-46232.html

最新回复(0)