[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;
}
