声明:我还没有退役,只是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}
floor}sum_{x_2=1}^{lfloor frac{m_2}{d}
floor}...sum_{x_n=1}^{lfloor frac{m_n}{d}
floor}{[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}
floor}sum_{x_2=1}^{lfloor frac{m_2}{d}
floor}...sum_{x_n=1}^{lfloor frac{m_n}{d}
floor}{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}
floor}{mu(t)}sum_{x_1=1}^{lfloor frac{m_1}{dt}
floor}sum_{x_2=1}^{lfloor frac{m_2}{dt}
floor}...sum_{x_n=1}^{lfloor frac{m_n}{dt}
floor}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}
floor}{mu(t)}sum_{x_1=1}^{lfloor frac{m_1}{dt}
floor}{{(m_1-d*t*x_1)}}sum_{x_2=1}^{lfloor frac{m_2}{dt}
floor}{(m_2-d*t*x_2)}...sum_{x_n=1}^{lfloor frac{m_n}{dt}
floor}{(m_n-d*t*x_n)}
\=sum_{d=1}^{m}{(^{d-1}_{c-2})}sum_{t=1}^{lfloor frac{m}{d}
floor}{mu(t)}prod_{i=1}^{n}{(m_i*lfloor frac{m_i}{dt}
floor-d*t*frac{lfloor frac{m_1}{dt}
floor*(lfloor frac{m_1}{dt}
floor+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}
floor-d*frac{lfloor frac{m_1}{d}
floor*(lfloor frac{m_1}{d}
floor+1)}{2})}}
]
至此
复杂度为(O(Tmn))算下来大概(1e8)
然而还可以再优化
设
[F(d)=prod_{i=1}^{n}{(m_i*lfloorfrac{m_i}{d}
floor-d*frac{lfloorfrac{m_1}{d}
floor*(lfloor frac{m_1}{d}
floor+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^3sqrt{m}))(处理(a)是(n^2),需要算(nsqrt{m})次,因为(d)有(nsqrt{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
",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;
}