【题解】[COCI2011-2012#1] SKAKAC

直接朴素模拟一下，时间复杂度是 \(\mathcal{O}(TN^2)\) ，可以通过 \(40\%\) 个测试点。

我们只需要记录 \(0/1\) 表示是否可达，直接状压，可以做到 \(\mathcal{O}(TN)\) 转移。但是有些格子在某些时间不可到达，处理这个限制仍然是 \(\mathcal{O}(TN^2)\) 。这个做法可以优化约 \(10\) 倍常数，可以通过 \(70\%\) 的测试点。

然后笔者天真的认为卡常可以过，卡了一年还是没有卡过去。

想了一下发现这道题非常的数论，所以还是考虑用数论的方法解决。

首选对于某个位置上的数 \(a_{i,j}\) ，我们可以把 \(a_{i,j}\) 倍数的时间对应的 \((i,j)\) 标记为可达，其余标记为不可达，然后每个时间与是否可达的状态按位与一下即可。

时间复杂度 \(\mathcal{O}(TN\log T)\) ，貌似卡一下就能过，但非常不幸的是空间复杂度是 \(\mathcal{O}(TN)\)，而我们只有 \(64M\) 。

考虑根号分治，对于 \(a_{i,j}\ge \sqrt{T}\) ，最多有 \(\sqrt T\) 个时间点该格子是开放的，暴力维护即可，时间复杂度 \(\mathcal{O}(N^2\sqrt T)\) 。

对于 \(a_{i,j}< \sqrt T\) ，最多有 \(168\) 个质数。

我们记录 \(f[i][j][k]\) 表示第 \(i\) 个质数 \(p_i\)，时间为 \(p_i^j\) 的倍数，第 \(k\) 行的状态，对每个 \(a_{i,j}\) 分解质因数即可。

然后对于每个时间 \(t\) ，我们对 \(t\) 分解质因数，然后对于每个质数的 \(f[i][j]\) 按位与即可。注意这里 \(j=0\) 的质数 \(p_i\) 也要计算进去。

直接计算的时间复杂度是 \(\mathcal{O}(|P|TN)\) ，\(|P|=168\) ，比朴素模拟还要劣。

但是观察到对于 \(t\) 合并答案的时候，\(t\) 最多只有 \(7\) 个不同的质因数，其余的质因数指数都是 \(0\) 。

所以我们可以预处理 \(g[l][r][k]=f[l][0][k]\ \&\ f[l+1][0][k]\ \&\ \cdots\ \&\ f[r][0][k]\) 。这部分的时间复杂度是 \(\mathcal{O}(|P|^2N)\) 。

所以总的时间复杂度是 \(\mathcal{O}(N^2\sqrt{T}+|P|N\log T+|P|^2N+TN)\) 。

#include<bits/stdc++.h>
#define rep(i,a,b) for(register int i=a;i<=b;i++)
#define pre(i,a,b) for(register int i=a;i>=b;i--)
#define N 31
#define M 1005
#define W 170
#define K 1000005
using namespace std;
int n,T,a[N],b[N],u[N][N];
int pm[W],v[M],tot;
void init(){
	rep(i,2,1000){
		if(!v[i])pm[++tot]=i;
		rep(j,1,tot){
			if(i*pm[j]>1000)break;
			v[i*pm[j]]=1;
			if(i%pm[j]==0)break;
		}
	}
}
typedef pair<int,int> Pr;
vector<Pr>c[K];
int f[W][N][N],g[W][W][N];
int h[K],nxt[K+K+(K>>1)],val[K+K+(K>>1)],idx;
void calc(){
	rep(i,0,n-1)rep(j,0,n-1){
		if(u[i][j]>=1000){
			for(int k=u[i][j];k<=T;k+=u[i][j])c[k].push_back(make_pair(i,j));
		}
		else{
			rep(k,1,tot){
				int cnt = 0;
				while(u[i][j]%pm[k]==0)cnt++,u[i][j]/=pm[k];
				f[k][cnt][i] |= 1<<j;
			}
		}
	}
	rep(i,1,tot)rep(j,1,20)rep(k,0,n-1)f[i][j][k]|=f[i][j-1][k];
	rep(i,1,tot){
		rep(j,0,n-1)g[i][i][j]=f[i][0][j];
		rep(j,i+1,tot)rep(k,0,n-1)g[i][j][k]=g[i][j-1][k]&f[j][0][k];
	}
	pre(i,tot,1){
		for(int j=pm[i];j<=T;j+=pm[i]){
			++idx;val[idx]=i;nxt[idx]=h[j];h[j]=idx;
		}
	}
}
int p[N],q[N],pre;
void maintain(int t){
	memset(p,0,sizeof(p));memset(q,~0,sizeof(q));
	for(int i=0;i<(int)c[t].size();i++)p[c[t][i].first]|=1<<c[t][i].second;
	pre=1;int now=0;
	for(int i=h[t];i;i=nxt[i]){
		now++;assert(now<=20);
		if(val[i]>pre)rep(j,0,n-1)q[j]&=g[pre][val[i]-1][j];
		pre=val[i]+1;int y=t,cnt=0;
		while(y%pm[val[i]]==0)y/=pm[val[i]],cnt++;
		rep(j,0,n-1)q[j]&=f[val[i]][cnt][j];
	}
	if(pre<=tot)rep(i,0,n-1)q[i]&=g[pre][tot][i];
	rep(i,0,n-1)q[i]=(q[i]|p[i])&((1<<n)-1);
	rep(i,0,n-1)a[i]&=q[i];
}
int main(){
	init();scanf("%d%d",&n,&T);
	int sx,sy;scanf("%d%d",&sx,&sy);
	a[sx-1]=1<<(sy-1);
	rep(i,0,n-1)rep(j,0,n-1)scanf("%d",&u[i][j]);
	calc();
	rep(t,1,T){
		rep(i,2,n-1)b[i]|=(a[i-2]<<1)|(a[i-2]>>1);
		rep(i,0,n-3)b[i]|=(a[i+2]<<1)|(a[i+2]>>1);
		rep(i,1,n-1)b[i]|=(a[i-1]<<2)|(a[i-1]>>2);
		rep(i,0,n-2)b[i]|=(a[i+1]<<2)|(a[i+1]>>2);
		rep(i,0,n-1)a[i]=b[i],b[i]=0;
		maintain(t);
	}
	int sum=0;
	rep(i,0,n-1)rep(j,0,n-1)if((a[i]>>j)&1)sum++;
	printf("%d\n",sum);
	rep(i,0,n-1)rep(j,0,n-1)if((a[i]>>j)&1)printf("%d %d\n",i+1,j+1);
	return 0;
}