BZOJ3434 [WC2014] 时空穿梭
Description
有一个$n$维的超立方体$[1, m_1]\times[1, m_2]\times\dots\times[1,m_n]$。
现在要再其内部选择$c$个整点，使得其每一维都是递增的（也即，$\forall i=1\dots c-1, j=1\dots n, x_{i+1, j}>x_{i,j}$），且这$c$个点共线。求方案数。
多组数据，至多$1000$组；$n\leqslant 11, c\leqslant 20, m_i\leqslant 10^5$。
Solution
首先，我们发现只要满足$x_i < y_i, i = 1,2,\dots,n$ ，（其中$(x_1, x_2,\dots,x_n),(y_1, y_2,\dots,y_n)$分别是第一个点和最后一个点的坐标）且点的不同排列算同一种即可。
那么，我们枚举这两个点之后，以其为端点的线段上共有（除去两端）$\gcd(y_1-x_1,y_2-x_2,\dots.y_n-x_n)$个点。其中选出$c-2$个点的方式有$C_{\gcd(y_1-x_1,y_2-x_2,\dots.y_n-x_n)}^{c-2}$种，所以总方案数是
$$ans=\sum_{\substack{1\leqslant x_i < y_i\leqslant m_i \\ i = 1,2,\dots,n}}C_{\gcd(y_1-x_1,y_2-x_2,\dots,y_n-x_n)}^{c-2}$$
我们发现，所有$p_i=y_i-x_i$固定之后，方案数恰有
$$C_{\gcd(p_1, p_2,\dots,p_n)}^{c-2}\prod_{i=1}^n(m_i-p_i)$$
种，这是因为$p_i$固定后$x_i$只能取$1,2,\dots,m_i-p_i$.
于是我们有
$$ans=\sum_{\substack{1\leqslant q_i < m_i \\ i = 1,2,\dots,n}}C_{\gcd(p_1, p_2,\dots,p_n)}^{c-2}\prod_{i=1}^n(m_i-p_i)$$
优先枚举$d=\gcd(p_1, p_2,\dots,p_n)$（下式中$m_{min}$表示$\min(m_1,m_2,\dots,m_n)$:
$$ans=\sum_{d=1}^{m_{min}-1}C_d^{c-2}\sum_{\substack{1\leqslant q_i < m_i \\ i = 1,2,\dots,n\\ \gcd(p_1, p_2,\dots,p_n)=d}}\prod_{i=1}^n(m_i-p_i)$$
令
$$f(d)=\sum_{\substack{1\leqslant q_i < m_i \\ i = 1,2,\dots,n\\ \gcd(p_1, p_2,\dots,p_n)=d}}\prod_{i=1}^n(m_i-p_i)$$
则有
$$\begin{aligned} \sum_{d\mid l}f(l)&=\sum_{\substack{1\leqslant q_i < m_i \\ i = 1,2,\dots,n\\d\mid\gcd(p_1, p_2,\dots,p_n)}}\prod_{i=1}^n(m_i-p_i)\\ &=\sum_{\substack{1\leqslant q_i'd < m_i \\ i = 1,2,\dots,n}}\prod_{i=1}^n(m_i-p_i'd)\\ &=\prod_{i=1}^n\sum_{p'=1}^{\left\lfloor\frac{m_i-1}d\right\rfloor}(m_i-p'd)\\ &=\prod_{i=1}^n\left[m_i\left\lfloor\frac{m_i-1}d\right\rfloor-\frac{d\left\lfloor\frac{m_i-1}d\right\rfloor\left(1+\left\lfloor\frac{m_i-1}d\right\rfloor\right)}2\right] \end{aligned}$$
最后一步是等差数列求和。
由莫比乌斯反演可知：
$$f(d)=\sum_{d\mid l}\mu\left(\frac ld \right)\prod_{i=1}^n\left[m_i\left\lfloor\frac{m_i-1}l\right\rfloor-\frac{l\left\lfloor\frac{m_i-1}l\right\rfloor\left(1+\left\lfloor\frac{m_i-1}l\right\rfloor\right)}2\right]$$
于是有
$$\begin{aligned} ans&=\sum_{d=1}^{m_{min}-1}C_d^{c-2}\sum_{d\mid l}\mu\left(\frac ld \right)\prod_{i=1}^n\left[m_i\left\lfloor\frac{m_i-1}l\right\rfloor-\frac{l\left\lfloor\frac{m_i-1}l\right\rfloor\left(1+\left\lfloor\frac{m_i-1}l\right\rfloor\right)}2\right]\\ &=\sum_{l=1}^{m_{min}-1}\prod_{i=1}^n\left[m_i\left\lfloor\frac{m_i-1}l\right\rfloor-\frac{l\left\lfloor\frac{m_i-1}l\right\rfloor\left(1+\left\lfloor\frac{m_i-1}l\right\rfloor\right)}2\right]\sum_{d\mid l}\mu\left(\frac ld \right)C_d^{c-2} \end{aligned}$$
易知$\prod_{i=1}^n\left[m_i\left\lfloor\frac{m_i-1}l\right\rfloor-\frac{l\left\lfloor\frac{m_i-1}l\right\rfloor\left(1+\left\lfloor\frac{m_i-1}l\right\rfloor\right)}2\right]$在将所有$\left\lfloor\frac{m_i-1}l\right\rfloor$作为常数看待的情况下是一个关于$l$的不超过$n$次（是因为$\bmod 10007$之后有可能次数更低）的多项式，而且引起至少一个$\left\lfloor\frac{m_i-1}l\right\rfloor$改变的$l$有$O(n\sqrt m)$个，于是可以对每一个所有$\left\lfloor\frac{m_i-1}l\right\rfloor$都不改变的段分别求和。由于此时它是关于$l$的多项式，我们只需预处理出所有的$S_{c,p,t}=\sum_{l=1}^tl^p\sum_{d\mid l}\mu\left(\frac ld \right)C_d^{c-2}$即可。至于如何在只改变一个因式时快速维护乘积中各项的系数，可以用线段树实现（代码中为zkw线段树）。
完。
Code
#include <algorithm>
#include <cstdio>
const int N = 12;
const int C = 21;
const int M = 100001;
const int mod = 10007;
int fac[mod], inv[mod];
int mu[M], f[C][M], S[C][M][N];
int p[N * 1000], m[N];
inline int CC(int a, int b) {
if (!b) return 1;
if (a % mod < b % mod) return 0;
return CC(a / mod, b / mod) * fac[a % mod] % mod * inv[fac[b % mod]] % mod * inv[fac[(a - b) % mod]] % mod;
}
int a[64][N], po[64];
int mdp[N];
int n, mm, c;
void upd(int x) {
int l = x * 2, r = x * 2 + 1;
po[x] = po[l] + po[r];
for (int i = 0; i <= po[x]; ++i) a[x][i] = 0;
for (int i = 0; i <= po[l]; ++i)
for (int j = 0; j <= po[r]; ++j)
a[x][i + j] = (a[x][i + j] + a[l][i] * a[r][j] % mod) % mod;
}
int sum(int l, int r) {
int ans = 0;
for (int i = 0; i <= po[1]; ++i)
ans = (ans + a[1][i] * (S[c][r][i] - S[c][l - 1][i]) % mod) % mod;
return ans;
}
int main() {
fac[0] = 1;
for (int i = 1; i < mod; ++i) fac[i] = fac[i - 1] * i % mod;
inv[1] = 1;
for (int i = 2; i < mod; ++i) inv[i] = mod - (mod / i) * inv[mod % i] % mod;
mu[1] = 1;
for (int i = 1; i < M; ++i)
for (int j = i * 2; j < M; j += i)
mu[j] -= mu[i];
for (int c = 2; c < C; ++c) {
for (int i = 1; i < M; ++i)
f[c][i] = 0;
for (int d = c - 1; d < M; ++d) {
int cc = CC(d - 1, c - 2);
for (int i = 1; i * d < M; ++i)
f[c][i * d] = (f[c][i * d] + cc * mu[i]) % mod;
}
for (int j = 0; j < N; ++j)
S[c][0][j] = 0;
for (int i = 1; i < M; ++i)
for (int j = 0, ij = 1; j < N; ++j, ij = ij * (i % mod) % mod)
S[c][i][j] = (S[c][i - 1][j] + f[c][i] * ij) % mod;
}
int T;
scanf("%d", &T);
while (T--) {
mm = 10000000;
scanf("%d%d", &n, &c);
for (int i = 0; i < n; ++i)
scanf("%d", &m[i]), mm = std::min(mm, m[i]);
if (mm < c) {
printf("0\n");
continue;
} if (n == 1) {
printf("%d\n", CC(m[0], c));
continue;
}
int t = 0;
p[t++] = 1;
for (int i = 0; i < n; ++i) {
int last = 1;
while (last < m[i] - 1) {
last = (m[i] - 1) / ((m[i] - 1) / last) + 1;
p[t++] = last;
}
}
std::sort(p, p + t);
t = (int)(std::unique(p, p + t) - p);
p[t] = mm;
int d = 1;
while (d < n) d <<= 1;
for (int i = 0; i < n; ++i) mdp[i] = m[i];
for (int i = 1; i <= d * 2; ++i) a[i][po[i] = 0] = 1;
int ans = 0;
for (int i = 0; p[i] < mm && i < t; ++i) {
int l = p[i], r = p[i + 1] - 1; //[l, r]
for (int j = 0; j < n; ++j)
if ((m[j] - 1) / l != mdp[j]) {
mdp[j] = (m[j] - 1) / l;
int k = d + j;
po[k] = 1;
a[k][0] = m[j] % mod * (mdp[j] % mod) % mod;
a[k][1] = -mdp[j] % mod * ((mdp[j] + 1) % mod) % mod * inv[2] % mod;
while (k /= 2)
upd(k);
}
ans = (ans + sum(l, r)) % mod;
}
printf("%d\n", (ans + mod) % mod);
}
return 0;
}