【XSY2524】唯一神 状压DP 矩阵快速幂 FFT

题目大意

  给你一个网格,每个格子有概率是\(1\)或是\(0\)。告诉你每个点是\(0\)的概率,求\(1\)的连通块个数\(\bmod d=0\)的概率。

  最开始所有格子的概率相等。有\(q\)次修改,每次修改一个格子的概率。要求输出初始时和每次修改后的概率。

  \(n\leq 200000,m\leq 3,d\leq 10,q\leq 1000\)

题解

  很容易想到状压DP:前\(i\)行在第\(i\)行的状态为\(j\)时连通块个数模\(d=k\)的概率。

  当\(m=3\)时每行状态有\(9\)个。

  然后用矩阵快速幂+线段树优化。

  显然会TLE。

  观察下状态转移方程:

\[f_{i+1,j',k'}+=f_{i,j,k}\times b
\]

  容易发现\(\Delta _k=k'-k\)与\(j\)无关。设

\[g_{i,j}=\sum_{k=0}^{d-1}f_{i,j,k}x^k
\]

  \(g\)的每次转移就等于乘上一个多项式。

  而且这是循环卷积(长度为\(d\))。

  我们可以在最开始先做一遍长度为\(d\)的DFT(直接暴力\(d^2\)),每次询问时把点值加起来做一遍IDFT,就可以得到多项式的常数项,也就是说我们想要的答案。

  时间复杂度:\(O(qs^3d\log n)\)

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<utility>
#include<cmath>
#include<functional>
#include<map>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
void sort(int &a,int &b)
{
if(a>b)
swap(a,b);
}
void open(const char *s)
{
#ifndef ONLINE_JUDGE
char str[100];
sprintf(str,"%s.in",s);
freopen(str,"r",stdin);
sprintf(str,"%s.out",s);
freopen(str,"w",stdout);
#endif
}
int rd()
{
int s=0,c;
while((c=getchar())<'0'||c>'9');
do
{
s=s*10+c-'0';
}
while((c=getchar())>='0'&&c<='9');
return s;
}
void put(int x)
{
if(!x)
{
putchar('0');
return;
}
static int c[20];
int t=0;
while(x)
{
c[++t]=x%10;
x/=10;
}
while(t)
putchar(c[t--]+'0');
}
int upmin(int &a,int b)
{
if(b<a)
{
a=b;
return 1;
}
return 0;
}
int upmax(int &a,int b)
{
if(b>a)
{
a=b;
return 1;
}
return 0;
}
const ll p=370137601;
const ll g=37;
ll fp(ll a,ll b)
{
ll s=1;
for(;b;b>>=1,a=a*a%p)
if(b&1)
s=s*a%p;
return s;
}
struct mat
{
ll a[10][10];
mat()
{
memset(a,0,sizeof a);
}
ll *operator [](int x)
{
return a[x];
}
};
int sz;
mat operator *(mat &a,mat &b)
{
mat c;
int i,j,k;
for(i=1;i<=sz;i++)
for(j=1;j<=sz;j++)
{
__int128 s=0;
for(k=1;k<=sz;k++)
s+=__int128(a[i][k])*b[k][j];
c[i][j]=s%p;
}
return c;
}
int level[200010];
mat w[11][20];
struct seg
{
struct node
{
mat *s;
int l,r,ls,rs;
};
node a[270000];
int cnt;
int rt;
void build(int &p,int l,int r,int k)
{
p=++cnt;
a[p].l=l;
a[p].r=r;
a[p].s=&w[k][level[r-l+1]];
if(l==r)
return;
int mid=(l+r)>>1;
build(a[p].ls,l,mid,k);
build(a[p].rs,mid+1,r,k);
}
void insert(int p,int x,mat &v)
{
if(a[p].l==a[p].r)
{
a[p].s=new mat;
*a[p].s=v;
return;
}
int mid=(a[p].l+a[p].r)>>1;
if(x<=mid)
insert(a[p].ls,x,v);
else
insert(a[p].rs,x,v);
a[p].s=new mat;
(*a[p].s)=(*a[a[p].ls].s)*(*a[a[p].rs].s);
}
};
seg a[11];
int n,m,k,q;
ll c[200010][4];
int dd[10][10];
int bb[10][10];
void get(ll *a,int x,int a1,int a2)
{
static ll b[10];
int i;
for(i=0;i<k;i++)
a[i]=0;
if(m==1)
{
b[1]=c[x][1];
b[2]=1-c[x][1];
a[dd[a1][a2]]=b[a2];
// for(l=0;l<k;l++)
// for(i=1;i<=sz;i++)
// for(j=1;j<=sz;j++)
// s[id(l,i)][id((l+dd[i][j])%k,j)]=b[j];
}
else if(m==2)
{
b[1]=c[x][1]*c[x][2]%p;
b[2]=(1-c[x][1])*c[x][2]%p;
b[3]=c[x][1]*(1-c[x][2])%p;
b[4]=(1-c[x][1])*(1-c[x][2])%p;
a[dd[a1][a2]]=b[a2];
// for(l=0;l<k;l++)
// for(i=1;i<=sz;i++)
// for(j=1;j<=sz;j++)
// s[id(l,i)][id((l+dd[i][j])%k,j)]=b[j];
}
else
{
b[1]=c[x][1]*c[x][2]%p*c[x][3]%p;
b[2]=(1-c[x][1])*c[x][2]%p*c[x][3]%p;
b[3]=c[x][1]*(1-c[x][2])%p*c[x][3]%p;
b[4]=c[x][1]*c[x][2]%p*(1-c[x][3])%p;
b[5]=(1-c[x][1])*(1-c[x][2])%p*c[x][3]%p;
b[6]=c[x][1]*(1-c[x][2])%p*(1-c[x][3])%p;
b[7]=(1-c[x][1])*(1-c[x][2])%p*(1-c[x][3])%p;
b[8]=b[9]=(1-c[x][1])*c[x][2]%p*(1-c[x][3])%p;
if(!bb[a1][a2])
a[dd[a1][a2]]=b[a2];
// for(l=0;l<k;l++)
// for(i=1;i<=sz;i++)
// for(j=1;j<=sz;j++)
// if(!bb[i][j])
// s[id(l,i)][id((l+dd[i][j])%k,j)]=b[j];
}
}
ll w1[200],w2[200];
void dft(ll *a)
{
static ll s[20];
int i,j;
for(i=0;i<k;i++)
{
s[i]=0;
for(j=0;j<k;j++)
s[i]=(s[i]+a[j]*w1[i*j])%p;
}
for(i=0;i<k;i++)
a[i]=s[i];
}
void idft(ll *a)
{
static ll s[20];
int i,j;
for(i=0;i<k;i++)
{
s[i]=0;
for(j=0;j<k;j++)
s[i]=(s[i]+a[j]*w2[i*j])%p;
}
ll inv=fp(k,p-2);
for(i=0;i<k;i++)
a[i]=s[i]*inv%p;
}
ll d[20];
void query()
{
int i,j;
for(j=1;j<=k;j++)
{
d[j-1]=0;
mat s=*a[j].a[a[j].rt].s;
for(i=1;i<=sz;i++)
d[j-1]=(d[j-1]+s[1][i])%p;
}
idft(d);
ll ans=d[0];
ans=(ans+p)%p;
printf("%lld\n",ans);
}
int main()
{
open("c");
scanf("%d%d%d%d",&n,&m,&k,&q);
int i,j,l,o;
if(m==1)
{
sz=2;
dd[1][2]=1;
}
else if(m==2)
{
sz=4;
dd[1][2]=dd[1][3]=dd[1][4]=1;
dd[2][3]=1;
dd[3][2]=1;
}
else
{
sz=9;
for(i=2;i<=7;i++)
dd[1][i]=1;
dd[1][8]=2;
dd[2][3]=dd[2][4]=dd[2][6]=dd[2][8]=1;
dd[3][2]=dd[3][4]=1;
dd[3][8]=2;
dd[4][2]=dd[4][3]=dd[4][5]=dd[4][8]=1;
dd[5][4]=dd[5][8]=1;
dd[6][2]=dd[6][8]=1;
dd[8][7]=k-1;
dd[8][3]=1;
dd[9][3]=1;
bb[1][9]=bb[2][9]=bb[3][9]=bb[4][9]=bb[5][9]=bb[6][9]=bb[8][9]=1;
bb[7][8]=bb[9][8]=1;
}
w1[0]=w2[0]=1;
w1[1]=fp(g,(p-1)/k);
w2[1]=fp(w1[1],p-2);
for(i=2;i<=k*k;i++)
{
w1[i]=w1[i-1]*w1[1]%p;
w2[i]=w2[i-1]*w2[1]%p;
}
int x,y;
ll p0,q0;
scanf("%lld%lld",&p0,&q0);
p0=p0*fp(q0,p-2)%p;
for(i=1;i<=n;i++)
for(j=1;j<=m;j++)
c[i][j]=p0;
for(i=1;i<=sz;i++)
for(j=1;j<=sz;j++)
{
get(d,1,i,j);
dft(d);
for(l=1;l<=k;l++)
w[l][0][i][j]=d[l-1];
}
for(l=1;l<=k;l++)
for(i=1;i<=18;i++)
w[l][i]=w[l][i-1]*w[l][i-1];
level[0]=-1;
for(i=1;i<=n;i++)
level[i]=level[i>>1]+1;
for(i=1;i<=k;i++)
a[i].build(a[i].rt,1,n,i);
query();
mat now[11];
for(i=1;i<=q;i++)
{
scanf("%d%d%lld%lld",&x,&y,&p0,&q0);
p0=p0*fp(q0,p-2)%p;
c[x][y]=p0;
for(j=1;j<=sz;j++)
for(l=1;l<=sz;l++)
{
get(d,x,j,l);
dft(d);
for(o=1;o<=k;o++)
now[o][j][l]=d[o-1];
}
for(j=1;j<=k;j++)
a[j].insert(a[j].rt,x,now[j]);
query();
}
// printf("sum=%d\n",sum);
return 0;
}
上一篇:SplitContainer如何实现左侧导航,正确显示和导航内容


下一篇:【BZOJ3590】[Snoi2013]Quare 状压DP