一个 \(B\) 进制 \(N\) 位数，每一位都非零。如果他能被自己的每个数位整除，则是一个“Digidivisible Numbers”。

例如，10进制的728，能够同时被7,2,8整除，所以它是一个Digidivisible Numbers。
又例如，8进制的264（等于10进制的180），能够同时被2,6,4整除，所以它也是一个Digidivisible Numbers。

给定 \(B,N\) ，即只考虑 \(B\) 进制 \(N\) 位数，进行 \(T\) 次询问，每次规定 \(1\sim B-1\) 中哪些数字可以使用，求这样的Digidivisible Numbers有多少个，答案对 \(999,999,001\) （是个质数）取模。

第一行两个整数 \(B,N\) 。

第二行一个整数 \(T\) 。

接下来 \(T\) 行，每行一个长度为 \(B\) 的01字符串 \(s_0\cdots s_{B-1}\) ，如果 \(s_i=1\) 则表示数字 \(i\) 可以使用。保证 \(s_0=0\) 且至少有一个数字可以使用。

每次询问输出一行一个整数答案。

样例解释

第一次询问，10进制的3位数Digidivisible Numbers共有56个。
第二次询问，每个数位都必须是偶数，共有17个：222, 224, 244, 248, 264, 288, 424, 444, 448, 488, 624, 648, 666, 824, 848, 864, 888。

数据范围

\(2\leq B\leq 10\)
\(1\leq N\leq 10^9\)
\(1\leq T\leq 2^{B-1}-1\)

数据由nodgd提供，原题链接

题意

一个 \(B\) 进制 \(N\) 位数，每一位都非零。如果他能被自己的每个数位整除，则是一个“Digidivisible Numbers”。

例如，10进制的728，能够同时被7,2,8整除，所以它是一个Digidivisible Numbers。
又例如，8进制的264（等于10进制的180），能够同时被2,6,4整除，所以它也是一个Digidivisible Numbers。

给定 \(B,N\) ，即只考虑 \(B\) 进制 \(N\) 位数，进行 \(T\) 次询问，每次规定 \(1\sim B-1\) 中哪些数字可以使用，求这样的Digidivisible Numbers有多少个，答案对 \(999,999,001\) （是个质数）取模。

数据范围： \(2\leq B\leq 10\) ， \(1\leq N\leq 10^9\) ， \(1\leq T\leq 2^{B-1}-1\) 。

时间限制： \(2\) 秒。

算法框架

设 \(M=\mathrm{lcm}(1,2,\cdots,B-1)\) ，最大值为 \(2520\) 。

设 \(P=999,999,001\) 。注意到 \(P\) 虽然不是NTT模数，但 \(2520\;\Big|\;(P-1)\) ，这个性质后面会用到。

第一步，计算 \(f[S][i]\) 表示只能用 \(S\) 中的数字且 \(\bmod P=i\) 的 \(N\) 位数有多少个，其中 \(S\subseteq\{1,2,\cdots,B-1\}\) ， \(i\in[0,M)\) 。

第二步，先子集容斥得到 \(f'[S][i]\) 表示只用 \(S\) 集合数字且每个数都至少用一次答案，对 \(\mathrm{lcm}(S)\;\Big|\;i\) 的项求和得到 \(ans'[S]\) ，再子集容斥(高维前缀和)得到 \(ans[S]\) ，即为答案数组。

第二步的时间复杂度为 \(O(B\cdot2^{B-1}\cdot M)\) ，关键在于第一步。

算法1

对于固定的集合 \(S\) ，考虑DP计算 \(f[S][0\sim M-1]\) 。

设 \(f_n[0\sim M-1]\) 是 \(n\) 位数的答案，从 \(f_n\) 转移到 \(f_{n+1}\) 只需枚举最低位是 \(x\) ： \(f_{n+1}[i\cdot B+x]\leftarrow f_n[i]\) 。

每一次DP的时间复杂度 \(O(NMB)\) ，总时间复杂度 \(O(2^{B-1}\cdot NMB)\) 。

可以通过子任务1,2,3,4，拿到 \(58\) 分。

算法2

用倍增的方式优化算法1中的状态转移： \(f_{2n}\leftarrow f_n\) 。

设 \(b=B^n\%M\) ，则转移方程 \(f_{2n}[i\cdot b+j]\leftarrow f_n[i]\times f_n[j]\) 。

每次倍增DP，时间复杂度 \(O(M^2\log N)\) ，总时间复杂度 \(O(2^{B-1}\cdot M^2\log N)\) 。

当 \(B\leq 6\) 时 \(M\leq60\) 才能跑出来，可以通过子任务 \(5\) 。

结合算法1可以拿到 \(71\) 分。

算法3

注意到算法2的状态转移可以写成循环卷积，可以强行MTT优化将 \(O(M^2)\) 变成 \(O(M\log M)\) 但常数很大。

或者注意到 \(P=999,999,001\) ，而 \(P-1\) 是 \(2520\) 的倍数，所以可以直接跑长度为 \(M\) 的DFT(或者叫做z变换)。

一次DFT的时间复杂度是 \(O(M\cdot质因数之和)\) ，当 \(M=2520\) 时 \(2+2+2+3+3+5+7=24\) 。

无论用MTT还是DFT，每次倍增 \(O(M\log M\log N)\) ，总时间复杂度 \(O(2^{B-1}M\log M\log N)\) 。

当 \(B\leq 8\) 即 \(M\leq 420\) 时可以跑得动，可以通过子任务6。结合算法1可以拿到 \(85\) 分。

对于 \(B=10,N=10^9\) 的极限数据，nodgd本机大约跑 \(14\sim 15\) 秒，估计CF上需要 \(20+\) 秒。（本来大多数情况下CF评测机应该比nodgd的电脑略快一点，但此题似乎CF评测机跑得非常慢，原因未知）。

算法4

观察倍增DP转移 \(f_{2n}[i\cdot b+j] \leftarrow f_n[i]\times f_n[j]\) ，算法3中会先执行 \(g[i\cdot b]\leftarrow f_n[i]\) ，再计算循环卷积 \(g\otimes f_n\) 。

因为DFT计算卷积时本身就是以长度为 \(M\) 进行循环卷积，所以对于连续的卷积 \(a\otimes b\otimes c\cdots\) ，可以分别执行DFT后点乘到一起然后只执行一次IDFT。

如果倍增过程中没有 \(g[i\cdot b]\leftarrow f_n[i]\) 这样的会打乱数组的操作，那么 \(f_n=f_1\otimes f_1\otimes \cdots\) ，就只需要先DFT然后对每个点值快速幂再IDFT，省略中间的所有DFT,IDFT操作，每次DP复杂度降低到 \(O(M\log M+M\log N)\) 。

为了避免倍增时的打乱数组的操作，可以将 \(N\) 个数位 \(B^0,B^1\cdots,B^{N-1}\) 按 \(\bmod M\) 的数值进行分组，同组一起计算就不会有打乱数组的操作。每组算完之后需要先执行打乱数组的操作，再循环卷积乘在一起。

这样做时间复杂度将会和具体分成多少组有关。显然 \(B^n\bmod M\) 是个 \(\rho\) 型，在 \(B=10\) 时共分成了 \(9\) 组，是最多的。具体如下表：

算法3在 \(N=2^{29}-1\) 时达到最坏情况，倍增需要执行 \(29\) 轮且每轮都要 \(2\) 次DFT和 \(2\) 次IDFT，合计 \(29\times 4=116\) 次。

如果换成这个跟组数相关的算法，每组依次是DFT、每个点值快速幂、IDFT、打乱、DFT，然后合并完所有组之后再进行一次IDFT，当 \(B=10\) 时合计 \(9\times 3+1=28\) 次，速度明显提升。

极限数据下，在nodgd本机测试大约需要 \(3.0\sim 3.5\) 秒，虽然仍然无法通过此题但已经很有希望，接下来进行亿点点常数优化。

算法4-优化1

FFT/NTT可以预处理单位根数组，可以写成“先换位置再蝶形变换”的非递归形式，显然DFT也可以做到。

因为每次DFT长度固定，所以可以先预处理每个数会被换到哪个位置。当 \(M=2520\) 时，假设蝶形变换按照质因数递增顺序执行，即 \([2,2,2,3,3,5,7]\) 的顺序，则 \(a_i\) 被换到 \(a'_j\) 的位置关系可以由变进制数的方式计算：将 \(i\) 写成 \([7,5,3,3,2,2,2]\) 顺序的变进制数 \(\overline{i_1i_2i_3i_4i_5i_6i_7}\) ，再把 \(\overline{i_7i_6i_5i_4i_3i_2i_1}\) 按 \([2,2,2,3,3,5,7]\) 的顺序识别得到 \(j\) 。

因为IDFT可以由先DFT，再乘 \(M\) 的逆元，最后 reverse(a+1,a+M) 来实现，所以只需要处理DFT的单位根。为了方便，可以对 \(7\) 轮蝶形变换分别预处理单位根。

在蝶形变换过程中，过多的取模运算是程序慢的主要原因。普通的写法在以 \(p(\in\{2,3,5,7\})\) 为基进行一轮蝶形变换时取模次数为 \(2p\cdot M\) ，其中一半在算 int+=int*int ，另一半在算单位根(或者预处理的单位根的数组下标)。对于 int+=int*int ，容易发现只会加 \(p\) 次，值域上限是 \(p\cdot P^2\) ，用可以 long long 全部加起来之后只取模一次。对于单位根的运算，设蝶形块长为 \(m(\leq M/p)\) ，这一轮蝶形变换需要用到的单位根是 \(\omega_{m\cdot p}^0\sim \omega_{m\cdot p}^{m\cdot p\cdot p}\) 。只要在预处理单位根预处理 \(p\) 个周期即可避免取模，这样预处理的数组大小为 \(蝶形轮数\times M\times p_{max}\) ，最大值为 \(7\times 2520\times 7\) ，几乎可以忽略。

采取这样“精打细算”的方式，可以将每一轮蝶形变换的取模次数精准的控制在 \(M\) 次。当 \(M=2520\) 时， 一次DFT执行的加法、乘法次数为 \(O(M\sum p)=O(24M)\) ，取模次数为 \(7M\) ，常数较小。

算法4-优化2

设 \(b=B^n\bmod M\) ，这样的数位有 \(t\) 个，可以用二元组 \((b,t)\) 来描述这个组。

因为 \(B^n\) 是 \(\rho\) 型结构，所以无论 \(N\) 多大，未进入循环的数位至多只会出现一次，即 \((b_k,t_k=1)\) ，此时用DFT来计算反而慢，而采用算法1暴力DP反而更快。

事实上可以设定一个阈值 \(tm\) （nodgd设的 \(10\) ），若这组的 \(t_k\leq tm\) 则用算法1暴力DP，否则使用DFT。

值得一提的是，只要设置了这个阈值就会有优化效果，和 \(tm\) 的具体值没有太大关系。因为当 \(N\) 较大时进入循环后的每组都会有大量的数位（循环节最大是 \(6\) 所以一组的数位有 \(\approx N/6\) 个），此时 \(tm=1\) 和 \(tm=10^5\) 都没有任何区别。nodgd设为 \(10\) 只是为了让前4个子任务大多数时候能使用暴力DP（说白了就是随便设的）。

当 \(B=10\) 时数位分成 \(9\) 组，优化前全部使用DFT/IDFT合计 \(28\) 次，优化后只有 \(6\) 组使用DFT/IDFT合计 \(6\times 3+1=19\) 次，而暴力DP的时间几乎可以忽略不计。

算法4-优化3

设当前集合为 \(S\) 。

计算每个组 \((b_k,t_k)\) 时，首次DFT的数组都是相同的，即 \(\forall x\in S,\ fa[x\%M]+=1\) 。每一组都会先对这个 \(fa\) 数组进行DFT，每个点值快速幂再IDFT。

不难想到，可以将这一步的DFT结果保存下来，供下一组 \((b_{k+1},t_{k+1})\) 使用。另外，这个 \(g\) 数组中只有至多 \(B-1\) 项非零，所以暴力DFT效率会优于普通DFT（实测要暴力快得多，耗时几乎可忽略）。此时，普通DFT/IDFT的次数从 \(19\) 次再次被降低到了 \(6\times 2+1=13\) 次。

更进一步， \(\rho\) 型结构进入循环后，每个组的数位数量 \(=(N-开头有几组)/周期长度\) ，向上或向下取整，也就是说，当 \(t_k>tm\) 时至多只有两种取值。当 \(t_k\) 相同时，每个点值快速幂再IDFT后的结果都是相同的。 因此，点值快速幂和IDFT也可以减少到只需要 \(2\) 次。此时，普通DFT/IDFT的次数从 \(13\) 次进一步降低到 \(2+6+1=9\) 次，每次 \(O(M\sum p)\) ；而每个点值进行快速幂也从 \(6\) 次降低到 \(2\) 次，每次 \(O(M\log N)\) 。

优化到这个程度，nodgd本机测试在不开O2的情况下需要 \(1.5\sim 1.6\) 秒，CF上可以 \(1900^+ms\) 勉强卡过。

算法4-优化4

对于每个不同的集合 \(S\) ，上面的算法对 \(\forall x\in S,\;fa[x\%M]+=1\) 的数组 \(fa\) 进行DFT+点值快速幂+IDFT。

如果将数组 \(fa\) 视为一个多项式 \(Fa(x)=\sum fa[i]\cdot x^i\) ，则这一步相当于计算 \(Fa(x)^{t_k}\) ，而循环卷积的本质是对 \(P(x)=x^M-1\) 取模。

显然 \(S\) 与 \(Fa(x)\) 是一一对应的，若 \(S\) 不相同则 \(Fa(x)\) 也不相同。但是可以考虑：哪些不同的集合 \(S_1,S_2\) ，他们的 \(Fa_1(x),Fa_2(x)\) 相似，使我们只要算出 \(Fb_1(x)=Fa_1(x)^{t_k}\) ，就可以快速算出 \(Fb_2(x)\) 。

多项式的幂运算只有单项式可以轻易的穿入穿出。也就是说，只有当 \(Fa_2(x)=x^r\cdot Fa_1(x)\) 时我们可以通过 \(Fb_2(x)=x^{r\cdot t_k}\cdot Fb_1(x)\) 来简化计算。对于循环卷积而言，就是对 \(fb\) 的数组下标增加一个偏移量。

什么时候 \(Fa_2(x)=x^r\cdot Fa_1(x)\) 呢？如果集合采用二进制表示，则充要条件是 \(S_2=2^rS_1\) 。

因此我们可以只对 \(1\in S\) 的集合执行“暴力DFT+点值快速幂+IDFT”的操作计算 \(fb\) 序列。而对于 \(1\notin S\) 的集合，借用 \(S_1=S/2^r\) 的 \(fb\) 序列来计算自己的 \(fb\) 序列。这一步可以和打乱操作同时执行从而几乎不增加运行时间： \(fc_2[(i+r\cdot t_k)\cdot b_k]\leftarrow fb_1[i]\) 。

这样可以将暴力DFT、点值快速幂、IDFT的次数减少一半。加上这个优化后，nodgd本机在不开O2的情况下需要 \(1.2\sim 1.3\) 秒，CF上 \(1400^+ms\) 。

算法4-优化5

其实还有一个最强的优化。

设IDFT结束后的数组是 \(fb[0\sim M-1]\) ，打乱过程是 \(fc[i\cdot b_k]\leftarrow fb[i]\) ，如果理解成多项式则 \(Fb(x)=Fa(x)^{t_k}\) ，而 \(Fc(x)=Fb(x^{t_k})\) 。

考虑优化掉“IDFT+打乱+DFT”的过程，即在 \(fb\) 的点值序列上进行打乱操作： \(Fc(\omega_M^i)=Fb(w_M^{i\cdot b_k})\) 。这样就完全避免了中途的DFT/IDFT操作，对于每个 \(S\) 只跑一个 \(fa\) 数组的暴力DFT和最终 \(f[S][0\sim M-1]\) 的IDFT。

将打乱操作放在点值序列上之后，如果要结合优化4，还需要将 \(fb\) 加偏移量的操作也放在点值序列上：

假设已经对于 \(1\in S_1\) 的集合 \(S_1\) 算出了 \(fb\) 的点值序列 \(Fb(\omega_M^0)\sim Fb(\omega_M^{M-1})\) 。

考虑 \(S_2=2^rS_1\) ，因为 \(Fb_2(x)=x^{r\cdot t_k}\cdot Fb_1(x)\) ，所以 \(Fb_2(\omega_M^i)=\omega_M^{i\cdot r\cdot t_k}\cdot Fb_1(\omega_M^i)\) 。

与打乱操作同时执行则是 \(Fc_2(\omega_M^i)=Fb_2(\omega_M^{i\cdot b_k})=\omega_M^{i\cdot b_k\cdot r\cdot t_k}\cdot Fb_1(\omega_M^{i\cdot b_k})\) 可以避免额外增加运行时间。

此时已将刚才的 \(9\) 次DFT/IDFT再次降低到 \(1\) 次，不再是唯一瓶颈。DFT/IDFT部分、点值快速幂、打乱，这三类操作的时间复杂度相平衡。

至此，本机不开O2跑 \(0.5\sim 0.6\) 秒，开启O2跑 \(0.3\) 秒，CF跑 \(900^+ms\) 。