这一整场比赛质量都很高，可以去看看，link
一个点 ((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}}) ，即每一次改变这一块的概率。
最容易想到的应该是算这个式子：
可惜我数学不好不会算。。。事实上这个是能算的，用二项式定理构造一玩意： (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) 的才过掉：
后来想到这玩意可以直接推通项，一个裸的一阶递推。
令 (A=1-2p,B=p) ，那么 (f_i=Af_{i-1}+B)
两边同时除掉 (A_i) ，然后搞一堆 (f_i-f_{i-1}) 的形式加起来（差分）发现是一个等比数列。推到最后应该是
但是注意边界是 (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;
}