高斯消元&线性基
由于期望DP要用到高斯消元，这几天我学习了一下高斯消元，终于差不多学会了
高斯消元模板
void gauss(){
int now=1,to;double t;
for(int i=1;i<=n;i++){
/*for(to=now;!a[to][i]&&to<=n;to++);//做除法时减小误差，可不写
if(to!=now)for(int j=1;j<=n+1;j++)swap(a[to][j],a[now][j]);*/
t=a[now][i];
for(int j=1;j<=n+1;j++)a[now][j]/=t;
for(int j=1;j<=n;j++)
if(j!=now){
t=a[j][i];
for(int k=1;k<=n+1;k++)a[j][k]-=t*a[now][k];
}
now++;
}
}
BZOJ1013 高斯消元
#include<cstdio>
int n;double f[21],a[21][21];
double sqr(double x){return x*x;}
void gauss(){
int now=1,to;double t;
for(int i=1;i<=n;i++){
t=a[now][i];
for(int j=1;j<=n+1;j++)a[now][j]/=t;
for(int j=1;j<=n;j++)
if(j!=now){
t=a[j][i];
for(int k=1;k<=n+1;k++)
a[j][k]-=t*a[now][k];
}
now++;
}
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)scanf("%lf",&f[i]);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++){
double t;
scanf("%lf",&t);
a[i][j]=2*(t-f[j]);
a[i][n+1]+=sqr(t)-sqr(f[j]);
}
gauss();
for(int i=1;i<=n-1;i++)printf("%.3lf ",a[i][n+1]);
printf("%.3lf\n",a[n][n+1]);
}
利用高斯消元解期望DP：可阅读《浅析竞赛中一类数学期望问题的解决方法》
BZOJ3143：
每个点经过次数的期望为：f[i]=1+Σf[j]*(1/du[j])(i=1)
Σf[j]*(1/du[j])(i=1)（j为从i可达点，注意不包括n）
而每条边经过的期望次数为边连的两个点的期望次数除以各自度数的和
#include<cstdio>
#include<algorithm>
#define N 502
#define M 202222
int n,m,i,x,y,du[N],n1[M],n2[M];
double a[N][N],w[M],ans;
void gauss(){
int now=1,to;double t;
for(int i=1;i<=n;i++){
t=a[now][i];
for(int j=1;j<=n+1;j++)a[now][j]/=t;
for(int j=1;j<=n;j++)
if(j!=now){
t=a[j][i];
for(int k=1;k<=n+1;k++)a[j][k]-=t*a[now][k];
}
now++;
}
}
int main(){
scanf("%d%d",&n,&m);
for(i=1;i<=m;i++){
scanf("%d%d",&n1[i],&n2[i]);
du[n1[i]]++;du[n2[i]]++;
}
for(i=1;i<=m;i++){
a[n1[i]][n2[i]]-=1.0/du[n2[i]];
a[n2[i]][n1[i]]-=1.0/du[n1[i]];
}
for(i=1;i<n;i++)a[i][i]=1,a[i][n]=0;
a[1][n]=1;n--;
gauss();n++;
for(i=1;i<=m;i++)w[i]=a[n1[i]][n]/du[n1[i]]+a[n2[i]][n]/du[n2[i]];
std::sort(w+1,w+m+1);
for(i=1;i<=m;i++)ans+=(m-i+1)*w[i];
printf("%.3lf",ans);
}
BZOJ2337按位期望DP，若边为1则a[x][y]-=1/du[x],a[x][n+1]-=1/du[x]否则a[x][y]+=1/du[x]，然后a[i][i]-=1,a[n][n]=1,a[n][n+1]=0
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define N 103
#define eps 1e-8
using namespace std;
int n,m,i,j,k,u,l,now,x[N*N],y[N*N],z[N*N],du[N];
double ans,a[N][N];
int main(){
scanf("%d%d",&n,&m);
for(i=1;i<=m;i++){
scanf("%d%d%d",&x[i],&y[i],&z[i]);
du[x[i]]++;if(x[i]!=y[i])du[y[i]]++;
}
for(u=0;u<32;u++){
memset(a,0,sizeof(a));
for(i=1;i<=m;i++){
if(x[i]!=n){
if((1<<u)&z[i])a[x[i]][n+1]-=(double)1/du[x[i]],a[x[i]][y[i]]-=(double)1/du[x[i]];
else a[x[i]][y[i]]+=(double)1/du[x[i]];
}
if(y[i]!=n&&x[i]!=y[i]){
if((1<<u)&z[i])a[y[i]][n+1]-=(double)1/du[y[i]],a[y[i]][x[i]]-=(double)1/du[y[i]];
else a[y[i]][x[i]]+=(double)1/du[y[i]];
}
}
for(i=1;i<n;i++)a[i][i]-=1;a[n][n]=1;a[n][n+1]=0;
for(now=1,i=1;i<=n;i++,now++){
for(l=now;l<=n;l++)if(fabs(a[l][i])>eps)break;
if(l>n)continue;
double t=a[l][i];
for(j=1;j<=n+1;j++)a[l][j]/=t;
for(j=1;j<=n;j++)if(l!=j){
double t=a[j][i];
for(int k=1;k<=n+1;k++)a[j][k]-=t*a[l][k];
}
}
ans+=a[1][n+1]*(double)(1<<u);
}
printf("%.3lf",ans);
}
BZOJ2707 tarjan+高斯消元，按拓扑序倒序堆每个强联通分量做一遍高斯消元得到答案。
#include<bits/stdc++.h>
#define N 10010
#define M 1001010
#define CL(a) memset(a,0,sizeof(a))
using namespace std;vector<int>V[N];double a[203][203],f[N];
int n,m,i,S,T,x,y,t,id,tot,cnt,fir[N],rk[N],q[N],du[N],low[N],dfn[N],bl[N],ne[M],la[M],v[N];
void ins(int x,int y){la[++tot]=y;ne[tot]=fir[x];fir[x]=tot;du[x]++;}
void tj(int x){
low[x]=dfn[x]=++id;v[x]=1;q[++t]=x;int y,i,o=0;
for(i=fir[x];i;i=ne[i])if(!dfn[y=la[i]])tj(y),low[x]=min(low[x],low[y]);else if(v[y])low[x]=min(low[x],dfn[y]);
if(dfn[x]==low[x])for(cnt++;o!=x;rk[o]=V[cnt].size(),V[cnt].push_back(o),bl[o]=cnt,v[o]=0)o=q[t--];
}
void dfs(int x){
v[x]=2;if(x==bl[T]){v[x]=1;return;}
for(int i=0,y,j,z;i<V[x].size();i++)
for(j=fir[y=V[x][i]];j;j=ne[j])if(z=bl[la[j]],z!=x){
if(!v[z])dfs(z);if(v[z]==1)v[x]=1;
}
}
void get(int x){
int i,j,k,y,z,n=V[x].size();double t;
for(i=0;i<n;i++)for(j=fir[y=V[x][i]];j;j=ne[j])if(z=bl[la[j]],z!=x&&!v[z])get(z);
for(CL(a),i=0;i<n;i++){
a[i][i]=1;y=V[x][i];if(y==T)continue;a[i][V[x].size()]=1;
for(j=fir[y];j;j=ne[j])if(z=la[j],x==bl[z])a[i][rk[z]]-=1.0/du[y];else a[i][V[x].size()]+=1.0/du[y]*f[z];
}
for(i=0;i<n;i++){
for(t=a[i][i],j=0;j<=n;j++)a[i][j]/=t;
for(j=0;j<n;j++)if(i!=j)for(t=a[j][i],k=0;k<=n;k++)a[j][k]-=t*a[i][k];
}
for(i=0;i<n;i++)f[V[x][i]]=a[i][n];v[x]=1;
}
int main(){
for(scanf("%d%d%d%d",&n,&m,&S,&T);m--;ins(x,y))scanf("%d%d",&x,&y);
for(i=1;i<=n;i++)if(!dfn[i])tj(i);CL(v);dfs(bl[S]);
for(i=1;i<=cnt;i++)if(v[i]==2)return puts("INF"),0;
CL(v);get(bl[S]);printf("%.3lf",f[S]);
}
BZOJ3640 每一层转移相同，高斯消元求逆，时间复杂度O(n^3+n^2hp)。
#include<bits/stdc++.h>
#define N 155
using namespace std;double t,ans,a[N][N],b[N][N],f[10100][N],e[N];
int n,m,H,i,j,k,x,y,u[N][N],v[N],d[N];
int main(){
for(scanf("%d%d%d",&n,&m,&H),i=1;i<=n;i++)scanf("%d",&v[i]);
for(i=1;i<=m;i++){
scanf("%d%d",&x,&y);u[x][y]++;d[x]++;
if(x!=y)u[y][x]++,d[y]++;
}
for(i=1;i<n;i++)for(j=1;j<=n;j++)if(!v[j])a[j][i]-=1.0*u[i][j]/d[i];
for(i=1;i<=n;i++)a[i][i]+=1,b[i][i]=1;
for(i=1;i<=n;i++){
for(t=a[i][i],j=1;j<=n;j++)a[i][j]/=t,b[i][j]/=t;
for(j=1;j<=n;j++)if(i!=j)for(t=a[j][i],k=1;k<=n;k++)a[j][k]-=t*a[i][k],b[j][k]-=t*b[i][k];
}
for(f[H][1]=1,k=H;k;k--){
for(i=1;i<n;i++)for(j=1;j<=n;j++)if(v[j]&&k+v[j]<=H&&d[i])f[k][j]+=f[k+v[j]][i]*u[i][j]/d[i];
for(i=1;i<=n;i++)e[i]=f[k][i];
for(i=1;i<=n;i++)if(!v[i])for(f[k][i]=0,j=1;j<=n;j++)if(v[j]||j==1)f[k][i]+=b[i][j]*e[j];
ans+=f[k][n];
}
printf("%.8lf",ans);
}
线性基
这几天一页一页地做BZOJ水题，发现好几道500B题可以用线性基解决，就学习了一下
线性基就是求出a[i]之间异或的所有结果
for(int i=1;i<=m;i++)
for(int j=63;j>=0;j--)
if((a[i]>>j)&1){
if(!ins[j]){ins[j]=a[i];break;}
else a[i]^=ins[j];
}
BZOJ2215：求出所有环的值后直接线性基
#include<cstdio>
#define N 55555
#define M 222222
int n,m,x,y,tot,last[M],next[M],first[N];long long z,value[M],d[N],a[M],ins[77];bool vis[N];
void insert(int x,int y,long long z){last[++tot]=y;value[tot]=z;next[tot]=first[x];first[x]=tot;}
void dfs(int x){
vis[x]=1;
for(int i=first[x];i;i=next[i]){
int y=last[i];
if(!vis[y])d[y]=d[x]^value[i],dfs(y);
else a[++m]=d[y]^value[i]^d[x];
}
}
void Ins(){
for(int i=1;i<=m;i++)
for(int j=63;j>=0;j--)
if((a[i]>>j)&1){
if(!ins[j]){ins[j]=a[i];break;}
else a[i]^=ins[j];
}
long long ans=d[n];
for(int i=63;i>=0;i--)
if((ans^ins[i])>ans)ans^=ins[i];
printf("%lld",ans);
}
int main(){
scanf("%d%d",&n,&m);
while(m--){
scanf("%d%d%lld",&x,&y,&z);
insert(x,y,z);insert(y,x,z);
}
dfs(1);
Ins();
}
BZOJ3105：从大往小选加入线性基
#include<cstdio>
#include<algorithm>
using namespace std;
long long n,i,j,a[111],ins[31],ans,tot;
int main(){
scanf("%lld",&n);
for(i=1;i<=n;i++)scanf("%lld",&a[i]);
sort(a+1,a+n+1);
for(i=1;i<=n;i++)tot+=a[i];
for(i=n;i;i--){
int t=a[i];
for(j=30;j>=0;j--)if((a[i]>>j)&1){
if(!ins[j]){ins[j]=a[i];break;}
else a[i]^=ins[j];
}
if(a[i])ans+=t;
}
if(ans)printf("%lld\n",tot-ans);else puts("-1");
}
BZOJ2460：从大往小选加线性基
#include<cstdio>
#include<algorithm>
using namespace std;
int n,i,j;long long ans,ins[70];
struct EG{long long x,y;}eg[1111];
bool cmp(EG a,EG b){return a.y>b.y;}
int main(){
scanf("%d",&n);
for(i=1;i<=n;i++)scanf("%lld%lld",&eg[i].x,&eg[i].y);
sort(eg+1,eg+n+1,cmp);
for(i=1;i<=n;i++){
for(j=63;j>=0;j--)if((eg[i].x>>j)&1){
if(!ins[j]){ins[j]=eg[i].x;break;}
else eg[i].x^=ins[j];
}
if(eg[i].x)ans+=eg[i].y;
}
printf("%lld",ans);
}
BZOJ4004 从小到大选。
#include<cstdio>
#include<algorithm>
#include<cmath>
#define eps 1e-5
using namespace std;
int n,m,i,j,k,ans,num,ins[1111];
double rate;
struct EG{double x[1111];int c;}a[1111];
bool cmp(EG a,EG b){return a.c<b.c;}
int main(){
scanf("%d%d",&n,&m);
for(i=1;i<=n;i++)
for(j=1;j<=m;j++)
scanf("%lf",&a[i].x[j]);
for(i=1;i<=n;i++)scanf("%d",&a[i].c);
sort(a+1,a+n+1,cmp);
for(i=1;i<=n;i++){
for(j=1;j<=m;j++)if(fabs(a[i].x[j])>eps){
if(ins[j]){
rate=a[i].x[j]/a[ins[j]].x[j];
for(k=j;k<=m;k++)a[i].x[k]-=rate*a[ins[j]].x[k];
}
else{ins[j]=i,ans+=a[i].c,num++;break;}
}
}
printf("%d %d",num,ans);
}
BZOJ2844 每个元素出现次数为2^(n-总出现次数)
#include<cstdio>
#include<algorithm>
#define N 100010
#define p 10086
using namespace std;
int n,m,i,j,k,ans,sum,a[N],ins[32],crs[32];
int pow(int a,int b){for(sum=1;b;a=a*a%p,b>>=1)if(b&1)sum=sum*a%p;return sum;}
int main(){
scanf("%d",&n);
for(i=1;i<=n;i++)scanf("%d",&a[i]);
sort(a+1,a+n+1);
for(i=n;i;i--){
for(j=31;j>=0;j--)
if((a[i]>>j)&1){
if(!ins[j]){ins[j]=a[i];break;}
else a[i]^=ins[j];
}
}
for(i=0;i<=31;i++)if(ins[i])crs[i]=m++;
scanf("%d",&k);
for(i=0;i<=31;i++)if((k>>i)&1)if(ins[i])ans+=(1<<crs[i]);
ans%=p;printf("%d",(ans*pow(2,n-m)+1)%p);
}
BZOJ3237&3563&3569 三倍经验的神题，先构造出一颗树，然后给每个非树边赋一个随机数，对于每条树边，其值为和他相邻所有边地异或值，这样如果选出删去的边中任意一个子集异或和为0就是不合法的，这可以用线性基判
#include<bits/stdc++.h>
#define N 100100
#define M 1001000
using namespace std;bool f;
int n,m,i,j,x,y,k,mt,tot,u,w,Q,d[33],fir[N],ne[M],la[M],a[M],g[N],v[N];
int Rd(){return (32767*(rand()))+rand();}
void ins(int x,int y){la[++tot]=y;ne[tot]=fir[x];fir[x]=tot;}
void dfs(int x,int fa){
v[x]=++mt;int i,y;
for(i=fir[x];i;i=ne[i])if(la[i]!=fa)if(!v[y=la[i]])dfs(y,x),a[i/2]=g[y],g[x]^=g[y];else if(v[x]>v[la[i]])a[i/2]=Rd(),g[x]^=a[i/2],g[y]^=a[i/2];
}
int main(){
for(srand(20000108),scanf("%d%d",&n,&m),tot=1;m--;ins(x,y),ins(y,x))scanf("%d%d",&x,&y);
for(dfs(1,0),scanf("%d",&Q);Q--;u+=f,puts(f?"Connected":"Disconnected"),memset(d,0,sizeof(d)))
for(f=1,scanf("%d",&k),i=1;i<=k;i++){
for(scanf("%d",&x),w=a[x^u],j=31;j>=0;j--)if(w>>j&1)
if(!d[j]){d[j]=w;break;}else w^=d[j];
if(!w)f=0;
}
}
BZOJ3811 好题。首先如果一个在集合里的数能被其它数xor出来，那么这个数没有存在的意义，不用鸟它。找出所有有用的数就可以了，如果k≥3，这样的数不会太多，暴力就可以了。如果k≤2，再特判掉。我懒得写k≤2的了。。
#include<bits/stdc++.h>
using namespace std;typedef unsigned long long LL;
int n,k,K,i,j,l;LL p,q,S,z1,z2,o,sm,a[100100];
int main(){
for(scanf("%d%d",&n,&K),i=1;i<=n;i++)scanf("%llu",&a[i]);
for(i=100,j=1;i>=0;i--){
for(p=0,k=j;k<=n;k++)if(a[k]>>i&1){p=k;break;}
if(!p)continue;swap(a[p],a[j]);
for(k=j+1;k<=n;k++)if(a[k]>>i&1)a[k]^=a[j];j++;
}
for(S=1<<j,i=0;i<S;i++){
for(p=sm=k=0,q=1;k<j;k++)if(i>>k&1)sm^=a[k+1];
for(k=0;k<K;k++)p*=sm,q*=sm,p+=q>>j,q&=((1<<j)-1);
z1+=p;z2+=q;//z1+=z2/S;z2%=S;
}
z1+=z2/S;z2%=S;
printf("%llu",z1);
if(z2)for(printf(".");z2;z2%=S)z2*=10,printf("%llu",z2/S);
}
这几题似乎都要拟阵证明，但蒟蒻并不太会，感觉就是排下序一个一个线性基
2015年6月05日 11:06
高斯消元需要学的啊？
你小学数学没学过？
2015年6月18日 16:00
黈力您很强嘛
2015年6月18日 16:01
黈黈黈？黈黈黈！