直接樸素模擬一下，時間複雜度是 \(\mathcal{O}(TN^2)\) ，可以通過 \(40\%\) 個測試點。

我們只需要記錄 \(0/1\) 表示是否可達，直接狀壓，可以做到 \(\mathcal{O}(TN)\) 轉移。但是有些格子在某些時間不可到達，處理這個限制仍然是 \(\mathcal{O}(TN^2)\) 。這個做法可以優化約 \(10\) 倍常數，可以通過 \(70\%\) 的測試點。

然後筆者天真的認為卡常可以過，卡了一年還是沒有卡過去。

想了一下發現這道題非常的數論，所以還是考慮用數論的方法解決。

首選對於某個位置上的數 \(a_{i,j}\) ，我們可以把 \(a_{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;
}