【日常训练】迪杂斯特
Problem
在 A 国，一年有 \(50!=1\times2\times\dots\times50\) 天，每天的编号从 \(1\) 到 \(50!\) 。
A 国的土地上经常发生自然灾害。据 A 国科学家的妍究， A 国有 \(n\) 种自然灾害，每种自然灾害的发生都有周期性。
具体地，第 \(i\) 种自然灾害可以用一个只包含 Y 和 N 组成的字符串组成。
设这个字符串的长度为 \(l\) 。如果这个字符串的第 \(j\) 个字符为 Y ，那么这就表示对于所有的 \(k\in\text{N}\) ，在第 \(k\times l+j\) 天会发生第 \(i\) 种自然灾害，否则不会发生。
一天同时发生的自然灾害种数决定了该天整个 A 国境内的危险程度。
所以，你需要对于所有的 \(0\le i\le n\) ，求出一年内有多少天同时发生的自然灾害种数恰好为 \(i\) 。
\(n \leq 30\)，字符串的长度不超过 \(50\)。
有几个部分分：
- 所有字符串长度的乘积不超过 \(1111111\)；
- 所有字符串长度两两互质或相等；
- 所有字符串长度不超过 \(48\)。
Solution
校内模拟赛的题，被锤爆了，深有体会。
部分分 1
将第 \(i\) 个字符串的长度称为 \(l_i\)。
将所有串的长度的 \(\mathrm{lcm}\) 记为 \(L\)。
那么我们可以考虑将所有串倍长到 \(L\) 的长度，然后就可以直接计算 \(1 \sim L\) 内的答案。
所有答案只要乘一个 \(\frac{50!}{L}\) 的系数即可。
部分分 2
我们考虑一下，对于一个 \(x\)，决定最终这一天发生的所有自然灾害数量实际上是所有的 \(x \bmod p(1 \leq p \leq 50)\) 的值。
如果我们反过来，去考虑每个 \(x \bmod l_i\) 的值，如果我们钦定了 \(x \equiv a_i \pmod{l_i}\)，那么就相当于列出一个同余方程组，并且这个同余方程组在模 \(L\) 意义下最多只有一个解。有解的充要条件是：
因为部分分里所有字符串长度两两互质或相等，互质就忽略上面的条件，相等就直接枚举。
如果忽略了上面的条件，那么只要设计一个 DP，\(f(i,j)\) 表示在 \(1 \sim \mathrm{lcm}(l_1,l_2,\dots ,l_i)\) 天中，前 \(i\) 个灾害，一共发生了 \(j\) 次的天数。
接下来的思路
像上面那样直接 DP 相当于没有限制下面这个方程一定有解。
考虑如何限制有解的条件。有一个思路是，取一个参数 \(K\)，令 \(t_i=\frac{l_i}{\mathrm{gcd}(l_i,K)}\)，使得不同的 \(t_i\) 两两互质。相当于加一个同余方程，使得方程组变成：
这样的目的是枚举第一个方程之后，后面的方程之间互相就不会矛盾了。这样我们就能通过将第一个方程与其他每个方程联立的方式，知道具体在原来的字符串是哪个位置，从而计算 DP 时需要的贡献：
于是现在 \(K\) 的选取就是关键的问题。再次注意 \(K\) 需要满足的条件：令 \(t_i=\frac{l_i}{\mathrm{gcd}(l_i,K)}\)，不同的 \(t_i\) 需要两两互质。
算法一：根号分类的思想
遇到这种小范围的数论题，可以考虑一下按根号将质因子分类的经典思路（套路）：[NOI2015]寿司晚宴
大于 \(\sqrt {50}\) 的质因子在一个 \(l_i\) 最多只会出现一次。不妨取 \(K=2^5\times 3^3 \times 5^2\times 7^2\) 先试试。
求出每个 \(l_i\) 唯一大于根号的质因子（没有就看成 \(1\)），即 \(t_i=\frac{l_i}{\mathrm{gcd}(l_i,K)}\)，然后根据 \(t_i\) 将所有 \(l_i\) 分类，这样不同的 \(t_i\) 显然是互质的。用和部分分 2 一样的 DP，\(t_i\) 相同的一起转移即可。最后同样要乘上一个系数
复杂度大概是 \(\mathcal O\left(K(m\cdot n^2+\sum t_i)\right)\)，其中 \(m\) 表示不同的 \(t_i\) 的个数。发现这样不太过得去。
我们考虑把 \(K\) 取成 \(K=2^5\times 3^3 \times 5^2\)，然后强制把 \(7|t_i\) 的 \(t_i\) 都设成 \(49\)（可以认为是将满足 \(7|t_i, t_i\neq 49\) 的那些 \(7\) 倍长），这样就把时间复杂度降低了。
实现的时候有一些技巧，可以避免使用中国剩余定理，可以看代码。
#include <bits/stdc++.h>
template <class T>
inline void read(T &x)
{
static char ch;
while (!isdigit(ch = getchar()));
x = ch - '0';
while (isdigit(ch = getchar()))
x = x * 10 + ch - '0';
}
template <class T>
inline void putint(T x)
{
static char buf[25], *tail = buf;
if (!x)
putchar('0');
else
{
for (; x; x /= 10) *++tail = x % 10 + '0';
for (; tail != buf; --tail) putchar(*tail);
}
}
const int MaxL = 55;
const int MaxN = 32;
const int mod = 1e4 + 7; //prime
const int M = 32 * 27 * 25; //Solution 中说的参数 K
int n, L = 50;
bool occ[MaxL];
int cnt[MaxL][MaxL];
std::vector<int> vec[MaxL];
int ans[MaxN];
int f[MaxN], g[MaxN];
bool vis[MaxL];
inline void add(int &x, const int &y)
{
x += y;
if (x >= mod)
x -= mod;
}
int main()
{
freopen("disaster.in", "r", stdin);
freopen("disaster.out", "w", stdout);
scanf("%d", &n);
for (int i = 1; i <= n; ++i)
{
static char s[MaxL];
scanf("%s", s);
int len = strlen(s);
occ[len] = true;
for (int j = 0; j < len; ++j)
if (s[j] == 'Y')
++cnt[len][j];
}
for (int i = 1; i <= L; ++i)
if (occ[i])
{
int t = i % 7 == 0 ? 49 : i / std::__gcd(i, M);
vec[t].push_back(i);
vis[t] = true;
}
vis[32] = vis[27] = vis[25] = true;
for (int tM = 0; tM < M; ++tM)
{
memset(f, 0, sizeof(f));
f[0] = 1;
for (int P = 1; P <= L; ++P)
if (!vec[P].empty())
{
static int tf[MaxN], cur[MaxN];
memset(tf, 0, sizeof(tf));
memset(g, 0, sizeof(g));
int im = vec[P].size();
for (int i = 0; i < im; ++i)
cur[i] = tM % vec[P][i];
for (int tim = 0; tim < P; ++tim)
{
int count = 0;
for (int i = 0; i < im; ++i)
{
int len = vec[P][i];
count += cnt[len][cur[i]];
cur[i] = (cur[i] + M) % len; //注意这里的处理
//备注：这里由于 K 取得不好，实现得不够精细，实际的复杂度不是上面所说的那样，比较好的实现可以看算法二
}
++g[count];
}
for (int i = 0; i <= n; ++i)
if (f[i])
{
for (int j = 0; i + j <= n; ++j)
add(tf[i + j], 1LL * f[i] * g[j] % mod);
}
for (int i = 0; i <= n; ++i)
f[i] = tf[i];
}
for (int i = 0; i <= n; ++i)
add(ans[i], f[i]);
}
for (int i = 0; i <= n; ++i)
{
for (int j = 1; j <= L; ++j)
if (!vis[j])
ans[i] = 1LL * ans[i] * j % mod;
putint(ans[i]);
putchar('\n');
}
return 0;
}
算法二：省去冗余枚举
我们发现算法一的 \(K\) 取的并不太好，因为我们发现，并没有必要让 \(t_i\) 一定是大于根号的质因子，如果是不超过根号的也是可以的，但是要满足不同的 \(t_i\) 互质的条件。
形式化地，我们考虑每个质因子 \(p_i\) 在 \(1\sim 50\) 范围内的最大可能的次数 \(q_i\)。取 \(K=\prod p_i^{q_i-1}=2^4\times 3^2\times 5\times 7\)，不难发现这个 \(K\) 取得恰到好处。
学习这种做法的时候顺便学到了一种简化的实现方法，我们将第 \(i\) 个串的长度倍长到 \(t_i\times K\)，即 \(\mathrm{lcm}(K,l_i)\)，并且记一个数组 \(cnt_{i,j}\) 表示倍长后所有 \(t_k=i\) 的串的第 \(j\) 位的总和。那么对相同的 \(t_i\) 转移的时候，直接利用对应位置的 \(cnt\) 转移即可。
时间复杂度是 \(\mathcal O\left(K(m\cdot n^2+\sum t_i)\right)\)。
//orz pen
//What a magical solution !
#include <bits/stdc++.h>
const int K = 16 * 9 * 5 * 7;
const int MaxN = K * 50;
const int MaxM = 20;
const int mod = 1e4 + 7;
int n, m;
char s[MaxN];
bool vis[MaxN];
int l[MaxN], pos[MaxN];
int ans[MaxN];
int f[MaxM][32];
int cnt[MaxM][MaxN], g[32];
inline void add(int &x, const int &y)
{
x += y;
if (x >= mod)
x -= mod;
}
inline int get_pos(int x)
{
if (!vis[x])
vis[l[pos[x] = ++m] = x] = true;
return pos[x];
}
inline int qpow(int x, int y)
{
int res = 1;
for (; y; y >>= 1, x = 1LL * x * x % mod)
if (y & 1)
res = 1LL * res * x % mod;
return res;
}
int main()
{
freopen("disaster.in", "r", stdin);
freopen("disaster.out", "w", stdout);
scanf("%d", &n);
for (int i = 1; i <= n; ++i)
{
scanf("%s", s);
int len = strlen(s);
int x = K * len / std::__gcd(K, len);
int pos = get_pos(x);
for (int j = 0; j < x; ++j)
cnt[pos][j] += s[j % len] == 'Y';
}
for (int k = 0; k < K; ++k)
{
f[0][0] = 1;
for (int i = 1; i <= m; ++i)
{
for (int j = 0; j <= n; ++j)
g[j] = 0;
for (int j = k; j < l[i]; j += K)
++g[cnt[i][j]];
for (int j = 0; j <= n; ++j)
{
f[i][j] = 0;
for (int p = 0; p <= j; ++p)
add(f[i][j], 1LL * f[i - 1][j - p] * g[p] % mod);
}
}
for (int i = 0; i <= n; ++i)
add(ans[i], f[m][i]);
}
int cur = 1;
for (int i = 1; i <= 50; ++i)
cur = 1LL * cur * i % mod;
for (int i = 1; i <= m; ++i)
cur = 1LL * cur * qpow(l[i] / K, mod - 2) % mod;
cur = 1LL * cur * qpow(K, mod - 2) % mod;
for (int i = 0; i <= n; ++i)
{
ans[i] = 1LL * ans[i] * cur % mod;
printf("%d\n", ans[i]);
}
return 0;
}
