一道比較繁瑣的DP題，細節比較多，寫了半小時調了一小時……

首先容易想到我們分開求騎士的位置和格子的通行情況，然後合併起來即可

前者比較容易，\(O(n^2\times t)\)的暴力DP很好想，然後在此基礎上把一行狀壓起來就可以做到\(O(n\times t)\)了

考慮如何處理某一時刻格子的通行情況，在涉及到倍數/因數的問題時，我們常常考慮按\(\sqrt n\)設閾值

首先當\(a_{i,j}\ge1000\)時，它存在的時刻數顯然\(\le1000\)，那麼可以直接暴力維護

考慮當\(a_{i,j}<1000\)時，我們可以考慮把在範圍內的質數全部搞出來，設計狀態\(g_{i,j,k}\)表示第\(k\)

然後顯然當\(t=p_1^{c_1}\times p_2^{c_2}\times\cdots\times p_m^{c_m}\)可以從所有的\(p_1^{c_1'}\times p_2^{c_2'}\times\cdots\times p_m^{c_m'}(c_1'\le c_1,c_2'\le c_2,\cdots,c_m'\le c_m)\)的狀態轉移來

然後我們給時刻分解質因數，將預處理好的\(g\)修改上去即可，然後我們注意到一段質數中間會出現一些次數為\(0\)的空隙，我們預處理下所有空隙的狀態即可

總體複雜度\(O(n\times t\log t)\)

#include<cstdio>
#include<algorithm>
#define RI register int
#define CI const int&
using namespace std;
const int N=35,M=500,S=1e6;
struct element
{
	int t,x,y;
	inline element(CI T=0,CI X=0,CI Y=0) { t=T; x=X; y=Y; }
	friend bool operator < (const element& A,const element& B)
	{
		return A.t<B.t;
	}
}a[N*N*1000],fc[M]; bool vis[S+5]; int n,t,sx,sy,tot,x,ansx[N*N],ansy[N*N],pt;
int f[2][N],g[500][N][N],sum[500][500][N],c[N],prime[S+5],cnt,num,d[S+5];
#define P(x) prime[x]
inline void init(CI n=S)
{
	RI i,j; for (i=2;i<=n;++i)
	{
		if (!vis[i]) P(++cnt)=i,d[i]=cnt;
		for (j=1;j<=cnt&&i*P(j)<=n;++j)
		{
			vis[i*P(j)]=1; d[i*P(j)]=j;
			if (i%P(j)==0) break;
		}
	}
	for (i=1;i<=cnt;++i) if (P(i)<=1000) num=i; else break;
}
inline void add1(int v,CI x,CI y)
{
	for (RI i=1;i<=num;++i)
	{
		int c=0; while (v%P(i)==0) ++c,v/=P(i);
		for (;c<=20;++c) g[i][c][x]|=(1<<y);
	}
}
inline void add2(CI v,CI x,CI y)
{
	for (RI i=1;i*v<=t;++i) a[++tot]=element(i*v,x,y);
}
inline void DP_sum(void)
{
	for (RI i=1,j,k,l;i<=num;++i) for (j=i;j<=num;++j) for (l=1;l<=n;++l)
	for (sum[i][j][l]=(1<<n+1)-1,k=i;k<=j;++k) sum[i][j][l]&=g[k][0][l];
}
inline void resolve(int x)
{
	RI i,j; for (pt=0;x!=1;)
	{
		int ct=0,t=d[x]; while (x%P(t)==0) x/=P(t),++ct;
		fc[++pt]=element(t,ct);
	}
	int lst=1; for (sort(fc+1,fc+pt+1),i=1;i<=pt&&fc[i].t<=num;++i)
	{
		for (j=1;j<=n;++j) c[j]&=g[fc[i].t][fc[i].x][j];
		if (lst<fc[i].t) for (j=1;j<=n;++j) c[j]&=sum[lst][fc[i].t-1][j]; lst=fc[i].t+1;
	}
	if (lst<num) for (j=1;j<=n;++j) c[j]&=sum[lst][num][j];
}
#undef P
int main()
{
	RI i,j,cur; for (scanf("%d%d%d%d",&n,&t,&sx,&sy),init(),i=1;i<=n;++i)
	for (j=1;j<=n;++j) if (scanf("%d",&x),x<1000) add1(x,i,j); else add2(x,i,j);
	/*for (cur=1;cur<=4;++cur) for (i=1;i<=n;++i) for (j=1;j<=n;++j)
	printf("%d ",g[cur][i][j]); putchar('\n');*/
	for (DP_sum(),sort(a+1,a+tot+1),i=cur=1,f[0][sx]|=(1<<sy);i<=t;++i)
	{
		int nw=i&1; for (j=1;j<=n;++j) f[nw][j]=0,c[j]=(1<<n+1)-1;
		for (j=1;j<=n;++j)
		{
			if (j>2) f[nw][j]|=(f[nw^1][j-2]>>1)|(f[nw^1][j-2]<<1);
			if (j+2<=n) f[nw][j]|=(f[nw^1][j+2]>>1)|(f[nw^1][j+2]<<1);
			if (j>1) f[nw][j]|=(f[nw^1][j-1]>>2)|(f[nw^1][j-1]<<2);
			if (j+1<=n) f[nw][j]|=(f[nw^1][j+1]>>2)|(f[nw^1][j+1]<<2);
		}
		/*RI k; for (printf("%d\n",i),j=1;j<=n;++j) for (k=1;k<=n;++k)
		printf("%d%c",(f[nw][j]>>k)&1," \n"[k==n]); putchar('\n');*/
		for (resolve(i);cur<=tot&&a[cur].t==i;++cur) c[a[cur].x]|=(1<<a[cur].y);
		/*for (j=1;j<=n;++j) for (k=1;k<=n;++k)
		printf("%d%c",(c[j]>>k)&1," \n"[k==n]); putchar('\n');*/
		for (j=1;j<=n;++j) f[nw][j]&=c[j];
		/*for (j=1;j<=n;++j) for (k=1;k<=n;++k)
		printf("%d%c",(f[nw][j]>>k)&1," \n"[k==n]); putchar('\n');*/
	}
	/*for (cur=1;cur<=4;++cur) for (i=1;i<=4;++i) for (j=1;j<=n;++j)
	printf("%d ",sum[cur][i][j]); putchar('\n');*/
	for (pt=0,i=1;i<=n;++i) for (j=1;j<=n;++j) if ((f[t&1][i]>>j)&1)
	ansx[++pt]=i,ansy[pt]=j; for (printf("%d\n",pt),i=1;i<=pt;++i)
	printf("%d %d\n",ansx[i],ansy[i]); return 0;
}