一道比较繁琐的DP题，细节比较多，写了半小时调了一小时……
首先容易想到我们分开求骑士的位置和格子的通行情况，然后合并起来即可
前者比较容易，(O(n^2 imes t))的暴力DP很好想，然后在此基础上把一行状压起来就可以做到(O(n imes t))了
考虑如何处理某一时刻格子的通行情况，在涉及到倍数/因数的问题时，我们常常考虑按(sqrt n)设阈值
首先当(a_{i,j}ge1000)时，它存在的时刻数显然(le1000)，那么可以直接暴力维护
考虑当(a_{i,j}<1000)时，我们可以考虑把在范围内的质数全部搞出来，设计状态(g_{i,j,k})表示第(k)行整除(p_i^j)的关系（其中(p_i)表示第(i)个质数）
然后显然当(t=p_1^{c_1} imes p_2^{c_2} imescdots imes p_m^{c_m})可以从所有的(p_1^{c_1'} imes p_2^{c_2'} imescdots imes p_m^{c_m'}(c_1'le c_1,c_2'le c_2,cdots,c_m'le c_m))的状态转移来
然后我们给时刻分解质因数，将预处理好的(g)修改上去即可，然后我们注意到一段质数中间会出现一些次数为(0)的空隙，我们预处理下所有空隙的状态即可
总体复杂度(O(n imes tlog 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('
');*/
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
",i),j=1;j<=n;++j) for (k=1;k<=n;++k)
printf("%d%c",(f[nw][j]>>k)&1,"
"[k==n]); putchar('
');*/
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,"
"[k==n]); putchar('
');*/
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,"
"[k==n]); putchar('
');*/
}
/*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('
');*/
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
",pt),i=1;i<=pt;++i)
printf("%d %d
",ansx[i],ansy[i]); return 0;
}