• 高斯消元总结

    A-Matrix Equation_第 45 届国际大学生程序设计竞赛（ICPC）亚洲区域赛（济南）（重现赛）

             自己写一个2维矩阵或者3维矩阵就可以发现对于每一列来说都是独立的，每一列的n个Cij都是都关系的，这就构成了一个n元一次方程组，其实这就是解一下这个方程组，但是他是问的有多少个矩阵，对于这个方程组构造出的矩阵来说，有多少自由元就说明有多少个Cij是可以任意取值的，也就是有1，2两种选择，不是自由元的就只能有1种；自由元是指系数都是0的未知数，自由元的个数就是高斯消元后0行的个数，所以用高斯消元就可以求出来；

            又因为有运算是模2的，所以也就可以看成是异或型矩阵，那么运算方式就要稍微的发生以下变化，比如

    上述图片来自2020icpc-济南站_琅歌的博客-CSDN博客 

    1. #include 
    2. using namespace std;
    3. #define endl '\n'
    4. #define int long long
    5. const int mod=998244353;
    6. const int inf=1e18;
    7. const int N = 2e5+100;
    8. const double eps=1e-8;
    9. int qpow(int a,int b)
    10. {
    11.     int res=1;
    12.     while(b)
    13.     {
    14.         if(b&1) res=res*a%mod;
    15.         a=a*a%mod;
    16.         b>>=1;
    17.     }
    18.     return res;
    19. }
    20. int getinv(int a){return qpow(a,mod-2LL);}
    21. int n;
    22. int a[220][220],b[220][220],d[220][220],fac[N];
    23. int gauss()
    24. {
    25.     int newline=1;
    26.     for(int k=1;k<=n;k++)
    27.     {
    28.         int p=newline;
    29.         for(int i=newline+1;i<=n;i++)
    30.         if(d[i][k]>d[p][k]) p=i;
    31.         if(d[p][k]==0) continue;
    32.         swap(d[newline],d[p]);
    33.         for(int i=newline+1;i<=n;i++)
    34.         {
    35.             if(d[i][k]==0) continue;
    36.             for(int j=1;j<=n;j++)
    37.             d[i][j]^=d[newline][j];
    38.         }
    39.         newline++;
    40.     }
    41.     newline--;
    42.     return n-newline;
    43. }
    44. void print()
    45. {
    46.     cout<<"--------------------\n";
    47.         for(int i=1;i<=n;i++)
    48.         {
    49.             for(int k=1;k<=n;k++)
    50.             cout<" ";
    51.             cout<
    52.         }
    53.         cout<<"--------------------\n";
    54. }
    55. signed main()
    56. {
    57.     //ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    58.     cin>>n;
    59.     fac[0]=1;
    60.     for(int i=1;i<=100000;i++) fac[i]=fac[i-1]*2LL%mod;
    61.     for(int i=1;i<=n;i++)
    62.     for(int j=1;j<=n;j++) cin>>a[i][j],d[i][j]=a[i][j];
    63.     for(int i=1;i<=n;i++)
    64.     for(int j=1;j<=n;j++) cin>>b[i][j];
    65.     int ans=1;
    66.     for(int j=1;j<=n;j++)
    67.     {
    68.         for(int i=1;i<=n;i++) 
    69.         for(int k=1;k<=n;k++) d[i][k]=a[i][k];
    70.         for(int i=1;i<=n;i++)
    71.         d[i][i]=a[i][i]^b[i][j];
    72.         //print();
    73.         int res=gauss();
    74.         //print();
    75.         //cout<
    76.         ans=ans*fac[res]%mod;
    77.         //cout<
    78.         //print();
    79.     }
    80.     cout<
    81.     system("pause");
    82.     return 0;
    83. }

    借助这个题又回顾了一下高斯消元的东西，高斯消元其实就是大一学的矩阵的行列变换使得矩阵变为行阶梯型，步骤为：

    1.设一变量newline代表进行到了第几行，最外层枚举列号k，内层枚举行i，找到最大的a[i][k]让其作为主元，将第i行与第newline行交换，如果最大的a[i][k]==0说明该变量是自由元，跳过去操作下一列；

    2.将主元下面的数都消掉，就是行变换

    3.操作结束后常数项/主元系数就是该未知数的解

    1. #include 
    2. using namespace std;
    3. #define endl '\n'
    4. #define int long long
    5. const int mod=998244353;
    6. const int inf=1e18;
    7. const int N = 2e5+100;
    8. const double eps=1e-8;
    9. int qpow(int a,int b)
    10. {
    11.     int res=1;
    12.     while(b)
    13.     {
    14.         if(b&1) res=res*a%mod;
    15.         a=a*a%mod;
    16.         b>>=1;
    17.     }
    18.     return res;
    19. }
    20. int sgn(double x)
    21. {
    22.     if(fabs(x)return 0;
    23.     else if(x<0) return -1;
    24.     else return 1;
    25. }
    26. int getinv(int a){return qpow(a,mod-2LL);}
    27. int n;
    28. double a[55][55];
    29. void gauss()
    30. {
    31.     int nl=1;
    32.     for(int k=1;k<=n;k++)
    33.     {
    34.         int p=nl;
    35.         for(int i=nl+1;i<=n;i++)
    36.         if(a[p][k]
    37.         //找一列中最大的数作为主元
    38.         if(sgn(a[p][k])==0) continue;
    39.         //如果最大的数是0那这个变量是一个自由元，跳过
    40.         swap(a[p],a[nl]);
    41.         //将主元所在的行换到上面
    42.         for(int i=1;i<=n;i++)
    43.         {
    44.             if(i==nl) continue;
    45.             double tmp=a[i][k]/a[nl][k];
    46.             for(int j=1;j<=n+1;j++)
    47.             a[i][j]-=tmp*a[nl][j];
    48.             //将该主元所在的列中的非零元素都消去，使得这一列中只有主元一个非零元素
    49.         }
    50.         nl++;
    51.     }
    52.     nl--;
    53.     if(nl
    54.     {
    55.         //如果行数不够n，说明含有自由元
    56.         for(int i=nl+1;i<=n;i++)
    57.         {
    58.             //如果常数项不为0说明是无解的
    59.             if(a[i][n+1]!=0){cout<<"-1\n";return;}
    60.         }
    61.         //如果常数项是0说明是无穷解
    62.         cout<<"0\n";return;
    63.     }
    64.     for(int i=1;i<=n;i++)
    65.     {
    66.         //常数项/主元系数就是未知数的解
    67.         cout<<"x"<"="<setprecision(2)<1]/a[i][i]<
    68.     }
    69. }
    70. signed main()
    71. {
    72.     ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    73.     cin>>n;
    74.     for(int i=1;i<=n;i++)
    75.     for(int j=1;j<=n+1;j++) cin>>a[i][j];
    76.     gauss();
    77.     system("pause");
    78.     return 0;
    79. }

    P2962 [USACO09NOV]Lights G - dfs+高斯消元

    这题也是异或方程组，可以发现假设每个灯按的次数是xi,那么对于每一个灯来说

    xj1^xj2^xj3^...^xi^...^xjm=1,j[1...m]是与i相连的灯，n个灯就构成了n元异或方程组，高斯消元后，如果有自由元的存在就要枚举这个灯开着优还是关掉优，最后取一个最优值即可

    题解 P2962 【[USACO09NOV]灯Lights】 - Youngsc_AFO 的博客 - 洛谷博客 (luogu.com.cn)

    1. #include 
    2. using namespace std;
    3. #define endl '\n'
    4. #define int long long
    5. const int mod=998244353;
    6. const int inf=1e18;
    7. const int N = 4e5+100;
    8. const double eps=1e-8;
    9. int qpow(int a,int b)
    10. {
    11.     int res=1;
    12.     while(b)
    13.     {
    14.         if(b&1) res=res*a%mod;
    15.         a=a*a%mod;
    16.         b>>=1;
    17.     }
    18.     return res;
    19. }
    20. int sgn(double x)
    21. {
    22.     if(fabs(x)return 0;
    23.     else if(x<0) return -1;
    24.     else return 1;
    25. }
    26. int getinv(int a){return qpow(a,mod-2LL);}
    27. int a[55][55],n,m;
    28. int vis[55],ans;
    29. void gauss()
    30. {
    31.     int nl=1;
    32.     for(int k=1;k<=n;k++)
    33.     {
    34.         int p=nl;
    35.         for(int i=nl+1;i<=n;i++) if(a[p][k]
    36.         if(a[p][k]==0) continue;
    37.         swap(a[p],a[nl]);
    38.         for(int i=1;i<=n;i++)
    39.         {
    40.             if(i==nl||a[i][k]==0) continue;
    41.             for(int j=1;j<=n+1;j++)
    42.             a[i][j]^=a[nl][j];
    43.         }
    44.         nl++;
    45.     }
    46.     nl--;
    47. }
    48. void dfs(int pos,int res)
    49. {
    50.     if(ans<=res) return;
    51.     if(pos==0)
    52.     {
    53.         ans=res;return;
    54.     }
    55.     int ma=n+1;
    56.     for(int i=pos;i<=n;i++)
    57.     {
    58.         if(a[pos][i]){ma=i;break;}
    59.     }
    60.     if(ma<=n)
    61.     {
    62.         int num=a[pos][n+1];
    63.         for(int i=ma+1;i<=n;i++)
    64.         {
    65.             if(a[pos][i]) num^=vis[i];
    66.         }
    67.         dfs(pos-1,res+num);
    68.     }
    69.     else
    70.     {
    71.         vis[pos]=0;
    72.         dfs(pos-1,res);
    73.         vis[pos]=1;
    74.         dfs(pos-1,res+1);
    75.         vis[pos]=0;
    76.     }
    77. }
    78. void print()
    79. {
    80.     cout<<"-----------------------------\n";
    81.     for(int i=1;i<=n;i++)
    82.     {
    83.         for(int j=1;j<=n+1;j++)
    84.         cout<
    85.         cout<
    86.     }
    87.     cout<<"------------------------------\n";
    88. }
    89. signed main()
    90. {
    91.     //ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    92.     //freopen("in.txt","r",stdin);
    93.         cin>>n>>m;
    94.         for(int i=1;i<=n;i++) a[i][i]=a[i][n+1]=1;
    95.         for(int i=1;i<=m;i++)
    96.         {
    97.             int u,v;cin>>u>>v;
    98.             a[u][v]=a[v][u]=1;
    99.         }
    100.         //print();
    101.         gauss();
    102.         ans=inf;
    103.         //print();
    104.         dfs(n,0);
    105.         cout<
    106.     system("pause");
    107.     return 0;
    108. }

    P5027 Barracuda - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

    根据题意列出方程然后枚举每一个方程是错误的，注意判断不合法的几个细节就可以，然后比较的时候一定要fabs，虽然都是正的，但是不加就是会wa,,,

    1. #include 
    2. using namespace std;
    3. #define endl '\n'
    4. #define int long long
    5. const int mod=998244353;
    6. const int inf=1e18;
    7. const int N = 4e5+100;
    8. const double eps=1e-10;
    9. int qpow(int a,int b)
    10. {
    11.     int res=1;
    12.     while(b)
    13.     {
    14.         if(b&1) res=res*a%mod;
    15.         a=a*a%mod;
    16.         b>>=1;
    17.     }
    18.     return res;
    19. }
    20. int sgn(double x)
    21. {
    22.     if(fabs(x)return 0;
    23.     else if(x<0) return -1;
    24.     else return 1;
    25. }
    26. bool integer(double s)
    27. {
    28.     if(abs(round(s)-s)return 1;
    29.     return 0;
    30. }
    31. int getinv(int a){return qpow(a,mod-2LL);}
    32. double a[105][105],b[105];
    33. int n;
    34. vector<int>v[105];
    35. int gauss()
    36. {
    37.     int nl=1;
    38.     for(int k=1;k<=n;k++)
    39.     {
    40.         int p=nl;
    41.         for(int i=nl+1;i<=n;i++) if(fabs(a[p][k])<fabs(a[i][k])) p=i;
    42.         //加上fabs才对，不知道为何，明明都是正的呀
    43.         if(sgn(a[p][k])==0) continue;
    44.         swap(a[p],a[nl]);
    45.         for(int i=1;i<=n;i++)
    46.         {
    47.             if(i==nl) continue;
    48.             double tmp=a[i][k]/a[nl][k];
    49.             for(int j=1;j<=n+1;j++)
    50.             a[i][j]-=tmp*a[nl][j];
    51.         }
    52.         nl++;
    53.     }
    54.     nl--;
    55.     if(nlreturn 0;
    56.     double res=0;
    57.     int ma=n+1;
    58.     for(int i=1;i<=n;i++)
    59.     {
    60.         double tmp=a[i][n+1]/a[i][i];
    61.         if(tmp<0||sgn(a[i][n+1])==0) return 0; 
    62.         int ztmp=tmp;
    63.         if(sgn(tmp-ztmp)!=0) return 0;
    64.         if(res
    65.         {
    66.             res=tmp;ma=i;
    67.         }
    68.     }
    69.     int cnt=0;
    70.     for(int i=1;i<=n;i++)
    71.     {
    72.         double tmp=a[i][n+1]/a[i][i];
    73.         if(sgn(res-tmp)==0) cnt++;
    74.     }
    75.     if(cnt>1) return 0;
    76.     return ma;
    77. }
    78. void print()
    79. {
    80.     cout<<"-----------------------------\n";
    81.     for(int i=1;i<=n;i++)
    82.     {
    83.         for(int j=1;j<=n+1;j++)
    84.         cout<
    85.         cout<
    86.     }
    87.     cout<<"------------------------------\n";
    88. }
    89. signed main()
    90. {
    91.     //ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    92.     //freopen("in.txt","r",stdin);
    93.     cin>>n;
    94.     for(int i=1;i<=n+1;i++)
    95.     {
    96.         int m;cin>>m;
    97.         for(int j=1;j<=m;j++)
    98.         {
    99.             int x;cin>>x;v[i].push_back(x);
    100.         }
    101.         cin>>b[i];
    102.     }
    103.     int ans=-1,mark=0;
    104.     for(int i=1;i<=n+1;i++)
    105.     {
    106.         for(int j=1;j<=n;j++)
    107.         for(int k=1;k<=n+1;k++) a[j][k]=0;
    108.         int cnt=0;
    109.         for(int j=1;j<=n+1;j++)
    110.         {
    111.             if(i==j)continue;
    112.             cnt++;
    113.             for(int k=0;ksize();k++)
    114.             {
    115.                 a[cnt][v[j][k]]=1;
    116.                 //cout<
    117.             }
    118.             a[cnt][n+1]=b[j];
    119.         }
    120.         //print();
    121.         int flag=gauss();
    122.         //print();
    123.         if(flag)
    124.         {
    125.             mark++;
    126.             ans=flag;
    127.             if(mark>1){ans=-1;break;}
    128.         }
    129.     }
    130.     if(ans==-1) cout<<"illegal\n";
    131.     else cout<
    132.     system("pause");
    133.     return 0;
    134. }

    P2447 [SDOI2010] 外星千足虫 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

    一开始就构造成m行n+1列的方程组直接高斯消元，记录最大用到的行就是最小的K，N太大就用bitset优化，有一篇题解是不需要bitset的，但是我不知道他的做法为什么对就没采用，还是老老实实的用bitset吧

    题解 P2447 【[SDOI2010]外星千足虫】 - YoungNeal 的博客 - 洛谷博客

    1. #include 
    2. using namespace std;
    3. #define endl '\n'
    4. #define int long long
    5. const int mod=998244353;
    6. const int inf=1e18;
    7. const int N = 4e5+100;
    8. const double eps=1e-10;
    9. int qpow(int a,int b)
    10. {
    11.     int res=1;
    12.     while(b)
    13.     {
    14.         if(b&1) res=res*a%mod;
    15.         a=a*a%mod;
    16.         b>>=1;
    17.     }
    18.     return res;
    19. }
    20. int sgn(double x)
    21. {
    22.     if(fabs(x)return 0;
    23.     else if(x<0) return -1;
    24.     else return 1;
    25. }
    26. bool integer(double s)
    27. {
    28.     if(abs(round(s)-s)return 1;
    29.     return 0;
    30. }
    31. int getinv(int a){return qpow(a,mod-2LL);}
    32. int n,m,b[2005],ans;
    33. char s[2005][1005];
    34. bitset<1005>a[2005];
    35. bool gauss()
    36. {
    37.     ans=0;
    38.     int nl=1;
    39.     for(int k=1;k<=n;k++)
    40.     {
    41.         int p=nl;
    42.         for(int i=nl;i<=m;i++) if(a[i][k]==1){p=i;break;}
    43.         if(a[p][k]==0) return 0;
    44.         ans=max(ans,p);
    45.         swap(a[p],a[nl]);
    46.         for(int i=1;i<=m;i++)
    47.         {
    48.  
    49.             if(i==nl||a[i][k]==0) continue;
    50.             // for(int j=k;j<=n+1;j++)
    51.             // a[i][j]^=a[nl][j];
    52.             a[i]^=a[nl];
    53.         }
    54.         nl++;
    55.     }
    56.     nl--;
    57.     return nl>=n;
    58. }
    59. void print()
    60. {
    61.     cout<<"-----------------------------\n";
    62.     for(int i=1;i<=n;i++)
    63.     {
    64.         for(int j=1;j<=n+1;j++)
    65.         cout<
    66.         cout<
    67.     }
    68.     cout<<"------------------------------\n";
    69. }
    70. signed main()
    71. {
    72.     ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    73.     //freopen("in.txt","r",stdin);
    74.     cin>>n>>m;
    75.     for(int i=1;i<=m;i++)
    76.     {
    77.         cin>>(s[i]+1);
    78.         for(int j=1;j<=n;j++) a[i][j]=(s[i][j]-'0');
    79.         int x;cin>>x;
    80.         a[i][n+1]=x;
    81.     }
    82.     if(gauss())
    83.     {
    84.         cout<
    85.         for(int i=1;i<=n;i++)
    86.         {
    87.             if(a[i][n+1]^a[i][i]) cout<<"Earth\n";
    88.             else cout<<"?y7M#\n";
    89.         }
    90.     }
    91.     else cout<<"Cannot Determine\n";
    92.     system("pause");
    93.     return 0;
    94. }

    P2973 [USACO10HOL]Driving Out the Piggies G - 高斯消元解方程组，概率

    假设期望是经过这个点几次才会爆炸，那么期望其实是和概率在数值上是一样的了，也就是把每走一步的代价看作1，也就相当于在求期望了，f[i]表示经过了多少次i点会爆炸，然后根据题意就可以列出

    f[i]=1+sum(f[j]*(1-P)*(1/in[j])) (i==1)

    f[i]=sum(f[j]*(1-P)*(1/in[j])) (i!=1)

    P是爆炸的概率，in[j]是j的度数，j可以到达i，j一共可以到达in[j]个点，j到达i的概率就是1/in[j],f[i]就可以通过j来转移，当i==1时因为一开始的起点就是1，所以1这个点一开始就经过了一次，所以要加上1

    1. #include 
    2. using namespace std;
    3. #define endl '\n'
    4. #define int long long
    5. const int mod=998244353;
    6. const int inf=1e18;
    7. const int N = 4e5+100;
    8. const double eps=1e-10;
    9. int qpow(int a,int b)
    10. {
    11.     int res=1;
    12.     while(b)
    13.     {
    14.         if(b&1) res=res*a%mod;
    15.         a=a*a%mod;
    16.         b>>=1;
    17.     }
    18.     return res;
    19. }
    20. int sgn(double x)
    21. {
    22.     if(fabs(x)return 0;
    23.     else if(x<0) return -1;
    24.     else return 1;
    25. }
    26. bool integer(double s)
    27. {
    28.     if(abs(round(s)-s)return 1;
    29.     return 0;
    30. }
    31. int getinv(int a){return qpow(a,mod-2LL);}
    32. int e[305][305],n,m,in[305];
    33. double p,q,a[305][305],P;
    34. void gauss()
    35. {
    36.     int nl=1;
    37.     for(int k=1;k<=n;k++)
    38.     {
    39.         int p=nl;
    40.         for(int i=nl+1;i<=n;i++) if(a[p][k]
    41.         if(sgn(a[p][k])==0) continue;
    42.         swap(a[p],a[nl]);
    43.         for(int i=1;i<=n;i++)
    44.         {
    45.             if(i==nl) continue;
    46.             double tmp=a[i][k]/a[nl][k];
    47.             for(int j=1;j<=n+1;j++)
    48.                 a[i][j]-=tmp*a[nl][j];
    49.         }
    50.         nl++;
    51.     }
    52.     for(int i=1;i<=n;i++)
    53.     {
    54.         double tmp=a[i][n+1]/a[i][i]*P;
    55.         cout<
    56.     }
    57. }
    58. signed main()
    59. {
    60.     //ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    61.     //freopen("in.txt","r",stdin);
    62.     cin>>n>>m>>p>>q;
    63.     P=p/q;
    64.     for(int i=1;i<=m;i++)
    65.     {
    66.         int u,v;
    67.         cin>>u>>v;
    68.         in[u]++;in[v]++;
    69.         e[u][v]=e[v][u]=1;
    70.     }
    71.     for(int i=1;i<=n;i++) a[i][i]=1;
    72.     for(int i=1;i<=n;i++)
    73.     {
    74.         for(int j=1;j<=n;j++)
    75.         {
    76.             if(e[i][j]) a[i][j]-=a[j][j]*(1.0-P)/in[j];
    77.         }
    78.     }
    79.     a[1][n+1]=1;
    80.     gauss();
    81.     system("pause");
    82.     return 0;
    83. }

  • 相关阅读:
    使用责任链模式实现登录风险控制
    python统计应用
    JAVA并发之谈谈你对AQS的理解
    大一新生常踩的5大坑，轻则损失钱财，重则影响毕业
    MATLAB——BP神经网络信号拟合程序
    API 与 SDK 之间的区别
    Naive UI 中使用message组件，在.vue文件之外使用
    数仓GreenPlum中数据实时同步的方式
    网络编程04-UDP的广播、组播
    TypeScript中使用superagent
  • 原文地址：https://blog.csdn.net/weixin_52621204/article/details/127707757