P7245 灯光效果

这一整场比赛质量都很高，可以去看看，link

P7245 灯光效果

一个点 ((i,j)) 被改变权值当且仅当 (x_{i_1}le ile x_{i_2},y_{j_1}le jle y_{j_2})

当它总共被改变奇数次就会对答案产生 (1) 的贡献。

显然全局只有 (m^2) 种本质不同的点，即这 (2m) 个点把整个“舞台屏幕”分成了 (m^2) 个矩形，每一个矩形内的点对答案的贡献是相同的。

一个非常naive的想法就是枚举这些点被改变了几次权值。这个就是在右下角选一个点，左上角选一个点，这两个值对于一块区域可以 (O(1)) 计算。

对于 ((x_{i-1},y_{j-1}),(x_i,y_j)) 这两个点组成的矩形，右下角有 ((n-i+1)(n-j+1)) 个点，左上角有 ((i-1)(j-1)) 个点。

每一次合法的选取两个关键点有 (dfrac{inom{m^2}{2}}{2}) 种方法（就是选两个块，然后去掉左下和右上选点的情况，根据对称性刚好是一半）

令 (p=dfrac{(n-i+1)(n-j+1)(i-1)(j-1)}{frac{m^2(m-1)^2}{4}}) ，即每一次改变这一块的概率。

最容易想到的应该是算这个式子：

[sum_{i=1}^{k}[i\%2=1]inom{k}{i}p^i(1-p)^{k-i} ]

可惜我数学不好不会算。。。事实上这个是能算的，用二项式定理构造一玩意： (dfrac{((1-p)+p)^k-((1-p)-p)^k}{2}) 就好了。

我这种菜鸡会垃圾做法，常数被上面那个高明做法吊锤了好几倍。

令 (f_i) 表示进行了 (i) 次变换之后这块区域为 (1) 的概率。

显然有 (f_i=(1-p)f_{i-1}+p(1-f_{i-1}))

这个式子可以矩乘。

我一开始推了个垃圾 (3 imes 3) 的矩阵，被卡常了。

后来推成 (2 imes 2) 的才过掉：

[egin{bmatrix} f_i,f_{i-1} end{bmatrix} * egin{bmatrix} 2-2p,1\ 2p-1,0 end{bmatrix} =egin{bmatrix} f_{i+1},f_{i} end{bmatrix} ]

后来想到这玩意可以直接推通项，一个裸的一阶递推。

令 (A=1-2p,B=p) ，那么 (f_i=Af_{i-1}+B)

两边同时除掉 (A_i) ，然后搞一堆 (f_i-f_{i-1}) 的形式加起来（差分）发现是一个等比数列。推到最后应该是

[f_n=dfrac{B(1-(frac{1}{A})^{n+1})}{1-frac{1}{A}}cdot A^{n} ]

但是注意边界是 (f_1=p) ，所以 (f_{k-1}) 才是答案。

原本以为常数会小很多，结果常数大了一倍。。。

同为 (log) ，你怎么被一个带 (4) 倍常数的矩乘吊起来打？？！！！

Code（矩乘）：

#define mod 998244353
const int N=1005;
int n,k,m,x[N],y[N],ans,ivm;
namespace math{
inline int qpow(int n,int k){int res=1;for(;k;k>>=1,n=1ll*n*n%mod)if(k&1)res=1ll*n*res%mod;return res;}
inline void fmod(int&x){x-=mod,x+=x>>31&mod;}
}
using math::qpow;
using math::fmod;

struct Matrix{
	int a[2][2];
	Matrix(){memset(a,0,sizeof(a));}
	Matrix operator * (const Matrix&t){
		Matrix res;
		// rep(i,0,1)rep(j,0,1)rep(k,0,1)fmod(res.a[i][j]+=1ll*a[i][k]*t.a[k][j]%mod);
		fmod(res.a[0][0]+=1ll*a[0][0]*t.a[0][0]%mod),fmod(res.a[0][0]+=1ll*a[0][1]*t.a[1][0]%mod);
		fmod(res.a[0][1]+=1ll*a[0][0]*t.a[0][1]%mod),fmod(res.a[0][1]+=1ll*a[0][1]*t.a[1][1]%mod);
		fmod(res.a[1][0]+=1ll*a[1][0]*t.a[0][0]%mod),fmod(res.a[1][0]+=1ll*a[1][1]*t.a[1][0]%mod);
		fmod(res.a[1][1]+=1ll*a[1][0]*t.a[0][1]%mod),fmod(res.a[1][1]+=1ll*a[1][1]*t.a[1][1]%mod);
		return res;
	}
	void e(){a[0][0]=a[1][1]=1,a[1][0]=a[0][1]=0;}
};
Matrix Matrix_qpow(Matrix a,int k){
	Matrix res;res.e();
	for(;k;k>>=1,a=a*a)if(k&1)res=res*a;
	return res;
}

int calc(int yx,int zs,int k){
	int tmp=1ll*yx*zs%mod*ivm%mod;
	Matrix tur,sta;
	sta.a[0][0]=tmp;
	tur.a[0][0]=mod+2-tmp*2%mod,tur.a[0][1]=1;
	tur.a[1][0]=(tmp*2-1)%mod,tur.a[1][1]=0;
	sta=sta*Matrix_qpow(tur,k-1);
	return sta.a[0][0];
}
signed main(){
	n=read(),m=read(),k=read(),ivm=qpow(1ll*(m*m)*(m*m-m*2+1)/4%mod,mod-2);
	for(int i=1;i<=m;++i)x[i]=read();
	for(int i=1;i<=m;++i)y[i]=read();
	for(int i=1;i<=m;++i)
		for(int j=1;j<=m;++j){
			int num=1ll*(x[i]-x[i-1])*(y[j]-y[j-1])%mod;
			int tmp=calc((m-i+1)*(m-j+1),(i-1)*(j-1),k);
			fmod(ans+=1ll*num*tmp%mod);
		}
	cout<<ans<<'
';
	return 0;
}

Code（通项）

#define mod 998244353
const int N=1005;
int n,k,m,x[N],y[N],ans,ivm;
namespace math{
inline int qpow(int n,int k){int res=1;for(;k;k>>=1,n=1ll*n*n%mod)if(k&1)res=1ll*n*res%mod;return res;}
inline void fmod(int&x){x-=mod,x+=x>>31&mod;}
}
using math::qpow;
using math::fmod;


int calc(int yx,int zs,int k){
	int tmp=1ll*yx*zs%mod*ivm%mod;
	int A=mod+1-tmp*2%mod,B=tmp,iA=qpow(A,mod-2),res;
	res=1ll*B*(1-qpow(iA,k)+mod)%mod*qpow(mod+1-iA,mod-2)%mod*qpow(A,k-1)%mod;
	return res;
}
signed main(){
	n=read(),m=read(),k=read(),ivm=qpow(1ll*(m*m)*(m*m-m*2+1)/4%mod,mod-2);
	for(int i=1;i<=m;++i)x[i]=read();
	for(int i=1;i<=m;++i)y[i]=read();
	for(int i=1;i<=m;++i)
		for(int j=1;j<=m;++j){
			int num=1ll*(x[i]-x[i-1])*(y[j]-y[j-1])%mod;
			if(num)fmod(ans+=1ll*num*calc((m-i+1)*(m-j+1),(i-1)*(j-1),k)%mod);
		}
	cout<<ans<<'
';
	return 0;
}