学习高斯消元

点击打开链接     习题  行列式

高斯消元问题类型:

用LCM 保持整型

1. 基本的高斯消元,裸模板 HDU3359

2. 开关问题,用^操作代替 -, 求x[i]时候一样用*  poj 1222 1830 1753

3. 枚举自由变元, return -1 是因为出现[0,0,0,0,a]这种情况,return 0 是唯一解,否则是有自由变元

4. 取模方程 (a1*x1+a2*x2...)%P=b在初等行变换后，每次进行消元的时候将所有值模P，最后求解回带的时候利用扩展欧几里得来对每一个ai求一个最小的可行解       

    3*a4 % P = t4，可以表示成3*a4 + K*P = t4  ,d=exgcd(a[i][i],mod,x,y); ans=ans/d*x;    poj 2947

5. 线性基问题,异或运算中求可组成的第K大数.每个数就是equ   二进制位数为var ,消元成线性基. 比如n个数字,m个线性基.组成2^(n-m)个数字.这个K的二进制形式每个位

    的权重刚好是线性基的顺序.如果有 0行,则0为第一小. 点击打开链接      SGU 275

6. 求电阻.依据KCL定理.每个节点电流流入流出相等.电路图是双向图,要从两个方面考虑.依据电流列等式,X[i]是电压.S流入为1,T流入为-1.

    电流流出为正,流入为负(相反也行),x[u][v]++,x[v][u]++,x[u][u]--,x[v][v]--; 表示u这个点流入流出和值为0电流是双向的,所以双向. HDU 5006 HDU 3976 点击打开链接 点击打开链接

    其实换个角度可以这样想,对于每个电阻其实是两个变量对应一个方程,所以n个变量只有n-1个方程(从变量的角度来看是n个),这样自然存在一个自由变量.
原方程不变,再给一个系数前加1 就表明这个系数相关的变量为0
(方程拆分多一个).其他所有人的式子其实都是这个意思.有的人是给x[]赋值,有的是a[s][n],a[t][n] 其


实是将x[]作为了增广.
如果不设定初始值,那么求出的再相减也是答案,只不过精度会出问题(如果你是先求出上三角形再回代会出问题求的是相对值,一样的)

求解线性方程的结果会出现三种情况：无解，多解和唯一解。看下图

        

当某一行出现（0,0，……0，a) 时，方程无解。因为x1*0+x2*0+……+xn*0=a（a不为0）等式无解

当从某一行开始至最后一行，出现（0,0，……0）说明方程存在自由元，则会出现多解

当矩阵的秩等于等式个数时，有唯一解

poj 1753

#include<cmath>#include<cstdio>#include<cstring>#include<iostream>#include<algorithm>using namespace std;#define inf 0x3f3f3f3f#define eps 1e-8const int MAXN=17;int N,M,T;int a[MAXN][MAXN],x[MAXN];int free_x[MAXN],x_num;int equ,var;int Gauss(int (*a)[MAXN]){    int row,col;    x_num=0;    for(row=0,col=0;row<equ && col<var;row++,col++){        int mxr = row;        for(int i=row+1;i<equ;i++)            if( fabs(a[i][col]) - fabs(a[mxr][col]) >eps )                mxr=i;        if(mxr != row)            for(int i=row;i<var+1;i++) swap(a[row][i],a[mxr][i]);///        if(fabs(a[row][col])< eps){            free_x[x_num++]=col;row--;continue;        }        for(int i=row+1;i<equ;i++) if(fabs(a[i][col])>eps){            for(int j=col;j<var+1;j++)                a[i][j] ^= a[row][j];        }    }    for(int i=row;i<equ;i++)        if(a[i][var]>0) return -1;    if(row<var) return var-row;    for(int i=var-1;i>=0;i--) if(a[i][i]){        x[i]=a[i][var];        for(int j=i+1;j<var;j++) x[i]^=(a[i][j] * x[j]);    }    return 0;}char G[4][4];int b[MAXN][MAXN];int solve(int (*a)[MAXN]){    int t=Gauss(a);    int ans=0;    if(t==-1) return inf;    if(t==0) {        for(int i=0;i<var;i++) if(x[i]) ans++;        return ans;    }    ///枚举自由变元    ans=inf;    for(int i=0;i<(1<<t);i++){        int cnt=0;        for(int j=0;j<t;j++)            if(i&(1<<j))                x[free_x[j]]=1,cnt++;            else x[free_x[j]]=0;        for(int j=var-t-1;j>=0;j-- ){            int idx;            for(idx=j;idx<var;idx++) if(a[j][idx])                break;            x[idx]=a[j][var];            for(int l=idx+1;l<var;l++)                x[idx] ^=a[j][l] * x[l];            cnt+=x[idx];        }        ans=min(ans,cnt);    }    return ans;}int main(){    equ=16,var=16;    for(int i=0;i<4;i++)       gets(G[i]);    int ans=inf;    memset(a,0,sizeof(a));    memset(b,0,sizeof(b));    for(int i=0;i<16;i++){        a[i][i]=b[i][i]=1;        if(G[i/4][i%4]=='b') a[i][16]=0,b[i][16]=1;        else a[i][16]=1,b[i][16]=0;    }    for(int i = 0;i < 4;i++)        for(int j = 0;j < 4;j++){            int t = i*4+j;            if(i > 0)b[(i-1)*4+j][t]=a[(i-1)*4+j][t] = 1;            if(i < 3)b[(i+1)*4+j][t]=a[(i+1)*4+j][t] = 1;            if(j > 0)b[i*4+j-1][t]=a[i*4+j-1][t] = 1;            if(j < 3)b[i*4+j+1][t]=a[i*4+j+1][t] = 1;        }//        for(int i=0;i<16;i++)printf("%d ",a[i][16]);cout<<endl;//        for(int i=0;i<16;i++)printf("%d ",b[i][16]);cout<<endl;    ans=solve(a);    ans=min(ans,solve(b));    if(ans==inf) printf("Impossible\n");    else printf("%d\n",ans);    return 0;}

poj 2947

#include<cmath>#include<cstdio>#include<cstring>#include<iostream>#include<algorithm>using namespace std;#define inf 0x3f3f3f3f#define eps 1e-8const int MAXN=400;int N,M,T;int a[MAXN][MAXN],x[MAXN];int free_x[MAXN],x_num;int equ,var;inline int gcd(int a,int b){    return b==0?a:gcd(b,a%b);}inline int lcm(int a,int b){    return a/gcd(a,b)*b;}int expgcd(int a,int b,int &x,int &y){    int q,tmp;    if(b==0){        q=a;x=1;y=0;        return q;    }    q=expgcd(b,a%b,x,y);    tmp=x;x=y;    y=tmp-(a/b)*y;    return q;}int Gauss(int (*a)[MAXN]){    int row,col;    x_num=0;    for(int i=0;i<=var;i++)        free_x[i]=1;    for(row=0,col=0;row<equ && col<var;row++,col++){        int mxr = row;        for(int i=row+1;i<equ;i++)            if( fabs(a[i][col]) - fabs(a[mxr][col]) >eps )                mxr=i;        if(mxr != row)            for(int i=row;i<var+1;i++) swap(a[row][i],a[mxr][i]);///        if(fabs(a[row][col])< eps){            free_x[x_num++]=col;row--;continue;        }        for(int i=row+1;i<equ;i++) if(fabs(a[i][col])>eps){            int LCM=lcm(abs(a[i][col]),abs(a[row][col]));            int ta=LCM/abs(a[i][col]);            int tb=LCM/abs(a[row][col]);            ///if(a[i][col]*a[row][col]<0) tb=-tb;            for(int j=col;j<var+1;j++)                a[i][j] = ((a[i][j]*ta-a[row][j]*tb)%7+7)%7;        }///用LCM的方法是因为要保持整型    }    for(int i=row;i<equ;i++)        if(a[i][var]>0) return -1;    if(row<var){///这样标记自由变元的方法避免了列移动        for(int k=row-1;k>=0;k--){            x_num=0;            int idx;            for(int j=0;j<var;j++)                if(a[k][j]>eps && free_x[j]) x_num++,idx=j;            if(x_num>1) continue;            int tmp=a[k][var];            for(int j=0;j<var;j++) if(a[k][j]>eps && j!=idx)                    tmp =((tmp - a[k][j]*x[j]%7)+7)%7;            x[idx]=(tmp/a[k][idx])%7;            free_x[idx]=0;        }        return var-row;    }    for(int i=var-1;i>=0;i--){        int tmp=a[i][var];        for(int j=i+1;j<var;j++)            tmp =((tmp - a[i][j]*x[j]%7)+7)%7;        while(tmp%a[i][i]) tmp+=7;        x[i]=(tmp/a[i][i])%7;        if(x[i]<3) x[i]+=7;    }    return 0;}int solve(int (*a)[MAXN]) {    int t=Gauss(a);    int ans=0;    if(t==-1) return inf;    if(t==0) {        for(int i=0;i<var;i++) if(x[i]) ans++;        return ans;    }    ///枚举自由变元    ans=inf;    for(int i=0;i<(1<<t);i++){        int cnt=0;        for(int j=0;j<t;j++)            if(i&(1<<j))                x[free_x[j]]=1,cnt++;            else x[free_x[j]]=0;        for(int j=var-t-1;j>=0;j-- ){            int idx;            for(idx=j;idx<var;idx++) if(a[j][idx])                break;            x[idx]=a[j][var];            for(int l=idx+1;l<var;l++)                x[idx] ^=a[j][l] * x[l];            cnt+=x[idx];        }        ans=min(ans,cnt);    }    return ans;}int tans(char *s){    if(s[0]=='M') return 1;    else if(s[0]=='W') return 3;    else if(s[0]=='F') return 5;    else if(s[0]=='T' && s[1]=='U') return 2;    else if(s[0]=='T' && s[1]=='H') return 4;    else if(s[0]=='S' && s[1]=='A') return 6;    else return 7;}int main() {    while(~scanf("%d%d",&N,&M) && N+M){        memset(a,0,sizeof(a));        for(int i=0;i<M;i++){            int n;            char s1[10],s2[10];            scanf("%d%s%s",&n,s1,s2);            a[i][N]=(tans(s2)-tans(s1)+8)%7 ;            for(int j=0;j<n;j++){                int t;                scanf("%d",&t);                a[i][--t]++;                a[i][t]%=7;            }        }        equ=M,var=N;        int ans=Gauss(a);        if(ans==-1) puts("Inconsistent data.");        else if(ans>0) puts("Multiple solutions.");        else {            for(int i=0;i<N;i++)                if(i==N-1) printf("%d\n",x[i]);                else printf("%d ",x[i]);        }    }    return 0;}

HDU3976

#include<stdio.h>#include<string.h>#include<math.h>#include<algorithm>const double eps = 1e-8 ;using namespace std;int T ,cas ,N,M;double mat[55][55] ;double ans[55] ;void gauss(){int row, col ;row = col = 0 ;for( ;row<N && col<N;row++,col++){int max_r = row ;for(int i=row+1;i<N;i++){if(fabs(mat[max_r][col]) < fabs(mat[i][col]) )max_r = i ;}if(max_r != row){for(int j=col;j<=N;j++){swap(mat[row][j] ,mat[max_r][j] );}}for(int i=row+1;i<N;i++){if(fabs(mat[i][col]) < eps)continue ;double a = - mat[i][col] / mat[row][col] ;for(int j=col;j<=N;j++){mat[i][j] += mat[row][j]*a ;}}}memset(ans , 0 ,sizeof(ans)) ;for(int i=row-2;i>=0;i--){double res = mat[i][N] ;for(int j=i+1;j<N;j++){res -= mat[i][j]*ans[j];}ans[i] = res / mat[i][i] ;}printf("%.2f\n",ans[0]-ans[N-1]);}int main(){scanf("%d",&T);cas = 0 ;while(T--){++cas  ;scanf("%d %d",&N,&M);memset(mat , 0 ,sizeof(mat));for(int i=0;i<M;i++){int a,b, c; scanf("%d %d %d",&a,&b,&c);a -- ; b-- ;mat[a][a] -= 1.0/(c*1.0) ;mat[b][b] -= 1.0/(c*1.0) ;mat[a][b] += 1.0/(c*1.0) ;mat[b][a] += 1.0/(c*1.0) ; }mat[0][N] = -1 ;mat[N-1][N] = 1 ;printf("Case #%d: ",cas);gauss() ;}return 0; }

HDU 5006

/**这个就是直接设定T电压1,S电压0从增广矩阵的系数直接的出**/#include<cstdio>#include<cstring>#include<algorithm>#include<iostream>using namespace std;#define N 500int n,m,C,cnt;int vis[N],q[N],ed[100010][2],p[N][N];double A[N][N];int br[10010],fa[10010];void clear(){    memset(fa,0,sizeof fa);    memset(vis,0,sizeof vis);    memset(A,0,sizeof A);    memset(p,0,sizeof p);    for (int i = 1;i <= n;i ++) fa[i] = i;    C = 0;    cnt = 0;}int find(int x){return fa[x] == x ? x : (fa[x] = find(fa[x]));}void merge(int u,int v){int U = find(u);int V = find(v);if (U == V) return;fa[U] = V;}void gauss(int n){    for (int i = 1;i <= n;i ++){        double tmp = A[i][i];        for (int j = i + 1;j <= n;j ++) if (A[j][i] < -1e-8 || A[j][i] > 1e-8){            double t = tmp / A[j][i];            for (int k = i;k <= n + 1;k ++){                A[j][k] = A[j][k] * t - A[i][k];            }            }        }        for (int i = n;i;i --){        if (A[i][i]<1e-8 && A[i][i] >-1e-8){        continue;        }        A[i][n + 1] = A[i][n + 1] / A[i][i];        double tmp = A[i][n + 1];        for (int j = i - 1;j;j --){        A[j][n + 1] -= A[j][i] * tmp;        A[j][i] = 0.0;        }    }}void bfs(int x){    int head = 0,tail = 0,now;    q[++ tail] = x;    vis[x] = 1;    memset(vis,0,sizeof vis);    while (head < tail){        now = q[++ head];        for (int i = 1;i <= C;i ++) if (!vis[i] && p[now][i]) q[++ tail] = i,vis[i] = 1;    }}int main(){    int T;    int u,v,c,s,t;    scanf("%d",&T);    //cout << T << endl;    while (T --){        scanf("%d%d",&n,&m);        scanf("%d%d",&s,&t);        clear();        for (int i = 1;i <= m;i ++) {            scanf("%d%d%d",&u,&v,&c);            if (c) ed[cnt][0] = u,ed[cnt ++][1] = v;            else merge(u,v);        }        int a;        for (int i = 1;i <= n;i ++) {            a = find(i);            if (a == i) br[a] = ++ C;        }        if (fa[s] == fa[t]) {cout << "0.000000" << endl;continue;}        for (int i = 1;i <= cnt;i ++) {            u = ed[i][0],v = ed[i][1];            if (fa[u] != fa[v]) p[br[fa[u]]][br[fa[v]]] ++,p[br[fa[v]]][br[fa[u]]] ++;        }        bfs(br[fa[s]]);        if (!vis[br[fa[t]]]) {            printf("inf\n");            continue;        }        for (int i = 1;i <= C;i ++) {            if (i == br[fa[s]]) {A[i][i] = 1,A[i][C + 1] = 1.0;continue;}            if (i == br[fa[t]]) {A[i][i] = 1,A[i][C + 1] = 0.0;continue;}///注意这里            for (int j = 1;j <= C;j ++) if (p[i][j] > 0) A[i][j] -= p[i][j],A[i][i] += p[i][j];        }        gauss(C);        int S = br[fa[s]];        double sumi = 0.0;        for (int i = 1;i <= C;i ++)if (p[S][i]) sumi += (1.0 - A[i][C + 1]) * p[S][i];        printf("%.6lf\n",1.0 / sumi);    }    return 0;}