UOJ #54. 【WC2014】时空穿梭
题目链接
题目大意
一个 \(n\) 维空间，第 \(i\) 维的值域为 \([1,m_i]\)，从空间中选出 \(c\) 个整点，满足各维坐标分别递增，且 \(c\) 个点共线。求方案数。
多测，\(T\leq 100\)，\(n\leq 11,c\leq 20,m_i\leq 100000\)
思路
枚举各维的极差 \(d_i\)，则对答案的贡献为 \(\displaystyle \prod_i (m_i-d_i)\binom{\gcd(d_1,..,d_n)-1}{c-2}\) 。
\[\begin{align}
&\sum_{d_1,d_2,...,d_n}\prod_{i=1}^n (m_i-d_i)\binom{\gcd(d_1,..,d_n)-1}{c-2} \\
=&\sum_{g}\sum_{\gcd(d_1,d_2,...,d_n)=1}\prod_{i=1}^n (m_i-gd_i)\binom{g-1}{c-2} \\
=&\sum_{g}\sum_{d|\gcd(d_1,d_2,...,d_n)}\mu(d)\prod_{i=1}^n(m_i-gd_i)\binom{g-1}{c-2}\\
=&\sum_g \binom{g-1}{c-2}\sum_d\mu(d)\sum_{d_1,...,d_n}\prod_{i=1}^n (m_i-gd\cdot d_i)\\
=&\sum_{g}\binom{g-1}{c-2}\sum_{d}\mu(d)\prod_{i=1}^n\frac{(2m_i-(\lfloor\frac{m_i}{gd}\rfloor+1) gd)\lfloor\frac{m_i}{gd}\rfloor}{2} \\
=&\sum_{g}\left(\sum_{d|g}\binom{d-1}{c-2}\mu(\frac{g}{d}) \right)\prod_{i=1}^n \frac{(2m_i-(\lfloor\frac{m_i}{g}\rfloor+1) g)\lfloor\frac{m_i}{g}\rfloor}{2}
\end{align}
\]
记 \(\displaystyle f(g)=\sum_{d|g}\binom{d-1}{c-2} \mu(\frac{g}{d})\)，\(f(g)\) 可以通过枚举 \(\displaystyle \mu(\frac{g}{d})\neq 0\) 的 \(\displaystyle\frac{g}{d}\) 然后调和级数求出，对于每一个 \(c\) 都要预处理一遍，从而 \(O(cm\log m)\) 。
djq：对于任意 \(n\)，\(n\) 以内都大概有 \(\displaystyle \frac{6n}{\pi^2}\) 的 \(\mu\neq 0\) 的数。
到这里就可以 \(O(Tnm)\) 算答案了，但是极度卡常。
注意到 \(n,m\) 规模差异悬殊，考虑时间复杂度向 \(n\) 倾斜，于是整除分块，设 \(k_i\) 为一段中的 \(\lfloor\frac{m_i}{g}\rfloor\)：
\[\begin{align}
=&\sum_{g=l}^r f(g)\prod_{i=1}^n \frac{k_i}{2} (2m_i-(k_i+1)g)\\
=&\sum_{g=l}^r f(g)\sum_{i=0}^n a_i g^i \\
=&\sum_{i=0}^n a_i\sum_{g=l}^r f(g)g^i \\
=& \sum_{i=0}^n a_i (S_i(r)-S_i(l-1))
\end{align}
\]
\(\displaystyle \sum_{i=0}^n a_i g^i\) 表示积式 \(\displaystyle \prod_{i=1}^n \frac{k_i}{2} (2m_i-(k_i+1)g)\) 形成的多项式，可以在整除分块时动态维护，每次 \(k_i\) 变化时 \(O(n)\) 乘上或除去一个 \(ax+b\) （遇到除 \(0\) 的情况需要精细处理）。
\(\displaystyle S_i(n)=\sum_{x=1}^n f(x)x^i\)，\(O(cnm)\) 预处理即可。
整除分块共 \(O(n\sqrt m)\) 段，每次需要 \(O(n)\) 维护多项式以及 \(O(n)\) 计算答案，于是询问复杂度 \(O(Tn^2\sqrt m)\) 。
综上，时间复杂度 \(O(cm(n+\log m)+Tn^2\sqrt m)\) 。
Code
\(O(Tnm)\) 的大力卡常代码，循环展开真好用
#include<iostream>
#define rep(i,a,b) for(int i = (a); i <= (b); i++)
#define per(i,b,a) for(int i = (b); i >= (a); i--)
#define mod 10007
#define inv2 5004
#define V 101000
#define N 12
#define M 21
#define ll long long
using namespace std;
int _C[V][M];
int mu[V], prime[V], siz, coef[V][M], inv[M];
bool vis[V];
ll n, c, m[N];
ll a[N], b[N];
void init(){
_C[0][0] = 1;
rep(i,1,V-1){
_C[i][0] = 1;
rep(j,1,min(M-1,i)) _C[i][j] = (_C[i-1][j-1] + _C[i-1][j]) % mod;
}
}
int C(int n, int m){ return m < 0 || n < m ? 0 : _C[n][m]; }
void euler(){
mu[1] = 1;
rep(i,2,V-1){
if(!vis[i]) mu[i] = -1, prime[++siz] = i;
for(int j = 1; j <= siz && prime[j] * i < V; j++){
int v = i * prime[j];
vis[v] = true, mu[v] = mu[i] * -1;
if(i%prime[j] == 0){ mu[v] = 0; break; }
}
}
}
void prework(){
rep(i,1,V-1) for(int j = 1; j*j <= i; j++) if(i%j == 0){
if(mu[i/j]) rep(c,1,M-1) (coef[i][c] += C(j-1, c-2) * mu[i/j]) %= mod;
if(j*j != i && mu[j]) rep(c,1,M-1) (coef[i][c] += C(i/j-1, c-2) * mu[j]) %= mod;
}
rep(i,1,V-1) rep(c,1,M-1) if(coef[i][c] < 0)
coef[i][c] += mod;
}
int main(){
init(), euler(), prework();
int T; cin>>T;
while(T--){
cin>>n>>c;
ll mn = V;
rep(i,1,n) cin>>m[i], mn = min(mn, m[i]);
if(n == 1){ cout<< C(mn, c) <<endl; continue; }
if(n == 2 && c == 3){
ll ans = 0;
for(int g = 1; g+7 <= mn; g += 8){
if(m[1]/g == m[1]/(g+7) && m[2]/g == m[2]/(g+7)){
ll a = (m[1]/g), b = (m[2]/g), c = m[1]%g, d = m[2]%g;
(ans += a * b % mod * (
coef[g][3] * (m[1]-g + c) * (m[2]-g + d) % mod +
coef[g+1][3] * (m[1]-g-1 + c-a) * (m[2]-g-1 + d-b) % mod +
coef[g+2][3] * (m[1]-g-2 + c-2*a) * (m[2]-g-2 + d-2*b) % mod +
coef[g+3][3] * (m[1]-g-3 + c-3*a) * (m[2]-g-3 + d-3*b) % mod +
coef[g+4][3] * (m[1]-g-4 + c-4*a) * (m[2]-g-4 + d-4*b) % mod +
coef[g+5][3] * (m[1]-g-5 + c-5*a) * (m[2]-g-5 + d-5*b) % mod +
coef[g+6][3] * (m[1]-g-6 + c-6*a) * (m[2]-g-6 + d-6*b) % mod +
coef[g+7][3] * (m[1]-g-7 + c-7*a) * (m[2]-g-7 + d-7*b) % mod
)) %= mod;
} else rep(i,g,g+7)
(ans += coef[i][3] * (m[1]-i + m[1]%i) * (m[1]/i) % mod * (m[2]-i + m[2]%i) * (m[2]/i) % mod) %= mod;
}
rep(g,(mn/8)*8+1,mn) (ans += coef[g][3] * (m[1]-g + m[1]%g) * (m[1]/g) % mod
* (m[2]-g + m[2]%g) * (m[2]/g) % mod) %= mod;
cout<< ans * inv2 * inv2 % mod <<endl;
continue;
}
int ans = 0;
for(int g = 1; g+3 <= mn; g += 4){
bool same = true;
rep(i,1,n) same &= m[i]/g == m[i]/(g+3), a[i] = m[i]/g, b[i] = m[i]%g;
if(same) rep(p,0,3){
ll val = coef[g+p][c];
rep(i,1,n) (val *= (m[i]-g-p + b[i]-p*a[i]) * a[i]) %= mod;
(ans += val) %= mod;
} else rep(p,g,g+3){
ll val = coef[p][c];
for(int i = 1; i+4 <= n; i += 5){
(val *= (m[i]-p + m[i]%p) * (m[i]/p) % mod) %= mod;
(val *= (m[i+1]-p + m[i+1]%p) * (m[i+1]/p) % mod) %= mod;
(val *= (m[i+2]-p + m[i+2]%p) * (m[i+2]/p) % mod) %= mod;
(val *= (m[i+3]-p + m[i+3]%p) * (m[i+3]/p) % mod) %= mod;
(val *= (m[i+4]-p + m[i+4]%p) * (m[i+4]/p) % mod) %= mod;
}
rep(i,(n/5)*5+1,n) (val *= (m[i]-p + m[i]%p) * (m[i]/p) % mod) %= mod;
(ans += val) %= mod;
}
}
rep(g,(mn/4)*4+1,mn){
ll val = coef[g][c];
rep(i,1,n) (val *= (m[i]-g + m[i]%g) * (m[i]/g) % mod) %= mod;
(ans += val) %= mod;
}
rep(i,1,n) (ans *= inv2) %= mod;
cout<< ans <<endl;
}
return 0;
}
\(O(Tn^2\sqrt m)\) ：
#include<iostream>
#include<fstream>
#define rep(i,a,b) for(int i = (a); i <= (b); i++)
#define per(i,b,a) for(int i = (b); i >= (a); i--)
#define mod 10007
#define N 12
#define M 21
#define V 101000
using namespace std;
int n, c, m[N];
int C[V][M], mu[V], inv[V], f[V][M], S[M][V][N];
int prime[V], siz;
bool vis[V];
void prework(){
inv[0] = inv[1] = 1;
rep(i,2,V-1) inv[i] = (mod-mod/i) * inv[mod%i] % mod;
C[0][0] = 1;
rep(i,1,V-1){
C[i][0] = 1;
rep(j,1,min(i,M-1)) C[i][j] = (C[i-1][j-1] + C[i-1][j]) % mod;
}
mu[1] = 1;
rep(i,2,V-1){
if(!vis[i]) mu[i] = -1, prime[++siz] = i;
for(int j = 1; j <= siz && i * prime[j] < V; j++){
int v = i * prime[j];
vis[v] = true, mu[v] = mu[i] * -1;
if(i % prime[j] == 0){ mu[v] = 0; break; }
}
}
rep(c,2,M-1) rep(d,1,V-1) if(mu[d])
for(int g = d; g < V; g += d)
(f[g][c] += C[g/d-1][c-2] * mu[d]) %= mod;
rep(c,2,M-1) rep(g,1,V-1){
if(f[g][c] < 0) f[g][c] += mod;
int pow = 1;
rep(i,0,N-1) S[c][g][i] = (S[c][g-1][i] + f[g][c] * pow) % mod, (pow *= g) %= mod;
}
}
int a[N], zero;
void mul(int x, int y){
x %= mod;
if(x == 0 && y == 0){ zero++; return; }
per(i,n,0){
if(i < n) (a[i+1] += a[i] * x) %= mod;
(a[i] *= y) %= mod;
}
}
void dinv(int x, int y){
x %= mod;
if(y == 0){
if(x == 0){ zero--; return; }
rep(i,0,n-1) a[i] = a[i+1]; a[n] = 0;
rep(i,0,n) (a[i] *= inv[x]) %= mod;
return;
}
rep(i,0,n){
(a[i] *= inv[y]) %= mod;
if(i < n) (a[i+1] += mod - a[i] * x % mod) %= mod;
}
}
int main(){
prework();
int T; cin>>T;
while(T--){
cin>>n>>c;
int mn = V;
rep(i,1,n) cin>>m[i], mn = min(mn, m[i]);
int ans = 0;
rep(i,0,n) a[i] = 0; a[0] = 1;
for(int l = 1, r; l <= mn; l = r+1){
r = mn;
rep(i,1,n){
r = min(r, m[i]/(m[i]/l));
if(l == 1 || m[i]/l != m[i]/(l-1)){
if(l > 1){
int k = m[i]/(l-1) % mod;
dinv(mod - k*(k+1) % mod * inv[2] % mod, k * m[i] % mod);
}
int k = m[i]/l % mod;
mul(mod - k*(k+1) % mod * inv[2] % mod, k * m[i] % mod);
}
}
rep(i,0,n) (ans += (zero ? 0 : a[i]) * (S[c][r][i] + mod - S[c][l-1][i])) %= mod;
}
cout<< ans <<endl;
}
return 0;
}
