官方题解只能做到 $O(4^p p^4 + Tp^2 \log^2_p n)$？玩的太差了。设 $q=(p-1)/2$、$S=\binom{3q}{q}$，$L$ 为询问中 $n$ 的最大 $p$ 进制位数。对奇素数 $p$，包含全部预处理的总时间为 $O(p^5S+p^6+Tp^2L\min(p^2,L))$，峰值空间为 $O(p^3S+p^6)$。对 $p=2$，读入后的计算只需要一次位运算。这些界按模 $p$ 的算术运算次数计算。
从图计数到数位递推
记 $R(n,d)$ 为 $n$ 个有标号顶点的简单无向 $d$-正则图数量。题目要求计算 $R(n,d)\bmod p$，其中 $n\le10^{18}$，原题的 $p\in\{2,3,5,7\}$；下面把 $p$ 当作参数分析，而不是把指数开销当作常数。
固定一个顶点，它有 $\binom{n-1}{d}$ 种邻居集合。对任意两个这样的集合，重新排列其余顶点，就得到两类合法图之间的双射。因此 $R(n,d)$ 是 $\binom{n-1}{d}$ 的整数倍。由 $(1+x)^{n-1}\equiv\prod_i(1+x^{p^i})^{(n-1)_i}\pmod p$，只要存在一位满足 $d_i>(n-1)_i$，答案就是 $0$。这只是必要条件，通过判零不代表答案非零。
将 $n$ 写成 $\sum_i n_ip^i$。对每个 $i$，划出 $n_i$ 个大小为 $p^i$ 的顶点块，并在各块内独立进行循环置换。这些置换组成一个阶为 $p$ 的幂的群。它作用在所有合法图上，轨道大小也是 $p$ 的幂，所以模 $p$ 计数时，只需要保留被整个群固定的图。这里数的仍是有标号图，不是在数同构类，也没有除以群的阶。分块只固定一次，不要乘上分块方案数。
两个不同块之间的边，在两边独立循环置换的作用下构成一个轨道。因此，这两个块之间要么没有边，要么是完全二分图。
先考虑奇素数。对大小为 $s=p^i$ 的块，内部每一种无向循环距离都贡献一个可选边集，选中以后所有顶点的度数增加 $2$。内部度数生成函数为 $F_s(x)=(1+x^2)^{(s-1)/2}$。令 $f(x)=(1+x^2)^q$、$f_e=[x^e]f(x)$，并把区间 $[0,p-1]$ 以外的 $f_e$ 定义为 $0$。由于 $(s-1)/2$ 的低 $i$ 个 $p$ 进制数位都是 $q$，有 $F_s(x)\equiv\prod_{j< i}f(x^{p^j})\pmod p$。因而，对 $0\le t< s$，系数 $[x^t]F_s(x)$ 就是 $t$ 的各位对应的 $f_e$ 的乘积。
考虑加入第 $i$ 位的块。记 $s=p^i$、$r=n\bmod s$、$t=d\bmod s$、$a=n_i$、$b=d_i$。已有 $r< s$ 个旧顶点，现在加入 $a$ 个大小为 $s$ 的新块。
当 $r>0$ 时，通过前面的判零可知 $t< r$。旧顶点在旧顶点集合内的度数必须为 $t$，而每连接一个新块，度数增加 $s$。因此，每个旧块都必须恰好连接 $b$ 个新块。这个要求与旧图内部的具体边集无关，所以这一位的扩展方案数可以直接乘到已有方案数上。当 $r=0$ 时没有旧顶点，旧部分取空图，仍可使用下面的生成函数。
用 $e_b(x_1,\ldots,x_a)$ 表示 $b$ 次初等对称多项式。一个大小为 $p^j$ 的旧块，选择连接哪些新块，对应 $e_b(x_1^{p^j},\ldots,x_a^{p^j})$。将所有旧块的贡献相乘，在模 $p$ 下正好得到 $e_b(\boldsymbol x)^r$。加上新块内部的边，得到 $P_i^{a,b}(\boldsymbol x)=e_b(\boldsymbol x)^r\prod_{j=1}^aF_s(x_j)$。
每个变量的次数最多为 $r+s-1< 2s$。要最终得到度数 $bs+t$，这个多项式里需要考虑的次数只有 $t$ 与 $t+s$。也就是说，每个新块只需要记录一个 $0/1$ 进位。
定义 $D_i^{a,b}(h)$：在 $P_i^{a,b}$ 中，把所有恰好有 $h$ 个 $1$ 的向量 $(c_1,\ldots,c_a)$ 对应的单项式 $\prod_jx_j^{t+sc_j}$ 的系数相加。多项式关于所有变量对称，所以只需记进位个数，状态数为 $a+1$。
再定义 $G_{a,b,h}$：在 $a$ 个有标号顶点的小图中，前 $h$ 个顶点的度数为 $b-1$，其余顶点的度数为 $b$ 的方案数。一个进位为 $c_j$ 的新块还需要连接 $b-c_j$ 个其他新块，所以本位扩展系数为
$$ E_i=\sum_{h=0}^a D_i^{a,b}(h)G_{a,b,h}. $$
$D_i^{a,b}(h)$ 已经对进位集合求和，不能再乘 $\binom ah$。最终答案是所有非空数位的 $E_i$ 的乘积。
接下来让所有参数对 $(a,b)$ 共用一次数位扫描。假设当前输入数位为 $m=n_i$、$v=d_i$。对于一组固定的参数 $(a,b)$，从 $P_i^{a,b}$ 到 $P_{i+1}^{a,b}$，新乘上的局部因子为 $e_b(\boldsymbol x)^m\prod_jf(x_j)$，其中变量在完整生成函数中应替换为 $x_j^s$。
若某个变量在 $e_b^m$ 中的次数为 $j$，旧进位为 $c$，新进位为 $c'$，从 $f$ 取得的次数为 $\ell$，就要求 $c+j+\ell=v+pc'$。定义 $M_{a,b}^{m,v}[h,k]$ 为固定一种大小为 $h$ 的旧进位集合后，对所有大小为 $k$ 的新进位集合求和所得的权值。于是有 $D_{i+1}^{a,b}(k)=\sum_hD_i^{a,b}(h)M_{a,b}^{m,v}[h,k]$。
初始化时，所有参数对都有 $D_0^{a,b}(0)=1$，其他状态为 $0$。每一位必须先计算扩展系数，再更新这些前缀状态。不能为每一位重新扫描它之前的数位，否则会产生不必要的 $\log^2 n$。
有两种情况可以直接处理。若 $b=0$，扩展系数是 $([x^t]F_s)^a$。若 $b=a$，旧块与新块全部相连，新块之间也全部相连，内部度数必须为 $s+t-r$；利用 $F_s$ 的系数对称性，它的系数等于 $[x^{r-1-t}]F_s$。此时 $r>0$，并且 $r-1-t$ 正好是 $n-1-d$ 的低位前缀。因此，分别维护 $d$ 和 $n-1-d$ 的各位 $f_e$ 的乘积即可。实际只需维护 $1\le b< a< p$ 的矩阵状态。
单指数预处理，以及进一步缩小的小图子过程
先计算转移矩阵。把 $e_b(x_1,\ldots,x_a)^m$ 看成一张 $m$ 行、$a$ 列的 $0/1$ 矩阵，每行恰好有 $b$ 个 $1$。现在逐列加入矩阵，维护每一行已经选中的位置数。
引入行变量 $\boldsymbol y=(y_1,\ldots,y_m)$，用 $u,z$ 分别记录旧、新进位数量。某列选了 $j$ 行，对应 $e_j(\boldsymbol y)$；该列四种进位组合的总权值为 $A_j(u,z)=f_{v-j}+uf_{v-1-j}+zf_{v+p-j}+uzf_{v+p-1-j}$。令 $\Phi=\sum_{j=0}^m e_j(\boldsymbol y)A_j(u,z)$，则
$$ M_{a,b}^{m,v}[h,k]=\binom ah^{-1}[y_1^b\cdots y_m^b u^h z^k]\Phi^a. $$
这里除以 $\binom ah$，是因为右侧把旧进位的位置也枚举了。由于 $a< p$，逆元存在。新进位集合本来就需要求和，不应再除以或乘以 $\binom ak$。
补集对称性可以把 $b$ 限制到 $b\le a/2$。将矩阵的所有位取反，同时把 $c,c',\ell$ 替换成 $1-c,1-c',p-1-\ell$，利用 $f_\ell=f_{p-1-\ell}$，得到 $M_{a,b}^{m,v}[h,k]=M_{a,a-b}^{m,m-v}[a-h,a-k]$。通过判零后，所有 $m>0$ 的实际输入数位都满足 $v\le m$，因此只需预处理这些转移。最终只需真正计算 $1\le b\le q$。
若 $m=0$，可以直接计算。$v=0$ 时转移是单位矩阵；$v>0$ 时所有新进位均为 $0$，旧进位必须全部等于 $v\bmod2$，对应权值为 $f_{v-(v\bmod2)}^a$。
由于各行可以任意交换，用 $c_r$ 记录当前度数为 $r$ 的行数即可。若当前允许的度数区间宽度为 $w$，共有 $\binom{m+w}{w}$ 个频数状态。状态里存储所有有标号行排列的系数之和，而不是某一种排列的系数。
频数压缩之后，还必须避免枚举所有类别各选多少行。用 $X_r$ 标记度数为 $r$ 的行，一次同时计算所有 $e_j(\boldsymbol y)F$，等价于代换 $X_r\mapsto X_r+zX_{r+1}$，其中 $z$ 记录选中行数。按 $r$ 从大到小执行这些代换。从某一类选出 $t$ 行，频数从 $(c_r,c_{r+1})$ 变成 $(c_r-t,c_{r+1}+t)$，权值乘 $\binom{c_r}{t}$。
这个顺序保证一行不会在同一列中被选中两次。各类别依次处理并合并相同的中间状态，因此不会出现各类别选择数的乘积。对全部 $N=\binom{m+w}{w}$ 个状态，有 $\sum_{\text{状态}}\sum_{r< w}c_r=wmN/(w+1)\le mN$，所以所有非零移动只有 $O(mN)$ 个。加上选中行数，以及 $u,z$ 的次数这三个辅助维度，一次加入列的代价是 $O(p^4N)$。
还应按剩余列数剪枝。处理完 $a$ 列后，后面最多再加入 $p-1-a$ 列。要在后续某个需要的幂次中达到每行度数 $b$，当前度数必须处于 $[\max(0,b-(p-1-a)),\min(a,b)]$。每次只为实际区间建立频数状态，区间下端统一减掉。固定 $(m,b,v)$ 后，连续乘上同一个 $\Phi$，就能同时得到所有需要的 $a$。
小图计数 $G$ 可以进一步降低指数底数：不要逐个加入顶点，而是从目标度数序列出发逐个删除顶点。
固定 $a,b$，把所有 $h$ 作为并行通道。初始时，$h$ 个顶点需要 $b-1$ 条边，其余顶点需要 $b$ 条边。若当前选择删除一个剩余度数为 $r$ 的顶点，先乘上该度数类别的顶点数并删除它，再从其他顶点中选择恰好 $r$ 个作为邻居，将它们的剩余度数减一。
邻居选择仍然使用上述逐类别二项式变换，只是把增加度数改成减少度数，并按相反的类别顺序处理。程序允许每一步选择任意一个剩余顶点，所以每张图被按全部 $a!$ 种删除顺序计数，最终除以 $a!$ 即可。每个 $h$ 通道初始化为一种固定的目标度数分配，所以这里不再除以 $\binom ah$。$a< p$，所用逆元存在。
删除 $k$ 个顶点以后，剩余度数一定属于 $[\max(0,b-1-k),\min(b,a-k-1)]$。一次转移的临时频数表中，若还剩 $m$ 个顶点、度数区间宽度为 $w$，则 $m+w\le a+1$。因此状态数不超过 $\binom{a+1}{m}$，对所有删除层求和只有 $O(2^a)$。使用并行通道和二项式变换，所有小图表可在 $O(p^5 2^p)$ 时间内计算；另一半参数用补图关系 $G_{a,b,h}=G_{a,a-b,a-h}$ 得到。
这个 $2^p$ 是小图子过程的界，不是整份算法的界。
总复杂度、特殊情况与实现
主转移最大的频数空间是 $S=\binom{3q}{q}$。还需要说明为什么逐个枚举各参数之后，总预处理时间是 $O(p^5S)$。
在固定 $b$ 的逐列计算中，第 $a$ 列使用的临时区间宽度为 $w=\min(a,b)-\max(0,b+a-1-2q)$。在所有 $1\le b\le q$、$1\le a\le2q$ 中，宽度 $w$ 出现 $4(q-w)+2$ 次。而对行数求和，有 $\sum_{m=1}^{2q}\binom{m+w}{w}=\binom{2q+w+1}{w+1}-1$。
后一个二项式系数在 $w=q$ 时小于 $3S$；当 $w$ 每减少 $1$，它至多缩小到原来的一半。因此，对全部 $(m,b,a)$ 求和后的工作状态数仍为 $O(S)$。再计入 $v$ 的 $O(p)$ 种取值，以及每个状态 $O(p^4)$ 的处理代价，就是 $O(p^5S)$。完整转移表有 $O(p^6)$ 个元素，其初始化和补集转换也明确计入复杂度。
对单次询问，只需维护该询问实际出现的参数对。若其数量为 $K$，则 $K\le\min(p^2,L)$。每对参数有 $O(p)$ 个状态，每次数位更新使用一个 $O(p^2)$ 大小的矩阵，查询时间为 $O(pL+p^2KL)$。状态在最后一次使用之后可以停止更新。
因此，包含小图计算、转移构造、完整表初始化和全部询问，总时间与空间为
$$ O\!\left(p^5\binom{3q}{q}+p^6+Tp^2L\min(p^2,L)\right),\qquad O\!\left(p^3\binom{3q}{q}+p^6\right). $$
由阶乘的渐近估计，$\binom{3q}{q}=\Theta((3\sqrt3/2)^p/\sqrt p)$。固定 $p$ 后，单次查询对 $n$ 的依赖为 $O(\log_p n)$。但关于 $p$，整体是 $O^*(2.598076\ldots^p)$。
对 $p=2$，每个非空数位只有一个块。大小为 $2^i$ 的循环块，其内部生成函数模 $2$ 为 $1+x+\cdots+x^{2^i-1}$，每个合法内部度数的系数都是 $1$。通过最初的数位判零后，每层扩展系数都为 $1$，所以答案就是 $[(d\mathbin{\&}(n-1))=d]$。
提交记录：https://qoj.ac/submission/2936376 ，QOJ 9 级题这么简单。