幻化成风

有一个长为 \(m\) 的序列 \(\{a_i\}\) 和一个 \(n\)，其中 \(n!\) 可以用下列方式表示：

\[n!=\prod {b_i}^{a_i} \]

其中 \(\{b_i\}\) 中的数两两不同。对不同表达方式计数。\(\bmod 10^9+7\)。\(1\le n\le 10^4,1\le m,\sum a\le 30\)。

两种表示方式不同当且仅当集合 \(\{(b_i, a_i)\}\) 不同。现在你需要对不同表示方式计数。答案对 \(10^9+7\) 取模。

考虑题中的两个限制怎么办。\(\{b_i,a_i\}\) 本质不同在 \(b_i\) 互不相同的情况下其实就是对于相同的 \(a_i\) 满足 \(b_i\) 无序。那我们直接当成无序做，然后除去 \(\prod cnt(a_i)\) 即可。

考虑 \(b_i\) 互不相同的限制。这让这道题非常难以直接下手。这个时候考虑使用 集合不相等容斥，枚举全集的划分 \(\mathcal S\)，满足 \(\mathcal S\) 内任意集合的所有元素全相等，那么答案就是 $\sum\limits_{\mathcal S}f(\mathcal S)\prod\limits_{s\in \mathcal S} (|s|-1)!\times (-1)^{|s|-1} $。其中 \(f(\mathcal S)\) 表示这种集合划分的对应答案。后面的是容斥系数。

直接枚举 \(\mathcal S\) 然后算当然相当寄，考虑我们计算 \(f(\mathcal S)\) 和后面的容斥系数需要什么：需要 \(\mathcal S\) 中所有元素 \(s\) 的 \(\sum\limits_{u\in s} a_u\) 和 \(|s|\)。而不需要关系 \(a\) 具体的划分情况。

这样的话，我们爆搜出来所有的本质不同的 \(\mathcal S\)，发现在 \(a_i=1\) 的极限情况下等价于 \(\mathrm{bell}(\sum a)\)，也即 \(\sum a\le 30\) 的拆分数个数。而我们有 \(\mathrm{bell}(\sum a)\le \mathrm{bell}(30)=5604\)，记 \(T=5604\)，现在对着每一种状态算原先的答案。

因为质因子互相独立，我们算出来 \(p_i\) 在 \(n!\) 内出现的次数 \(lim_i\)，那么就相当于一个完全背包，\(\mathcal O(nm)\) 做一遍即可。时间复杂度 \(\mathcal O(nmT)\)。

#include <bits/stdc++.h>
#define pb emplace_back
#define fir first
#define sec second

using i64 = long long;
using u64 = unsigned long long;
using pii = std::pair<int, int>;

constexpr int maxn = 1e4 + 5, maxm = 30 + 5;
constexpr int mod = 1e9 + 7;
constexpr u64 base = 19260817;
void add(int& x, int y) { if ((x += y) >= mod) x -= mod; return ; }
void sub(int& x, int y) { if ((x -= y) < 0) x += mod; return ; }
int inc(int x, int y) { return (x + y) >= mod ? (x + y - mod) : (x + y); }
int dec(int x, int y) { return (x - y) < 0 ? (x - y + mod) : (x - y); }
int power(int x, int y) {
	int ans = 1;
	for (; y; y >>= 1, x = (i64)x * x % mod) if (y & 1) ans = (i64)ans * x % mod;
	return ans;
}
int fac[maxm], ifac[maxm], ans, a[maxm], n, m, prime[maxn], cnt, lim[maxn], f[maxn], buc[maxm];
bool flg[maxn];
u64 pw[maxm << 1];
int C(int n, int m) {
	if (n < 0 || m < 0 || n < m) return 0;
	return (i64)fac[n] * ifac[n - m] % mod * ifac[m] % mod;
}
void init(int n) {
	for (int i = fac[0] = 1; i <= n; ++i)
		fac[i] = (i64)fac[i - 1] * i % mod;
	ifac[n] = power(fac[n], mod - 2);
	for (int i = n - 1; ~i; --i)
		ifac[i] = (i64)ifac[i + 1] * (i + 1) % mod;	
	for (int i = pw[0] = 1; i <= n * 2; ++i)
		pw[i] = (u64)pw[i - 1] * base;
	return ;
}
void sieve(int n) {
	for (int i = 2; i <= n; ++i) {
		if (!flg[i]) prime[++cnt] = i;
		for (int j = 1; j <= cnt && i * prime[j] <= n; ++j) {
			flg[i * prime[j]] = true;
			if (i % prime[j] == 0) break ;
		}
	}
	for (int i = 1; i <= cnt; ++i) {
		int x = n;
		while (x) lim[i] += x / prime[i], x /= prime[i];
	}
	return ;
}
int sgn(int x) { return (x & 1) ? (mod - 1) : 1; }

struct info {
	pii a[maxm];
	u64 hash;
	info() { memset(a, 0, sizeof(a)); hash = 0; }
	void gethash() {
		std::sort(a, a + m, std::greater<pii>()); hash = 0;
		for (int i = 0; i < m && a[i].fir; ++i)
			hash += a[i].fir * pw[i] + a[i].sec * pw[i + 30];
		return ;
	}
	bool operator < (const info& rhs) const {
		return hash < rhs.hash;
	}
} tmp;
std::map<info, int> dp[2];
std::map<std::vector<int>, int> coef;

int main() {
	freopen("count.in", "r", stdin);
	freopen("count.out", "w", stdout);
	scanf("%d %d", &n, &m);
	init(30); sieve(n);
	for (int i = 1; i <= m; ++i)
		scanf("%d", &a[i]);
	std::sort(a + 1, a + 1 + m);
	int now = 0; dp[now][info()] = 1;
	for (int i = 1; i <= m; ++i) {
		now ^= 1; dp[now].clear();
		for (auto& [key, val] : dp[now ^ 1]) {
			for (int k = 0; k < m && (k == 0 || key.a[k - 1].fir > 0); ++k) {
				tmp = key; tmp.a[k].fir += a[i];
				++tmp.a[k].sec; tmp.gethash();
				add(dp[now][tmp], val);
			}
		}
	}
	for (auto& [key, val] : dp[now]) {
		std::vector<int> g; int w = val; 
		for (int k = 0; k < m && key.a[k].fir; ++k)
			g.pb(key.a[k].fir), w = (i64)w * fac[key.a[k].sec - 1] % mod * sgn(key.a[k].sec - 1) % mod;
		add(coef[g], w); g.clear();
	}
	for (auto& [key, val] : coef) {
		memset(f, 0, sizeof(f)); f[0] = 1;
		for (auto& v : key)
			for (int j = v; j <= n; ++j)
				add(f[j], f[j - v]);
		int res = val;
		for (int i = 1; i <= cnt; ++i)
			res = (i64)res * f[lim[i]] % mod;
		add(ans, res);
	}
	for (int i = 1; i <= m; ++i)
		++buc[a[i]];
	for (int i = 1; i <= 30; ++i)
		if (buc[i]) ans = (i64)ans * ifac[buc[i]] % mod;
	printf("%d\n", ans);
	return 0;
}