高斯消元

来源:互联网 发布:宾馆网络电视怎么打开 编辑:程序博客网 时间:2024/06/05 20:22

介绍这个算法:
1.通过一遍消元可以使得矩阵成为上三角矩阵(主对角线以下都是0的矩阵)。
2.再消一遍可以得到对角矩阵。
我们通过分析第一次消元的过程就能知道第二次怎么消除。
步骤:
1.首先我们枚举要用哪行消除。
2.如果在正常情况下,即在消元的过程中,i,aii0,那么我们可以一直用第i个方程去消去[i+1,n]的所有系数。
3.然而有可能在消去的过程中,第i个方程的第i个系数为0,那么我们就找到随意一个k>i使得aki0,把第k行和第i行对换,然后用当前的第i行消去以后的方程即可。(当前的第i行就是之前找到的那个第k行)。
4.如果找不到aki0的方程,那么我们就知道第i个元是自由元。
5.通过这样的方法把矩阵消成上三角矩阵。


再看怎么得到对角矩阵:
1.我们有了上三角矩阵。
2.我们可以用第i行消去[1,i - 1]的系数。
3.得到对角矩阵。


这样,我们可以计算n阶行列式的值:
对于任意的上(下)三角行列式,我们有:

D=i=1naii

这样我们就可以在O(n3)的时间内解方程组,计算行列式的值。


于是开始刷第一只板子:
题目:bzoj1013

#include <cstdio>#include <cstring>#include <cmath>#include <algorithm>#define Rep(i,n) for(int i = 1;i <= n;i ++)#define RD(i,l,r) for(int i = l;i <= r;i ++)#define Dwn(i,n) for(int i = n;i;i --)using namespace std;const int N = 25;int n;const double eps = 1e-7;double pos[N][N],a[N][N],ans[N];double sqr(double x){return x * x;}int sgn(double x){return fabs(x) < eps ? 0 : x > 0 ? 1 : -1;}void Init(){    Rep(i,n + 1)Rep(j,n)a[i][j] = -2.0 * pos[i][j];    int p = n + 1;    Rep(i,n + 1)Rep(j,n)a[i][n + 1] -= sqr(pos[i][j]);    Rep(i,n)Rep(j,n + 1)a[i][j] -= a[p][j];}void Gauss(){    double x;int j;    Rep(i,n)    {        for(j = i;j <= n;j ++)if(sgn(a[j][i]))break;        if(j > n)continue;        if(i != j)Rep(k,n + 1)swap(a[i][k],a[j][k]);        RD(j,i + 1,n)        {            x = a[j][i] / a[i][i];            if(!sgn(x))continue;            RD(k,i,n + 1)a[j][k] -= x * a[i][k];        }    }    Dwn(i,n)    {        ans[i] = a[i][n + 1];        RD(j,i + 1,n)ans[i] -= ans[j] * a[i][j];        ans[i] /= a[i][i];    }}int main (){    scanf("%d",&n);    Rep(i,n + 1)Rep(j,n)scanf("%lf",&pos[i][j]);    Init();    Gauss();    Rep(i,n - 1)printf("%.3f ",ans[i]);printf("%.3f\n",ans[n]);    return 0;}

第二题:
1770: [Usaco2009 Nov]lights 燈
貝希和她的閨密們在她們的牛棚中玩遊戲。但是天不從人願,突然,牛棚的電源跳閘了,所有的燈都被關閉了。貝希是一個很膽小的女生,在伸手不見拇指的無盡的黑暗中,她感到驚恐,痛苦與絕望。她希望您能夠幫幫她,把所有的燈都給重新開起來!她才能繼續快樂地跟她的閨密們繼續玩遊戲! 牛棚中一共有N(1 <= N <= 35)盞燈,編號為1到N。這些燈被置於一個非常複雜的網絡之中。有M(1 <= M <= 595)條很神奇的無向邊,每條邊連接兩盞燈。 每盞燈上面都帶有一個開關。當按下某一盞燈的開關的時候,這盞燈本身,還有所有有邊連向這盞燈的燈的狀態都會被改變。狀態改變指的是:當一盞燈是開著的時候,這盞燈被關掉;當一盞燈是關著的時候,這盞燈被打開。 問最少要按下多少個開關,才能把所有的燈都給重新打開。 數據保證至少有一種按開關的方案,使得所有的燈都被重新打開。
首先每个灯只会要么被开一次要么不开。
所以我们求的是xi^xv0...k(i>kNULL)的异或方程组的解。
这样就可以高斯消元来求了。
然而还有自由元,我们暴力枚举所有自由元,加上剪枝即可。

#include <cstdio>#include <cstring>#include <algorithm>#include <cmath>#include <bitset>#define Rep(i,n) for(int i = 1;i <= n;i ++)#define Dwn(i,n) for(int i = n;i;i --)#define RD(i,l,r) for(int i = l;i <= r ;i ++)using namespace std;const int N = 55;int n,m,head[N],cnt,tot,mn,ans[N];bitset<N>a[N];void dfs(int x){    if(mn <= tot)return;    if(!x)return (void)(mn = tot);    if(a[x][x])    {        int t = a[x][n + 1];        RD(i,x + 1,n)if(a[x][i])t ^= ans[i];        ans[x] = t;        if(t)tot ++;        dfs(x - 1);        if(t)tot --;    }    else {ans[x] = 1;tot ++;dfs(x - 1),tot --;ans[x] = 0;dfs(x - 1);}}void Gauss(){    int j;    Rep(i,n)    {        for(j = i;j <= n;++ j)if(a[j][i])break;        if(j > n)continue;        if(i != j)swap(a[i],a[j]);        RD(k,i + 1,n)if(a[k][i])a[k] ^= a[i];    }}int main (){    memset(head,-1,sizeof(head));    scanf("%d%d",&n,&m);    mn = n + 1;    Rep(i,n)a[i][i] = 1,a[i][n + 1] = 1;    Rep(i,m){int b,bb;scanf("%d%d",&bb,&b);a[b][bb] = a[bb][b] = 1;}    Gauss();    dfs(n);    printf("%d\n",mn);    return 0;}

然后我们可以刷线性基的题了QAQ
然后来刷线性基的一个板子QAQ
这两个blog讲的很清楚(其中一个是来我们学校的大爷刚刚写的blog真是害怕QAQ

DMute的博客
线性基 复习总结
题目是hdu3949 XOR


我们来说什么是线性空间的基:
就是极大的能表示线性空间的线性无关组。
我们考虑一个n个方程的k个未知数的方程组:
如果这n个k维向量线性相关:
即:a⃗ V a⃗ =kb⃗  or a⃗ =b⃗ +c⃗ 
也就是说,存在至少一个”没有用”的向量。
比如(1,1) 和(2,2)这种东西就完全没有用。
我们可以知道,对于k维空间,用k个k维线性无关的向量组就可以表示。
我们可以把向量组写成矩阵的形式。
我们对矩阵进行高斯消元,消成对角阵,这样的话,我们就可以直接用这个矩阵里非0的向量来表示这个线性空间中的向量。
我们消元消出的这个极大的线性无关组就叫做线性基。


比如二维向量,它的其中一组线性基就是(1,0)(0,1)。
当然,(1,0)(0,7)这种也是。
所以,线性基有很多组。
但是线性基的大小是一定的。


线性基的大小等于描述这个线性空间的矩阵的秩。
行秩:横行的极大线性无关组的大小。
列秩:竖列的极大线性无关组的大小。
性质:矩阵的行秩等于列秩。
用矩阵乘法可以很愉悦地证明。


介绍一下极简线性基和上三角线性基。
1.极简线性基就趋近于对角矩阵的形式,即在消元时,假设有一个数的最高位的1在第k位,那么其他的数不存在最高位的1也在第k位。
2.上三角就是和高斯消元一样的消法。
好像大家还是比较喜欢消成第一种QAQ
不过注意一点,就是虽然说第一种消法在矩阵满秩的情况下是可以消掉的,但是矩阵欠秩的话就消不掉了。
比如[5,3]这种线性基,就一点办法都没有了QAQ


hdu XOR
求在N个数中最大的异或和。
考虑消出它的极简线性基,我们可以保证消元后得到的极简线性基异或起来之后得到最大值。

先贴个代码一会再说QAQ

#include <cstdio>#include <cstring>#include <algorithm>#include <cmath>#include <bitset>#define Rep(i,n) for(int i = 1;i <= n;i ++)#define RD(i,l,r) for(int i = l;i <= r ;i ++)using namespace std;typedef long long LL;const int N = 10005;int n;LL a[N];int Gauss(){    int cnp = 1;    for(int d = 63;d >= 0;d --)    {        RD(i,cnp,n)            if(a[i] >> d & 1)            {                swap(a[i],a[cnp]);                Rep(j,n)if(a[j] >> d & 1 && j != cnp)a[j] ^= a[cnp];                cnp ++;                break;            }    }    return -- cnp;}void print(int id){    puts("");    Rep(i,id)printf("@@@%lld ",a[i]);    puts("");}int main (){    int T,now,QQQ;LL K;bool p = 0;    for(now = 1,scanf("%d",&T);now <= T;now ++)    {        printf("Case #%d:\n",now);        scanf("%d",&n);        Rep(i,n)scanf("%lld",&a[i]);        int id = Gauss();        if(id != n)p = 1;        else p = 0;        scanf("%d",&QQQ);        while(QQQ --)        {            scanf("%lld",&K);            LL ans = 0;int LG = 0;            K -= p;if(!K){printf("0\n");continue;}            if((1ll << id) <= K){puts("-1");continue;}            for(int i = 0;K >> i;i ++)                if(K >> i & 1)ans ^= a[id - i];            printf("%lld\n",ans);        }    }    return 0;}

题目【BZOJ 2844】 albus就是要第一个出场
中文版题意:
在可重升序异或空间中寻找q的第一个出现的位置。
消元得出线性基,然后考虑每个线性基可不可选,显然如果这个线性基选了之后还小于q的话,那么就可选,并结算不可选的贡献。
对于消元后得到的0向量的个数k,我们把答案乘2k.

#include <cstring>#include <cstdio>#include <cmath>#include <algorithm>#define Rep(i,n) for(int i = 1;i <= n;i ++)#define Dwn(i,n) for(int i = n;i;i --) #define RD(i,l,r) for(int i = l;i <= r ;i ++)typedef long long LL;const int N = 100005;int n,a[N];using namespace std;int Gauss(){    int id = 1;    for(int d = 31;~ d;d --)    {        RD(i,id,n)            if(a[i] >> d & 1)            {                swap(a[i],a[id]);                Rep(j,n)if(a[j] >> d & 1 && j != id)a[j] ^= a[id];                id ++;break;            }    }    return -- id;}int main (){    scanf("%d",&n);     Rep(i,n)scanf("%d",&a[i]);    int id = Gauss(),q,now = 0;    scanf("%d",&q);    LL ans = 0,mod = 10086;    Rep(i,id)    {        if((a[i] ^ now) <= q) //如果异或上这个仍然小于q的话,那么这个基就可以选。并且加上这个基不选择的贡献。         {            now ^= a[i];            ans = (ans + (LL)(1ll << (id - i)) % mod) % mod;        }    }    Rep(i,n - id)ans <<= 1ll,ans %= mod; //有n个向量,id个线性基,那么则每个数字会出现2 ^ (n - id)次。     ans ++;    printf("%lld\n",ans);    return 0;}

然后我们就可以去水一水WC的题了QAQ
WC[2011]XOR
题意:无向图,寻找一条1到n的路径使得经过的所有的权值异或起来的结果最大。
高斯消元打错调了十五分钟。。。QAQ
题解:
0.我们任意选一条路线到达S - > T
1.对S - > T的任意路线,我们可以不走某一段而改走其简单环的补集。
2.这样我们可以求一个dfs树,对于非树边(u,v),构成了(u,v)的简单环。
3.容易知道,如果这个简单环不在S - >T的路径上,那么这个简单环就可以通过从某个点出发,再返回这个点的方法直接选取。
4.如果这个简单环在S - >T的选择的这条路径上,我们就需要在考虑选取这个简单环的同时,把之前选取的在简单环上的点那些边权都消掉。
5.异或满足消去律,我们相当于对于每一个环考虑选或者不选,再把这个环的权值和当前的权值异或起来,取最大的那个即可。
(实际上3和4没有本质区别)

#include <cstdio>#include <cstring>#include <algorithm>#define Rep(i,n) for(int i = 1;i <= n;i ++)#define Dwn(i,n) for(int i = n;~ i;i --)#define RD(i,l,r) for(int i = l;i <= r;i ++)#define v edge[i].to#define wl edge[i].w#define RepG(i,x) for(int i = head[x];~ i;i = edge[i].next)using namespace std;const int M = 200005;const int N = 50005;typedef long long LL;int n,m,cnt,head[N],tot;LL a[M + N],d[N];bool vis[N];struct Edge{int next,to;LL w;}edge[M];void save(int u,int b,LL c){    edge[cnt] = (Edge){head[u],b,c};head[u] = cnt ++;    edge[cnt] = (Edge){head[b],u,c};head[b] = cnt ++;}void dfs(int x){    vis[x] = 1;    RepG(i,x)        if(!vis[v])d[v] = d[x] ^ wl,dfs(v);        else a[++ tot] = d[v] ^ d[x] ^ wl;}int Gauss(){    int id = 1;    Dwn(dd,63)    {        RD(i,id,tot)            if((a[i] >> dd) & 1)            {                swap(a[i],a[id]);                Rep(j,tot)if(((a[j] >> dd) & 1) && j != id)a[j] ^= a[id];                ++ id;                break;            }    }    return -- id;}int main (){    memset(head,-1,sizeof(head));    scanf("%d%d",&n,&m);    int u,b;LL c;    Rep(i,m)scanf("%d%d%lld",&u,&b,&c),save(u,b,c);    dfs(1);    int id = Gauss();    LL ans = d[n];    Rep(i,id)ans = max(ans ^ a[i],ans);    printf("%lld\n",ans);    return 0;}

bzoj 3105 新nim游戏
肯定要给第二个人留一个线性基。
然后这样做的话,实际上就考虑如何使得第一个人的花费尽量少。
可以尽量让权值大的石子作为线性基来减少花费,这样就应该没什么问题了。(具体证明要用到我根本不会的拟阵233333)
排序之后,我们可以动态维护线性基,即每次向线性空间里加入一个向量,判断加入之后它是不是一个线性无关组。
如果无关它就是一个线性基。
然而实际上呢,不需要动态维护,我们直接排序好了之后进行高斯消元应该也是可以的。(待考证TAT)

#include <cstdio>#include <cstring>#include <algorithm>#define Rep(i,n) for(int i = 1;i <= n;i ++)#define Dwn(i,n) for(int i = n;~ i;i --)#define RD(i,l,r) for(int i = l;i <= r;i ++)#define v edge[i].tousing namespace std;const int N = 105;typedef long long LL;int n,a[N],w[N],b[N],vis[N];bool cmp(int a,int b){return a > b;}LL Gauss(){    sort(a + 1,a + 1 + n,cmp);    Rep(i,n)b[i] = a[i];//,printf("%d ",a[i]);puts("");    LL ans = 0;    Rep(p,n)    {        int now = a[p];        Dwn(d,31)            if(!w[d] && a[p] >> d & 1){w[d] = p;break;}            else if(a[p] >> d & 1)a[p] ^= a[w[d]];        if(a[p])ans += now,printf("%d\n",p); //如果是线性基的话,那么这一组就不用拿走     }/*  Dwn(d,31)    {        Rep(i,n)            if(!vis[i] && !w[d] && a[i] >> d & 1)            {                w[d] = i;vis[i] = 1;                Rep(j,n)                    if(a[j] >> d & 1 && j != i)                        a[j] ^= a[i];                break;            }    }    Rep(i,n)if(vis[i])ans += b[i];//,printf("%d ",i);*/    return ans;}int main (){    scanf("%d",&n);    LL ans = 0;    Rep(i,n)scanf("%d",&a[i]),ans += a[i];    LL res = Gauss();    if(!res)puts("-1");    else printf("%lld\n",ans - res);    return 0;   }

bzoj 4004: [JLOI2015]装备购买
垃圾题目卡精度。。。
然而自己还是变量没有赋初值23333
这个只是在一般线性空间里的线性基。和上题一样。

#include <cstdio>#include <cstring>#include <algorithm>#include <cmath>#define Rep(i,n) for(int i = 1;i <= n;i ++)#define RD(i,l,r) for(int i = l;i <= r ;i ++)using namespace std;struct vec{int c;long double a[1005];}v[1005];const double eps = 1e-5;int n,m,w[1005],ans,id;bool by[1005];int sgn(double x){fabs(x) < eps ? 0 : x > 0 ? 1 : -1;}bool cmp(vec a,vec b){return a.c < b.c;}bool zo(vec vv){Rep(i,m)if(sgn(vv.a[i]))return 1;return 0;}void Gauss(){    Rep(i,n)    {        Rep(j,m)            if(fabs(v[i].a[j]) > eps)            {                if(!w[j])                {                    w[j] = i;                    id ++,ans += v[i].c;                    break;                }                else if(w[j])                {                    long double x = v[i].a[j] / v[w[j]].a[j];                    RD(k,j,m)v[i].a[k] -= x * v[w[j]].a[k];                }            }    }}int main (){    scanf("%d%d",&n,&m);    Rep(i,n)Rep(j,m)scanf("%Lf",&v[i].a[j]);Rep(i,n)scanf("%d",&v[i].c);    sort(v + 1,v + 1 + n,cmp);    Gauss();    printf("%d %d\n",id,ans);    return 0;}
0 0