1007

题意：有一个 \(k\) 维空间，第 \(k\) 维的范围为 \([0,N_k]\)。一个人在 \((0,\cdots,0)\)，要走到 \((N_1,\cdots,N_k)\)。他的移动方式有 \(n\) 种，第 \(i\) 种用 \((y_{i,1},\cdots,y_{i,k})\) 表示，代表他在第 \(j\) 维上前进 \(y_{i,j}\) 个单位长。（当然他不能走出去）另外还给定了 \(m\) 个位置不能走。数据范围是 \(2\le k\le 16,1\le n \le 10^5,1\le m \le 1000\)。注意PDF的数据范围错了，应该是 \(2\le k \le16\)。

看到这个问题，我首先想到了一个类似的CF题。那个题是有一个巨大的棋盘，可以向上向右走，也是给定了 \(m\) 个位置不能走，问从左下角到右上角的方案数。这个题的解法是容斥：用 \(dp_i\) 表示到达第 \(i\) 个不能走的位置，此前不经过别的这种位置的方案数。那么 \(dp_i=ways((0,0),v_i)-\sum_{j=1}^{i-1}dp_j\cdot ways(v_j,v_i)\)。\(v_i\) 表示第 \(i\) 个不能走的点，\(ways(u,v)\) 表示从 \(u\) 走到 \(v\) 的方案数。把终点视作最后一个不能走的位置，跑这个 DP，就可以 \(O(m^2)\) 过这道题。这个题的题号是 CF559C。

回到 1007，我们发现这个做法可以照搬下来，那么就剩下求 \(ways\)。题目给了总点数 \(N=\prod_{i=1}^k(N_i+1)\le 2\cdot 10^5\)，稍微扩大一下，我们来求 \((0,...,0)\) 到各个点的方案数。

如果熟悉多项式方面的事情，应该可以想到用生成函数处理。定义

\[Y_i (x_1,\cdots,x_k)=x_1^{y_{i,1}}\cdots x_k^{y_{i,k}} \]

\[Y=\sum_{i=1}^nY_i \]

那么 \(Y^j\) 中项 \(x_1^{c_1}\cdots x_k^{c_k}\) 的系数就表示移动 \(j\) 次到达 \((c_1,c_2,\cdots,c_k)\) 的方法数。我们的目标函数就是 \(\sum_{i\ge 0}Y^i=\frac{1}{1-Y}\)，它的各项系数是我们要求的答案。

做一个 \(k\) 维的多项式求逆多少是有点逆天的，所以换一种描述方式。对所有点按照字典序重编号，那么每个 \(Y_i\) 其实都可以描述成一个下标的加法，只是有些位置因为越界了不能加。在多项式中对于越界，常常会选择“直接丢掉”的方式。但是如果只有这一个数的刻画也没法丢。

我对于这道题的思考到这里就进行不下去了。（所以我会认为前面都很简单doge

再来回顾一下刚才的定义：点 \((x_1,\cdots,x_k)\) 对应数 \(x_1+x_2n_1+\cdots+x_kn_1\cdots n_{k-1}\)。（记 \(n_j=N_j+1\)）这时直接对多项式求逆会出问题，也就是多项式乘法出了问题。

这时题解给出了一个相当妙的操作：增加一个变元 \(t\)，它在 \(x^i\) 项的次数是 \(\lambda(i)=\lfloor \frac{i}{n_1} \rfloor+\lfloor \frac{i}{n_1n_2} \rfloor+\cdots+\lfloor \frac{i}{n_1n_2\cdots n_{k-1}} \rfloor\)。在做多项式乘法时，如果 \(x^i\) 乘以 \(x^j\) 能够贡献给 \(x^{i+j}\)，就要求 \(\lambda(i)+\lambda(j)=\lambda(i+j)\)。观察到 \(\lfloor \frac{i}{n} \rfloor+\lfloor \frac{j}{n} \rfloor-\lfloor \frac{i+j}{n} \rfloor \in \{-1,0\}\)（高斯函数的常用性质），所以 \(\lambda(i)+\lambda(j)-\lambda(i+j) \in \{-k+1,\cdots,0\}\)。那么我们只要做 \(\bmod (t^k-1)\) 的多项式乘法，保留合法的项就可以了。

具体到实现细节上，假设要对这样的两个多项式 \(f,g\) 卷积，设 \(F_{i,j}=[\lambda(j)==i]f_jx^j\)，\(G\) 类似。然后对各个 \(F_{i},G_{i}\) DFT 转成点值。要在 \(t\) 这一维上处理这些点值，因为 \(k\) 很小，我们直接对 \(t\) 这个变元暴力卷积，也就是 \(F_{i_1,j}\cdot G_{i_2,j}\)（点值）累加到 \(F_{(i_1+i_2)\bmod k,j}\)。再把每个 \(F_i\) IDFT。最后只保留 \(f_i=F_{\lambda(i),i}\) 作为卷积的结果。

这样我们就可以对这个多项式求逆了！这样我们就彻底解决这个题了！

讲题后UPD：增加一个数学意义上的严谨解释。

我们重载了多项式乘法：\(f(x)=\sum_{i=0}^{n-1} a_ix^i\) 与 \(g(x)=\sum_{j=0}^{m-1} b_jx^j\) 的乘积是 \(f(x)g(x)=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1} a_ib_jx^{i+j}[\lambda(i)+\lambda(j)==\lambda(i+j)]\)。这里对 \(\lambda\) 函数的所有运算都在 \(\bmod k\) 意义下。

为此定义

\[f(x,t)=\sum_{i=0}^{n-1} a_ix_it^{\lambda(i)},g(x,t)=\sum_{j=0}^{m-1} b_jx_jt^{\lambda(j)} \]

再定义

\[F_p(x)=\sum_{i=0}^{n-1}a_ix^i[\lambda(i)==p],G_q(x)=\sum_{j=0}^{n-1}b_jx^j[\lambda(j)==q] \]

则

\[f(x,t)=\sum_{p=0}^{k-1}F_p(x)t^p,g(x,t)=\sum_{q=0}^{k-1}G_q(x)t^q \]

从而

\[\begin{aligned} f(x,t)g(x,t)&=\sum_{p=0}^{k-1}\sum_{q=0}^{k-1}F_p(x)G_q(x)t^{p+q \bmod k}\\&=\sum_{r=0}^{k-1}\sum_{p=0}^{k-1}\sum_{q=0}^{k-1}F_p(x)G_q(x)[p+q\equiv r \pmod k]t^r \end{aligned}\]

我们通过 DFT 求出 \(F_p,G_q\) 在各个原根位置的点值，要求 \(\sum_{p=0}^{k-1}\sum_{q=0}^{k-1}F_p(x)G_q(x)[p+q\equiv r \pmod k]\) 在各个原根位置的点值，那太容易了，暴力就行。
卷完之后只保留需要的项的系数就行。

#include<bits/stdc++.h>
using namespace std;
const int K=22,N=1e6+5,P=998244353,_G=3,invG=(P+1)/3;
int T,k,b[K],facb[K],n,m,lm[N];
int Y[N],ways[N],tr[N];
int qpow(int a,int b=P-2){int c=1;for(;b;b>>=1,a=1ll*a*a%P)if(b&1)c=1ll*c*a%P;return c;}
void tpre(int n){
    for(int i=0;i<n;i++)
        tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *f,int flag,int n){
    for(int i=0;i<n;i++)if(i<tr[i])swap(f[i],f[tr[i]]);
    for(int p=2;p<=n;p<<=1){
        int len=p>>1;
		int tG=qpow(flag==1?_G:invG,(P-1)/p);
        for(int k=0;k<n;k+=p){
            int buf=1;
            for(int l=k;l<k+len;l++){
                int tmp=1ll*buf*f[len+l]%P;
                f[len+l]=(f[l]-tmp+P)%P;f[l]=(f[l]+tmp)%P;
                buf=1ll*buf*tG%P;
            }
        }
    }
}
int F[K][N],G[K][N],H[N],res[K];
void mul(int*A,int*B,int n,int type,int *ans){
    int len=1;for(;len<=n*1.5;len<<=1);tpre(len);
    for(int i=0;i<k;i++)
		for(int j=0;j<len;j++){F[i][j]=0;if(type)G[i][j]=0;}
    for(int i=0;i<n;i++){F[lm[i]][i]=A[i];if(type)G[lm[i]][i]=B[i];}
    for(int i=0;i<k;i++){NTT(F[i],1,len);if(type)NTT(G[i],1,len);}
    for(int i=0;i<len;i++){
        for(int j=0;j<k;j++)res[j]=0;
        for(int j=0;j<k;j++)
            for(int l=0;l<k;l++){
                int c=(j+l)%k;
                res[c]=(res[c]+1ll*F[j][i]*G[l][i]%P)%P;
            }
        for(int j=0;j<k;j++)F[j][i]=res[j];
    }
    for(int i=0;i<k;i++)NTT(F[i],-1,len);
    int invn=qpow(len);
    for(int i=0;i<n;i++)ans[i]=1ll*F[lm[i]][i]*invn%P;
}
void inv(int*a,int n,int*ans){
    if(n==1){ans[0]=qpow(a[0]);return;}
    inv(a,(n+1)>>1,ans);mul(a,ans,n,1,H);mul(H,ans,n,0,H);
    for(int i=0;i<n;i++)ans[i]=(2ll*ans[i]-H[i]+P)%P;
}
int p[1005],dp[1005];
bool cmp(int u,int v){
    for(int i=1;i<=k;i++)
        if(u/facb[i-1]%b[i]>v/facb[i-1]%b[i])return false;
    return true;
}
int main(){
    facb[0]=1;
    scanf("%d",&T);
    while(T--){
        scanf("%d",&k);
        for(int i=1;i<=k;i++)
            scanf("%d",b+i),++b[i],facb[i]=facb[i-1]*b[i];
        scanf("%d%d",&n,&m);
        while(n--){
            int tmp=0;
            for(int i=1,x;i<=k;i++)
                scanf("%d",&x),tmp+=x*facb[i-1];
            --Y[tmp];
        }
        Y[0]=1-P;
        for(int i=0;i<facb[k];i++)Y[i]=(Y[i]+P)%P;
        memset(lm,0,sizeof(lm));
        for(int i=0;i<facb[k];i++){
            for(int j=1;j<k;j++)lm[i]+=i/facb[j];
            lm[i]%=k;
        }
        memset(ways,0,sizeof(ways));
        inv(Y,facb[k],ways);
        memset(Y,0,sizeof(Y)); 
        memset(p,0,sizeof(p));
        for(int i=1;i<=m;i++)
            for(int j=1,x;j<=k;j++)
                scanf("%d",&x),p[i]+=x*facb[j-1];
        sort(p+1,p+m+1);
        p[++m]=facb[k]-1;
        memset(dp,0,sizeof(dp));
        for(int i=1;i<=m;i++){
            dp[i]=ways[p[i]];
            for(int j=1;j<i;j++)if(cmp(p[j],p[i]))
                dp[i]=(dp[i]-1ll*ways[p[i]-p[j]]*dp[j]%P+P)%P;
        }
        printf("%d\n",dp[m]);
    }
    return 0;
}