[WC2014]时空穿梭

声明:我还没有退役,只是Typora太好用,懒得写博客了而已,比较好的题还是会放上来

鼓起勇气写了一个冬令营题

有了[BZOJ3518] 点组计数的基础

显然答案为

\[\sum_{x_1=1}^{m_1}\sum_{x_2=1}^{m_2}...\sum_{x_n=1}^{m_n}(^{gcd(x_1,...x_n)-1}_{c-2})*\prod_{i=1}^{n}(m_i-x_i) \]

然而这次\(gcd\)函数在组合数里面

并不能用\(\varphi\)反演

所以这次用\(\mu\)反演

\[\sum_{x_1=1}^{m_1}\sum_{x_2=1}^{m_2}...\sum_{x_n=1}^{m_n}(^{gcd(x_1,...x_n)-1}_{c-2})*\prod_{i=1}^{n}{(m_i-x_i)} \\=\sum_{d=1}^{m}{(^{d-1}_{c-2})}\sum_{x_1=1}^{\lfloor \frac{m_1}{d}\rfloor}\sum_{x_2=1}^{\lfloor \frac{m_2}{d}\rfloor}...\sum_{x_n=1}^{\lfloor \frac{m_n}{d}\rfloor}{[gcd(x_1,...x_n)=1]*\prod_{i=1}^{n}{(m_i-d*x_i)}} \\=\sum_{d=1}^{m}{(^{d-1}_{c-2})}\sum_{x_1=1}^{\lfloor \frac{m_1}{d}\rfloor}\sum_{x_2=1}^{\lfloor \frac{m_2}{d}\rfloor}...\sum_{x_n=1}^{\lfloor \frac{m_n}{d}\rfloor}{\sum_{t|gcd(x_1,...x_n)}{\mu(t)}\prod_{i=1}^{n}{(m_i-d*x_i)}} \\=\sum_{d=1}^{m}{(^{d-1}_{c-2})}\sum_{t=1}^{\lfloor \frac{m}{d}\rfloor}{\mu(t)}\sum_{x_1=1}^{\lfloor \frac{m_1}{dt}\rfloor}\sum_{x_2=1}^{\lfloor \frac{m_2}{dt}\rfloor}...\sum_{x_n=1}^{\lfloor \frac{m_n}{dt}\rfloor}\prod_{i=1}^{n}{(m_i-d*t*x_i)} \\=\sum_{d=1}^{m}{(^{d-1}_{c-2})}\sum_{t=1}^{\lfloor \frac{m}{d}\rfloor}{\mu(t)}\sum_{x_1=1}^{\lfloor \frac{m_1}{dt}\rfloor}{{(m_1-d*t*x_1)}}\sum_{x_2=1}^{\lfloor \frac{m_2}{dt}\rfloor}{(m_2-d*t*x_2)}...\sum_{x_n=1}^{\lfloor \frac{m_n}{dt}\rfloor}{(m_n-d*t*x_n)} \\=\sum_{d=1}^{m}{(^{d-1}_{c-2})}\sum_{t=1}^{\lfloor \frac{m}{d}\rfloor}{\mu(t)}\prod_{i=1}^{n}{(m_i*\lfloor \frac{m_i}{dt}\rfloor-d*t*\frac{\lfloor \frac{m_1}{dt}\rfloor*(\lfloor \frac{m_1}{dt}\rfloor+1)}{2})} \\=\sum_{d=1}^{m}\sum_{d'|d}{(^{d'-1}_{c-2})*\mu(\frac{d}{d'})*\prod_{i=1}^{n}{(m_i*\lfloor \frac{m_i}{d}\rfloor-d*\frac{\lfloor \frac{m_1}{d}\rfloor*(\lfloor \frac{m_1}{d}\rfloor+1)}{2})}} \]

至此

复杂度为\(O(Tmn)\)算下来大概\(1e8\)

然而还可以再优化

设

\[F(d)=\prod_{i=1}^{n}{(m_i*\lfloor\frac{m_i}{d}\rfloor-d*\frac{\lfloor\frac{m_1}{d}\rfloor*(\lfloor \frac{m_1}{d}\rfloor+1)}{2})} \]

显然\(F(d)\)是关于\(d\)的一个\(n\)次函数

可以\(O(n^2)\)计算其所有系数

设这些系数是\(a_0,a_1,...a_n\)

原式即可变为

\[\sum_{d=1}^{m}\sum_{d'|d}{(^{d'-1}_{c-2})*\mu(\frac{d}{d'})*\sum_{i=0}^{n}{a_i*d^i}} \\=\sum_{d=1}^{m}\sum_{i=0}^{n}{a_i*(d^i*\sum_{d'|d}{(^{d'-1}_{c-2})*\mu(\frac{d}{d'}}))} \]

后面那部分可以在\(O(mnc)\)内预处理出前缀和

每次询问\(O(n^3\sqrt{m})\)(处理\(a\)是\(n^2\),需要算\(n\sqrt{m}\)次,因为\(d\)有\(n\sqrt{m}\)种取值)

代码如下

#include<bits/stdc++.h>

using namespace std;

#define gc c=getchar()
#define r(x) read(x)
#define ll long long 

template<typename T>
inline void read(T&x){
    x=0;T k=1;char gc;
    while(!isdigit(c)){if(c=='-')k=-1;gc;}
    while(isdigit(c)){x=x*10+c-'0';gc;}x*=k;
}

const int p=10007;
const int N=1e5+7;
const int S=21;
const int W=12;

inline int add(int a,int b){
	a+=b;
	if(a>=p)a-=p;
	return a;
}

int C[N][S],P[N][W],f[N][S],sum[N][S][W];

int tot;
int pri[N],mu[N];
bool mark[N];

inline void pre(int n){
	mu[1]=1;
	for(int i=2;i<=n;++i){
		if(!mark[i])pri[++tot]=i,mu[i]=-1;
		for(int j=1,tmp;j<=tot&&(tmp=i*pri[j])<=n;++j){
			mark[tmp]=1;
			if(i%pri[j]==0){
				mu[tmp]=0;
				break;
			}
			mu[tmp]=-mu[i];
		}
	}
	C[0][0]=1;
	for(int i=1;i<=n;++i){
		C[i][0]=1;
		for(int j=1;j<S;++j){
			C[i][j]=add(C[i-1][j-1],C[i-1][j]);
		}
	}
	for(int i=0;i<=n;++i){
		P[i][0]=1;
		for(int j=1;j<W;++j){
			P[i][j]=P[i][j-1]*i%p;
		}
	}
	for(int i=1;i<=n;++i){
		for(int j=1;i*j<=n;++j){
			for(int k=2;k<S;++k){
				f[i*j][k]=add(f[i*j][k],add(p,C[i-1][k-2]*mu[j]));
			}
		}
	}
	for(int i=1;i<=n;++i){
		for(int j=2;j<S;++j){
			for(int k=0;k<W;++k){
				sum[i][j][k]=add(sum[i][j][k],add(sum[i-1][j][k],f[i][j]*P[i][k]%p));
			}
		}
	}
}

int a[W];

struct query{
	int n,c,M;
	int m[W];

	inline void init(int d){
		memset(a,0,sizeof(a));
		a[0]=1;
		for(int i=1;i<=n;++i){
			ll t=m[i]/d;
			ll x=m[i]*t%p;
			ll y=p-(t*(t+1)>>1)%p;
			for(int j=n;j;--j){
				a[j]=(a[j]*x+a[j-1]*y)%p;
			}
			a[0]=a[0]*x%p;
		}
	}

	inline void solve(){
		ll ans=0;
		for(int i=1,nex;i<=M;i=nex+1){
			nex=M;
			for(int j=1;j<=n;++j){
				nex=min(nex,m[j]/(m[j]/i));
			}
			init(i);
			for(int j=0;j<=n;++j){
				ans+=a[j]*add(sum[nex][c][j],p-sum[i-1][c][j]);
			}
		}
		printf("%lld\n",ans%p);
	}

}Q[1005];

int main(){
//	freopen(".in","r",stdin);
//	freopen(".out","w",stdout);
	int T,MaxM=0;r(T);
	for(int i=1;i<=T;++i){
		r(Q[i].n),r(Q[i].c);Q[i].M=N;
		for(int j=1;j<=Q[i].n;++j){
			r(Q[i].m[j]);
			MaxM=max(MaxM,Q[i].m[j]);
			Q[i].M=min(Q[i].M,Q[i].m[j]);
		}
	}
	pre(MaxM);
	for(int i=1;i<=T;++i)Q[i].solve();
	return 0;
}