【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;
}