zoukankan      html  css  js  c++  java
  • 【XSY2524】唯一神 状压DP 矩阵快速幂 FFT

    题目大意

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

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

      (nleq 200000,mleq 3,dleq 10,qleq 1000)

    题解

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

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

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

      显然会TLE。

      观察下状态转移方程:

    [f_{i+1,j',k'}+=f_{i,j,k} imes 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^3dlog 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
    ",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
    ",sum);
    	return 0;
    }
    
  • 相关阅读:
    决战72hours
    学习中的十七条建议
    数学建模终结篇
    数学建模(7)建模开始
    ASP升级程序
    为blog挑选logo
    Mysql源代码分析系列(4): 主要调用流程(续)转载
    AS学习步骤
    什么是敏捷软件测试[转]
    Mysql源代码分析(6): Plugin架构介绍(续)转载
  • 原文地址:https://www.cnblogs.com/ywwyww/p/8513354.html
Copyright © 2011-2022 走看看