高斯消元
高斯消元,方程组求解
高斯消元的题目一般有以下特点:
可以列出方程
一般没有阶段性和次序性(不同于dp或搜索)
由若干个状态确定一个已知状态
最基础的,我们用double类型进行gauss
但是这种gauss有一个小缺陷:精度损失(毕竟gauss里面会有大量的除法)
以前有一道题就被卡精度,所以需要把系数化为一放在外面
做法总结如下:
首先我们枚举未知量ii,同时记录一个nownow,表示当前已经使用过(已经用来消元)的方程数
在剩下的方程中找到一个未知量ii的系数不为0的方程
为了便于之后的消元,我们把ta放到nownow位置上
之后再其余所有的方程中找到未知量ii的系数不为0的方程jj进行消元
a[j][k]=a[j][k]−a[j][i]/a[now][i]a[j][k]=a[j][k]−a[j][i]/a[now][i]
a[j][i]/a[now][i]a[j][i]/a[now][i]就是要把方程jj的未知量ii的系数消掉需要的倍数
void gauss(int n) {
int now=1,to; //now 处理到第几个方程
for (int i=1;i<=n;i++) //枚举未知量
{
for (to=now;to<=n;to++)
if (fabs(a[to][i])>eps)
break;
if (to>n) continue;
if (to!=now)
for (int j=1;j<=n+1;j++)
swap(a[to][j],a[now][j]);
for (int j=1;j<=n;j++) //回代
if (j!=now) {
double t=a[j][i]/a[now][i];
for (int k=1;k<=n+1;k++)
a[j][k]-=t*a[now][k];
}
now++; //处理的方程数递增
}
for (int i=1;i<=n;i++) //系数化为一
a[i][n+1]/=a[i][i];
}
这种gauss的一种经典应用就是求解概率
在我们的印象中,概率都是用dp解决
但是dp的前提就是没有后效性,也就是说只能在DAG(或树)之类的结构上进行
如果是在有环的图上,一个点的期望受各种制约,就不方便直接dp了
这种情况下,我们可以通过期望的计算式得到若干个方程,用高斯消元求解即可
经典例题:高斯消元+概率(一)
经典例题:高斯消元+概率(二)
以上只是高斯消元的一个形式,有题目要求解异或方程
没有精度问题,没有回代求值,相比较而言比较简单
这里需要强调一点,if (to>n) continue;
,说明未知量ii是一个自由元
结束gauss之后,方程组会存在n−numn−num个自由元,如果自由元个数不为零,该方程组没有唯一解
void gauss(int n) {
int now=1,to;
for (int i=1;i<=n;i++) {
for (to=now;to<=n;to++)
if (a[to][i]) break;
if (to>n) continue;
if (to!=now)
for (int j=1;j<=n+1;j++)
swap(a[to][j],a[now][j]);
for (int j=1;j<=n;j++)
if (j!=now&&a[j][i]) //a[j][i] 只有i的系数不为0的方程需要XOR
for (int k=1;k<=n+1;k++)
a[j][k]^=a[now][k];
now++;
}
}
上面都是小case
我认为有点难度的是求解同余方程组
然而归根结底还是高斯消元的思路,只不过把除法变成了乘法逆元,多了专门的回代求值
我试了一下去掉专门回代的方法,答案还是正确的(本质没有区别)
但是这样时间就是铁铮铮的O(n3)O(n3),容易TLE(可能和取模有关)
在O(松)的启迪下,我们还是老老实实把回代写在外面吧
在这里说一下无解的判断:元素ii的值!=0!=0,但是系数=0=0
bool gauss(int n) {
int now=1,to;
for (int i=1;i<=n;i++) { //枚举未知量
for (to=now;to<=n;to++)
if (abs(a[to][i]))
break;
if (to>n) continue;
if (to!=now)
for (int j=1;j<=n+1;j++)
swap(a[to][j],a[now][j]);
for (int j=now+1;j<=n;j++)
if (a[j][i]) //只有i的系数不为0的方程需要消元
{
int t=((a[j][i]*inv(a[now][i])%p)%p+p)%p;
//除法变乘法逆元
for (int k=1;k<=n+1;k++) {
a[j][k]-=t*a[now][k];
a[j][k]=(a[j][k]%p+p)%p;
}
}
now++;
}
for (int i=n;i>=1;i--) {
if (a[i][i]==0&&a[i][n+1]!=0) return 0; //无解
else if (a[i][i]==0&&a[i][n+1]==0) continue; //自由元
ans[i]=((a[i][n+1]*inv(a[i][i])%p)%p+p)%p; //系数化为一
//回代
for (int j=i-1;j>=1;j--)
if (a[j][i]) {
a[j][n+1]-=a[j][i]*ans[i];
a[j][n+1]=(a[j][n+1]%p+p)%p;
}
}
return 1;
}
如果模数不是质数,我们就不能用费马小定理计算逆元了
这个时候就需要辗转相除的gauss
在矩阵树定理(辗转相除的gauss)的讲解中有较详细的介绍
矩阵树定理中的高斯消元,只消成了上三角形
ll gauss(int n) {
int now,to;
int opt=0;
for (int i=1;i<=n;i++) {
for (to=now;to<=n;to++)
if (a[to][i]) break;
if (to>n) return 0; //无解
if (to!=now) {
for (int j=1;j<=n;j++)
swap(a[to][j],a[now][j]);
opt^=1;
}
for (int j=i+1;j<=n;j++)
while (a[j][i]!=0) {
int x=a[j][i]/a[now][i];
if (x!=0) {
for (int k=i;k<=n;k++) {
a[j][k]-=(ll)x*a[now][k];
a[j][k]=(a[j][k]%p+p)%p;
}
}
else {
for (int k=i;k<=n;i++) //系数比较小
swap(a[now][k],a[j][k]);
opt^=1;
}
}
now++;
}
ll ans=1;
for (int i=1;i<=n;i++)
ans=(ans*a[i][i])%p;
if (opt) ans=mod-ans;
return ans;
}
线性积
线性积详解
在OI中,通常用线性积解决异或和的问题
线性积有两种写法(具体视情况而定)
而本质上都是高斯消元(这就是为什么我要把这两种算法放在一起)
值得注意的是,如果一个数成功插入了线性积中,要及时breakbreak
上三角型
代码简单,解决最大异或和等问题比较方便
void prepare() {
for (int i=1;i<=n;i++)
for (int j=63;j>=0;j--) //int 31
if (a[i]>>j&1) {
if (!b[j]) {
b[j]=a[i];
break; //成功插入线性积,返回
}
else a[i]^=b[j];
}
}
对角线型
得到的线性积比较规范,解决第k大等问题比较方便
void prepare() {
for (int i=1;i<=n;i++)
for (int j=63;j>=0;j--)
if (a[i]>>j&1) {
if (b[j]) a[i]^=b[j];
else {
b[j]=a[i];
for (int k=i-1;k>=0;k--)
if (b[k]&&(b[j]>>k&1)) b[j]^=b[k];
for (int k=j+1;k<=63;k++)
if (b[k]>>j&1) b[k]^=b[j];
break; //成功插入线性积,返回
}
}
}
线性积有一个经典应用:求第k大(以及逆操作)
第k大异或和:
解决这种问题,对角线型比较方便
求出线性积,然而对于i=0−63i=0−63的每一位,b[i]b[i]不一定有值
所以我们把b[i]!=0b[i]!=0的线性积紧凑复制到BB数组中
i 0 1 2 3 4 5 6
b 1 2 0 8 16 32 0
B 1 2 8 16 32
之后将k二进制差分,如果第xx为是1,那么+B[x]+B[x]
这里再强调一个小问题:异或和能不能等于0
设cntcnt是BB数组的大小,nn是元素个数
如果cnt=ncnt=n,说明每个元素至少给线性积贡献一个1,这种情况下异或和不能等于0
这样我们就可以解决这个问题了
不过我们这样得到的是不计重复的第k大异或和
如果我们允许重复异或和怎么办?
这就需要通过线性积的性质来理解了
如果一个元素能够插入到在线性基中,则说明ta一定能为线性基贡献一个1
如果一个元素不能插入线性基,则说明ta一定能用若干个线性基组合而成
那么n−cntn−cnt就是这些自由元素(可以直接视为0)的个数
那么每种线性积的重复次数即为2n−cnt2n−cnt
询问异或和排名:
逆运算,我们还是求出对角线型的线性积,得到B数组
对于询问的异或值XX
我们先计算小于等于X−1X−1的数字数量,加上1就得到了XX的第一个出现位置
从高位到低位枚举,如果(X−1) Xor B[i](X−1) Xor B[i],排名k+2ik+2i(0什么的特判一下就好了)
同样,这种方法得到的是不重复计数的排名
如果我们允许重复,重复次数即为2n−cnt2n−cnt
当然,图上的路径最大异或和也可以用线性积解决,只要记住:
任意一条1到nn的路径的异或和,都可以由任意一条11到nn路径的异或和与图中的一些环的异或和来组合得到
如果我们要合并两个线性积怎么办?
直接暴力合并
经典例题:LCA+暴力合并线性积
最后:
姜爷的线性积最强例题