诈骗题说是，但话说真的有快速阶乘和算法说是，建议 $n\le 10^{18}$。

注意到 $i\ge 10000$ 时，$i!$ 必然能被 $10000$ 整除，后面的部分不用算即可，复杂度为 $O(1)$。

fhx 跟我说这题是简单题我就一眼秒了，没想到意思发现了标解以外的在线算法。

假设 $n,q,V$ 同阶。

对 $a$ 分块，块长为 $B$，散块部分直接暴力计算，复杂度为 $O(B)$。

对于整块的部分，考虑对询问的 $k$ 进行阈值分治，令阈值也为 $B$，查询的整块区间为 $L\sim R$ 块。

若 $k\le B$，预处理 $f_{i,j}$ 表示 $1\sim i$ 块中所有数除以 $k$ 下取整之和，则这部分就是 $f_{R,k}-f_{L-1,k}$。

若 $k> B$，枚举 $v\in [0,\lfloor\frac{V}{B}\rfloor]$ 表示存在多少个 $\lfloor\frac{a_i}{k}\rfloor=v$ 在整块区间内，转化为查询多少个值域在 $[vB,(v+1)B)$ 的数在整块区间里面，预处理 $g_{i,j}$ 表示前 $i$ 块 $\le j$ 的数的个数，也可以 $O(1)$ 求。

直接取 $B=\sqrt{n}$ 可以得到理论最优复杂度为 $O(n\sqrt{n})$，空间为 $O(n\sqrt{n})$。实际上对于第二问转化出来的询问离线处理也可以做到空间线性。

别看我，代码是 ds 写的。

分析

构造方法是：

1.  对于每个点，随机选择一个比其标号小的点作为父亲，生成一棵**递增标号树**。

2. 对上述树的标号打乱。

不难想到，第一步产生的数则种类有 $\prod_{i=2}^{N} (i-1)$ 个，第二步的排列数有 $N!$ 个可能。所以，总方案数为：

$$\prod_{i=2}^{N} (i-1) \times N!$$

关键点来了：有多少种情况可以生成这个树。反向考虑这棵树能怎么生成，也就是存在一种情况得到了一个递增树和一个排列，可以得到这个树。显然，如果一个递增树可以加上某个排列形成这棵树，那么这个排列是唯一确定的。那么，任意一个与输入的树结构相同的递增标号树都可以得到这个树，做一个贡献。如何计算这个贡献？首先，由定义可知 `1` 号节点一定是根，令其 $k$ 个子节点的子树大小分别为 $s_1, s_2, \dots, s_k$，则每个子树分配点集方案树为

$$C_{n-1}^{s_1} \times C_{n-1-s_1}^{s_2} \times \dots \times C_{n-1-\sum_{i=1}^{k-1} s_i}^{s_k}$$

化简可得

$$\frac{N!}{s_1! \times s_2! \times \dots \times s_k!}$$

每个子树同样满足这个标号规则，令 $f(T)$ 为以 $T$ 为根的子树的方案，则：

$$f(T) = \frac{N!}{\prod s_i!} \times \prod f(子树i)$$

展开，化简可得：

$$f(T) = \frac{N!}{\prod sz_u}$$

其中，$sz_u$ 表示树上每个节点的子树大小。这样子我们可以求出来一个节点为根的方案树，接下来我们需要换根。

令 $f_u = \frac{N!}{\prod sz_x},\quad f_v = \frac{N!}{\prod sz'_x}$，从 $u$ 到 $v$ 只会有以下树改变：

$$sz_u = N,\quad sz_v = sz_v$$

$$sz_u' = N-sz_v,\quad sz_v'=N$$

分母产生的变化

$$D_u=(\prod_{x \ne u,v} sz_x) \times N \times sz_v$$

$$D_v = (\prod_{x\ne u,v} sz_x) \times (N-sz_v) \times N$$

比一下，

$$\frac{f_v}{f_u} = \frac{D_u}{D_v} = \frac{N \times sz_v}{(N-sz_v) \times N} = \frac{sz_v}{(N-sz_v)}$$

所以，$f_v = f_u \times \frac{sz_v}{(N-sz_v)}$。

之后我们将其求和，得到总答案数，再逆元除以总情况数即可。

参考代码

#include <cstring>
#include <iostream>
using namespace std;
const int N = 2e5 + 5, mod = 1e9 + 7;
int hd[N], nxt[N * 2], to[N * 2], num = 1;
void add(int u, int v) {
    num++;
    to[num] = v;
    nxt[num] = hd[u];
    hd[u] = num;
}
int sz[N];
void dfs(int u, int fa) {
    sz[u] = 1;
    for (int i = hd[u]; i; i = nxt[i]) {
        int v = to[i];
        if (v == fa) continue;
        dfs(v, u);
        sz[u] += sz[v];
    }
}
long long fpow(long long x, long long y) {
    long long res = 1;
    while (y) {
        if (y & 1) res = res * x % mod;
        x = x * x % mod;
        y >>= 1;
    }
    return res;
}
int n;
long long f[N];
void dfs2(int u, int fa) {
    for (int i = hd[u]; i; i = nxt[i]) {
        int v = to[i];
        if (v == fa) continue;
        f[v] = f[u] * sz[v] % mod * fpow(n - sz[v], mod - 2) % mod;
        dfs2(v, u);
    }
}
int main() {
    int t;
    cin >> t;
    while (t--) {
        cin >> n;
        num = 1;
        memset(hd, 0, sizeof(hd));
        for (int i = 1; i <= n - 1; i++) {
            int u, v;
            cin >> u >> v;
            add(u, v);
            add(v, u);
        }
        dfs(1, 0);
        long long fz = 1;  // 分子
        for (int i = 1; i <= n; i++) {
            fz = fz * i % mod;
        }
        long long fm = 1;  // 分母
        for (int i = 1; i <= n; i++) {
            fm = fm * sz[i] % mod;
        }
        f[1] = fz * fpow(fm, mod - 2) % mod;
        dfs2(1, 0);
        long long ans = 0;
        for (int i = 1; i <= n; i++) {
            ans += f[i];
            ans %= mod;
        }
        fm = 1;
        for (int i = 1; i <= n - 1; i++) fm = fm * i % mod;
        for (int i = 1; i <= n; i++) fm = fm * i % mod;
        cout << ans * fpow(fm, mod - 2) % mod << endl;
    }
    return 0;
}

一个组合计数比较好的题。

恰好一次都不出现的方案数，其实可以容斥成“出现至少 $k$ 次”的方案数（下文表示为 $f_k$）,具体的，可以表示为 $\sum_{i=0}(-1)^if_{i}$，简单容斥不在阐述。

对于至少出现 $k$ 次的方案数，其实可以先生成 $k$ 个段，然后插入到最终的序列里，把每个段缩成一个位置，原本有 $n$ 个位置，$4k$ 个位置被缩进段了，$k$ 个段缩成了一个位置，还剩下 $n-4k+k=n-3k$ 个位置，要选 $k$ 个位置，就是 $C_{n-3k}^k$。

对于剩下的部分，其实就是求出将四种数填满 $n-4k$ 个位置，且剩下每种数都不超过 $a/b/c/d-k$ 次。下文方便讨论，$n\gets n-4k,a,b,c,d$ 各自减去 $k$。

对于 $a+b+c+d=n$ 的情况，实际上是多重集的排列数，即 $\frac{n!}{a!b!c!d!}$。

对于 $a+b+c+d\ge n$ 的情况，实际上可以枚举 $a'\le a,b'\le b,c'\le c,d'\le d$，满足 $a'+b'+c'+d'=n$，然后产生 $\frac{n!}{a'!b'!c'!d'!}$ 的方案，实际上做四次卷积即可解决，这个做法需要外层枚举一个 $k$，内层做卷积，复杂度为 $O(n^2\log n)$。

实际上一个更优的做法，枚举 $x=a'+b'$，则 $n-x=c'+d'$，枚举 $y=a'$，则 $x-y=b'$，注意到 $a'\in [0,a],b'\in [0,b]$，则对于 $y$ 合法取值的限制实际上是和 $x$ 相关的某个区间，可以 $O(1)$ 算法，对于 $c'+d'$ 的部分同理，所以只需要枚举一个 $k$，枚举一个 $x$，然后预处理出组合数的前缀和即可 $O(1)$ 回答，时间复杂度为 $O(n^2)$。

//n ^ 3
#include<iostream>
using namespace std;

#define int long long
const int M = 998244353;
const int N = 1005;
int C[N][N];
int f[5][N];
int n, a, b, c, d;

void init(){
	for(int i = 0;i < N;i ++){
		C[i][0] = 1;
		for(int j = 1;j <= i;j ++){
			C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % M;
		}
	}
}

signed main(){
	cin.tie(0) -> ios::sync_with_stdio(0);
	cin >> n >> a >> b >> c >> d;
	init();
	int ans = 0;
	
	int limit = min(a, min(b, min(c, min(d, n / 4))));
	
	for(int i = 0;i <= limit;i ++){
		int rev = n - 4 * i;
		int lim[5] = {0, a - i, b - i, c - i, d - i};
		for(int j = 0;j <= 4;j ++){
			for(int k = 0;k <= rev;k ++){
				f[j][k] = 0;
			}
		}
		f[0][0] = 1;
		for(int k = 1;k <= 4;k ++){
			for(int j = 0;j <= rev;j ++){
				for(int x = 0;x <= min(j, lim[k]);x ++){
					f[k][j] = (f[k][j] + f[k - 1][j - x] * C[j][x] % M) % M;
				}
			}
		}
		int way = C[n - 3 * i][i] * f[4][rev] % M;
		if(i & 1) ans = (ans - way + M) % M;
		else ans = (ans + way) % M;
	}
	cout << ans << '\n';
	return 0;
}
//n ^ 2 log n
#include<iostream>
#include<vector>
#include<algorithm>

using namespace std;

#define ll long long
const int MOD = 998244353;
const int MAXN = 4005;
const int G = 3;
const int Gi = 332748118;

ll n, num[4];
ll fact[MAXN], inv[MAXN];
int rev[MAXN << 2];
ll A[MAXN << 2], B[MAXN << 2], Poly[MAXN << 2];

ll qpow(ll a, ll b){
	ll res = 1;
	while(b){
		if(b & 1) res = res * a % MOD;
		a = a * a % MOD;
		b >>= 1;
	}
	return res;
}

void NTT(ll A[], int len, int op){
	for(int i = 0;i < len;i ++) if(i < rev[i]) swap(A[i], A[rev[i]]);
	for(int mid = 1;mid < len;mid <<= 1){
		ll Wn = qpow(op == 1 ? G : Gi, (MOD - 1) / (mid << 1));
		for(int i = 0;i < len;i += (mid << 1)){
			ll wk = 1;
			for(int j = 0;j < mid;j ++, wk = wk * Wn % MOD){
				ll x = A[i + j], y = wk * A[i + j + mid] % MOD;
				A[i + j] = (x + y) % MOD;
				A[i + j + mid] = (x - y + MOD) % MOD;
			}
		}
	}
	if(op == -1){
		ll invL = qpow(len, MOD - 2);
		for(int i = 0;i < len;i ++) A[i] = A[i] * invL % MOD;
	}
}

void NTT_init(int len){
	int L = 0; while((1 << L) < len) L ++;
	for(int i = 0;i < len;i ++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (L - 1));
}

void init(){
	fact[0] = 1;
	for(int i = 1;i < MAXN;i ++) fact[i] = fact[i - 1] * i % MOD;
	inv[MAXN - 1] = qpow(fact[MAXN - 1], MOD - 2);
	for(int i = MAXN - 2;i >= 0;i --) inv[i] = inv[i + 1] * (i + 1) % MOD;
}

ll C(int n, int m){
	if(m < 0 || m > n) return 0;
	return fact[n] * inv[m] % MOD * inv[n - m] % MOD;
}

ll NFT(int k){
	int m = n - 4 * k;
	if(m < 0) return 0;
	int len = 1;
	while(len <= 4 * m) len <<= 1;
	NTT_init(len);
	for(int i = 0;i < len;i ++) Poly[i] = 1;
	for(int i = 0;i < 4;i ++){
		int limit = num[i] - k;
		for(int j = 0;j < len;j ++) A[j] = (j <= limit && j <= m) ? inv[j] : 0;
		NTT(A, len, 1);
		for(int j = 0;j < len;j ++) Poly[j] = Poly[j] * A[j] % MOD;
	}
	NTT(Poly, len, -1);
	return C(n - 3 * k, k) * Poly[m] % MOD * fact[m] % MOD;
}

int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);
	init();
	cin >> n >> num[0] >> num[1] >> num[2] >> num[3];
	ll ans = 0;
	int lim_k = min({(ll)n / 4, num[0], num[1], num[2], num[3]});
	for(int k = 0;k <= lim_k;k ++){
		ll res = NFT(k);
		if(k & 1) ans = (ans - res + MOD) % MOD;
		else ans = (ans + res) % MOD;
	}
	cout << ans << '\n';
	return 0;
}