高斯消元模板

来源:互联网 发布:中国的地位知乎 编辑:程序博客网 时间:2024/06/08 08:46

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

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

//整型版#include <string.h>#include <cstdio>#define ll long longusing 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 longusing 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.
原创粉丝点击