uoj54-bzoj3434-时空穿梭

题意

在一个 (n) 维空间中，求一个点可以用一个 (n) 维向量 ((x_1,x_2,dots x_n)) 表示。现在要选出 (c) 个点，有三个限制：

• 设 (x_i) 表示任意一个点的第 (i) 个分量，那么 (x_iin [1,m_i],x_iin Z) 。
• 选出的所有点在同一条直线上。

问方案数。 (Tle 100,nle 11,cle 20,m_ile 10^5) 。

分析

好题！！！

看到这种感觉比较难入手的统计方案数问题，我们可以先尝试枚举一些东西，把问题变简单再计算。如果是很复杂的问题（比如说这里的 (n) 维），可以先从简单的情况开始推理。这里我们从 (n=2) 开始想，第一步枚举这条线的两端，或者枚举每个分量上的长度再枚举起点：

[sum _{x=1}^{m_1}sum _{y=1}^{m_2}sum _{i=1}^{m_1-x}sum _{j=1}^{m_2-y}inom {gcd(x,y)-1}{c-2} ]

即枚举每个分量上的长度，那么这个线段上的整点个数就是 (gcd(x,y)-1) ，分别为 ((frac{x}{gcd(x,y)},frac{y}{gcd(x,y)}),(frac{2x}{gcd(x,y)},frac{2y}{gcd(x,y)}),dots ,(frac{(gcd(x,y)-1)x}{gcd(x,y)},frac{(gcd(x,y)-1)y}{gcd(x,y)})) ，统计取出 (c-2) 个的方案数（ ((0,0)) 和 ((x,y)) 已经取了）。用经典的莫比乌斯函数方法对这个式子实施一些变形。设 (m_1le m_i) ：

[egin{aligned} ans&=sum _{d=1}^{m_1}inom {d-1}{c-2} sum _{x=1}^{lfloorfrac{m_1}{d} floor}sum _{y=1}^{lfloorfrac{m_2}{d} floor}[gcd(x,y)=1](m_1-dx)(m_2-dy) \ &=sum _{d=1}^{m_1}inom {d-1}{c-2} sum _{x=1}^{lfloorfrac{m_1}{d} floor}sum _{y=1}^{lfloorfrac{m_2}{d} floor}(m_1-dx)(m_2-dy)sum _{e|x,e|y}mu (e) \ &=sum _{d=1}^{m_1}inom {d-1}{c-2}sum _{e=1}^{lfloorfrac{m_1}{d} floor}mu (e)sum _{x=1}^{lfloorfrac{m_1}{de} floor}sum _{y=1}^{lfloorfrac{m_2}{de} floor} (m_1-xde)(m_2-yde) end{aligned} ]

令 (f(m,k)=sum _{i=1}^{lfloorfrac{m}{k} floor}(m-ik)) ，那么就有

[egin{aligned} ans&=sum _{d=1}^{m_1}inom {d-1}{c-2}sum _{e=1}^{lfloorfrac{m_1}{d} floor}mu (e)prod _{i=1}^nf(m_i,de) \ &=sum _{a=1}^{m_1}(prod _{i=1}^nf(m_i,a))sum _{d|a}mu (frac{a}{d})inom {d-1}{c-2} end{aligned} ]

注意，这里我们已经把结论推广到 (n) 维的情况了——找到每一维独立的东西，按同样的方式加入。再看看 (f) 的形式：

[egin{aligned} f(m,k)&=sum _{i=1}^{lfloorfrac{m}{k} floor}(m-ik) \ &=mlfloorfrac{m}{k} floor-kfrac{lfloorfrac{m}{k} floor(lfloorfrac{m}{k} floor+1)}{2} end{aligned} ]

这里面带有 (k) 的底式很难处理。而注意到数据范围比较小，我们可以直接对于一种 (lfloorfrac{m}{k} floor) 把 (f(m,k)) 多项式直接乘出来。设得到的多项式系数向量为 (b) ，那么：

[egin{aligned} ans&=sum _{a=down}^{up}sum _{j=0}^nb_ja^jsum _{d|a}mu (frac{a}{d})inom {d-1}{c-2} end{aligned} ]

如果我们预先处理出 (g(a,c,k)=a^ksum _{d|a}mu (frac{a}{d})inom {d-1}{c-2}) 的前缀和，那么我们就可以分段快速计算答案了。

接下来分析复杂度。预处理的复杂度为 (O(mcln m+nmc)) ，前一个为枚举倍数的调和级数计算 (g(a,c,0)) ，后面的是递推得到 (g(a,c,k)) 的前缀和。

整个查询会被分成最多 (2nsqrt m) 段，每段中需要暴力乘出多项式，这个复杂度为 (n^2) ，所以总复杂度为 (O(mcln m+nmc+n^3sqrt m)) 。可能我的常数比较大，所以bzoj上过不了。uoj总用时4秒。

有一个需要注意的地方，由于这题的模数是 (10^4+7<m) ，我们求组合数就不能用阶乘来求，由于这题的数据范围比较小，可以直接递推出所有的。在线性筛 (mu) 函数的时候记得对 (mu) 取模。

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long giant;
inline char nchar() {
	static const int bufl=1<<20;
	static char buf[bufl],*a,*b;
	return a==b && (b=(a=buf)+fread(buf,1,bufl,stdin),a==b)?EOF:*a++;
}
inline int read() {
	int x=0,f=1;
	char c=nchar();
	for (;!isdigit(c);c=nchar()) if (c=='-') f=-1;
	for (;isdigit(c);c=nchar()) x=x*10+c-'0';
	return x*f;
}
const int maxn=12;
const int maxc=21;
const int maxm=1e5+1;
const int q=1e4+7;
const int inf=1e9+7;
inline int Plus(int x,int y) {return ((giant)x+(giant)y)%q;}
inline int Sub(int x,int y) {return Plus(x,q-y);}
inline int Multi(int x,int y) {return (giant)x*y%q;}
inline void Max(int &x,int y) {x=max(x,y);}
inline void Min(int &x,int y) {x=min(x,y);}
bool np[maxm];
int p[maxm],ps=0,mu[maxm],g[maxc][maxn][maxm],m[maxn],n,c,mi,com[maxm][maxc];
inline int C(int n,int m) {return n>=m?com[n][m]:0;}
struct poly {
	int a[maxn],e;
	inline int& operator [] (const int x) {return a[x];}
	void clear() {memset(this,0,sizeof(poly));}
	poly () {clear();}
} b[maxn];
poly operator * (poly a,poly b) {
	poly c;
	for (int i=0;i<=a.e;++i) for (int j=0;j<=b.e;++j) c[i+j]=Plus(c[i+j],Multi(a[i],b[j]));
	c.e=a.e+b.e;
	return c;
}
inline int solve() {
	if (mi<n) return 0;
	int ret=0;
	for (int i=1,j;i<=mi;i=j+1) {
		j=inf;
		for (int k=1;k<=n;++k) {
			int tmp=m[k]/i;
			Min(j,m[k]/(tmp));
			b[k].clear();
			b[k].e=1;
			b[k][0]=Multi(m[k],tmp);
			b[k][1]=q-((giant)tmp*(tmp+1)/2ll%q);
		}
		for (int k=2;k<=n;++k) b[1]=b[1]*b[k];
		for (int k=0;k<=n;++k) ret=Plus(ret,Multi(b[1][k],Sub(g[c][k][j],g[c][k][i-1])));
	}
	return ret;
}
int main() {
#ifndef ONLINE_JUDGE
	freopen("test.in","r",stdin);
#endif
	com[0][0]=1;
	for (int i=1;i<maxm;++i) {
		com[i][0]=1;
		for (int j=1;j<maxc;++j) com[i][j]=Plus(com[i-1][j],com[i-1][j-1]);
	}
	mu[1]=1;
	for (int i=2;i<maxm;++i) {
		if (!np[i]) p[++ps]=i,mu[i]=q-1; 
		for (int j=1,tmp;j<=ps && (tmp=i*p[j])<maxm;++j) {
			np[tmp]=true;
			if (i%p[j]==0) break;
			mu[tmp]=q-mu[i];
		}
	}
	for (int x=2;x<maxc;++x) {
		int (*h)[maxm]=g[x];
		for (int i=1;i<maxm;++i) for (int j=1,tmp;(tmp=i*j)<maxm;++j) h[0][tmp]=Plus(h[0][tmp],Multi(mu[j],C(i-1,x-2)));
		for (int i=1;i<maxn;++i) for (int j=1;j<maxm;++j) h[i][j]=Multi(h[i-1][j],j);
		for (int i=0;i<maxn;++i) for (int j=1;j<maxm;++j) h[i][j]=Plus(h[i][j-1],h[i][j]);
	}
	int T=read();
	while (T--) {
		n=read(),c=read(),mi=inf;
		for (int i=1;i<=n;++i) Min(mi,m[i]=read());
		int ans=solve();
		printf("%d
",ans);
	}
	return 0;
}