【2020杭电多校round6】HDU6833 A Very Easy Math Problem
题目大意
给定\(k,x\)。\(t\)次询问，每次给出一个\(n\)。请你求：
\(\gcd(a_1,a_2,\ldots ,a_n)\)是\(a_1,a_2,...,a_n\)的最大公约数。
函数\(f(x)\)定义为：如果存在一个整数\(k\) (\(k>1\))使得\(k^2\)是\(x\)的约数，则\(f(x)=0\)，否则\(f(x)=1\)。
答案对\(10^9+7\)取模。
数据范围：\(1\leq k,x\leq 10^9\)，\(1\leq t\leq 10^4\)，\(1\leq n\leq 2\times 10^5\)。
本题题解
先考虑单组数据的情况。
可以计算每个\(\gcd\)值的贡献。不过精确考虑“\(\gcd\)等于”几有些困难，我们可以先对每个\(d\)求出“\(\gcd\)是\(d\)的倍数”的贡献，再容斥回去：
// F[d] 表示gcd是d的倍数,对答案的贡献
// 容斥成 G[d] 表示gcd就是d,对答案的贡献
int ans=0;
for(int d=n;d>=1;--d){
G[d]=F[d];
for(int j=d+d;j<=n;j+=d){
G[d]-=G[j];
}
ans+=G[d]*f[d]*d;
}
考虑钦定了\(\gcd\)是\(d\)的倍数，如何求出\(\sum_{a_1=1}^{n}\sum_{a_2=1}^{n}\ldots \sum_{a_x=1}^{n}\left (\prod_{j=1}^{x}a_j^k\right )\)，也就是代码片段里的F[d]。此时的限制，就是所有\(a_i\)都要是\(d\)的倍数。可以先令所有\(a_i\)都除以\(d\)，那剩下部分的取值，就只有\([1,\lfloor\frac{n}{d}\rfloor]\)，每个\(a_i\)都能取这些值，那它们的乘积之和，就等于：\((1^k+2^k+\dots+\lfloor\frac{n}{d}\rfloor^k)^x\)。再把所有的\(d\)乘回去，则\(F[d]=d^{kx}(1^k+2^k+\dots+\lfloor\frac{n}{d}\rfloor^k)^x\)。求出F[]后，再用上面的代码容斥求出答案。
但这只是单组数据时的做法。
对多组数据，每次独立求答案是不现实的。我们考虑从\(n-1\)向\(n\)递推答案，每次只需要计算这两者之间微小的变化。
具体来说，要使\(n\)对新答案有影响，则\(d\)一定要是\(n\)的约数（否则\(\lfloor\frac{n}{d}\rfloor=\lfloor\frac{n-1}{d}\rfloor\)，求出来的\(F\)值显然还是一样的）。那么我们每次枚举\(n\)的约数，时间复杂度相当于\(1\dots n\)的约数数量之和，约等于\(O(n\log n)\)（可以考虑每个数能作为哪些数的约数，显然是\(i,2i,3i\dots ,\lfloor\frac{n}{i}\rfloor\cdot i\)，所以数量是\(\frac{n}{i}\)。\(1\dots n\)的总数就是\(\sum_{i=1}^{n}\frac{n}{i}\)，也就是调和级数）。
现在，我们虽然能快速求出哪些\(F\)有变化，也能得到新的\(F\)值，但如果每次重新做一遍容斥，复杂度还是无法接受的。我们考虑把这个容斥转化一下，改成求出每个\(F[d]\)在容斥里的贡献系数\(c[d]\)。这样不管\(F[d]\)怎么变化，我们只需要往答案里加上变化值\(\times c[d]\)即可。
\(c[d]\)怎么求呢？可以先初始化，令\(c[d]=f(d)\cdot d\)（这里的\(f\)是题面中定义的函数，千万不要与我们定义的大\(F\)搞混）。然后：
for(int i=1;i<=n;++i)
c[i]=f[i]*i;
for(int i=1;i<=n;++i){
for(int j=i+i;j<=n;j+=i){
c[j]-=c[i];
}
}
这是因为每个\(F\)值会在它的约数处被减掉。可以把这个代码对照上面的容斥理解，上面的容斥（求出的ans）就等价于\(\sum_{i=1}^{n}F[i]\cdot c[i]\)。这就是我们求出这个系数\(c\)数组的妙处。
于是我们就能\(O(1)\)修改并更新答案了。
总时间复杂度\(O(n\log n+t)\)。
参考代码：
//problem:
#include <bits/stdc++.h>
using namespace std;
#define pb push_back
#define mk make_pair
#define lob lower_bound
#define upb upper_bound
#define fi first
#define se second
#define SZ(x) ((int)(x).size())
typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
template<typename T>inline void ckmax(T& x,T y){x=(y>x?y:x);}
template<typename T>inline void ckmin(T& x,T y){x=(y<x?y:x);}
const int MOD=1e9+7;
inline int mod1(int x){return x<MOD?x:x-MOD;}
inline int mod2(int x){return x<0?x+MOD:x;}
inline void add(int& x,int y){x=mod1(x+y);}
inline void sub(int& x,int y){x=mod2(x-y);}
inline int pow_mod(int x,int i){int y=1;while(i){if(i&1)y=(ll)y*x%MOD;x=(ll)x*x%MOD;i>>=1;}return y;}
const int MAXN=2e5;
int K,X,n;
int pw[MAXN+5],sum[MAXN+5],pwX_of_pwK[MAXN+5],pwX_of_sum[MAXN+5];
int F[MAXN+5],coef[MAXN+5],ans[MAXN+5];
vector<int>factors[MAXN+5];
int p[MAXN+5],cnt_p,f[MAXN+5];
bool v[MAXN+5];
void sieve(int lim){
f[1]=1;
for(int i=2;i<=lim;++i){
if(!v[i]){
p[++cnt_p]=i;
f[i]=1;
}
for(int j=1;j<=cnt_p && i*p[j]<=lim;++j){
v[i*p[j]]=1;
if(i%p[j]==0){
break;
}
if(f[i])
f[i*p[j]]=1;
}
}
}
void init(){
sieve(MAXN);
for(int i=1;i<=MAXN;++i){
pw[i]=pow_mod(i,K);
sum[i]=mod1(pw[i]+sum[i-1]);
pwX_of_pwK[i]=pow_mod(pw[i],X);
pwX_of_sum[i]=pow_mod(sum[i],X);
}
}
int main() {
int T;cin>>T>>K>>X;
init();
/*
//单组数据的简单容斥:
cin>>n;
for(int g=1;g<=n;++g){
F[g]=(ll)pow_mod(sum[n/g],X)*pow_mod(pw[g],X)%MOD;
}
int ans=0;
for(int i=n;i>=1;--i){
G[i]=F[i];
for(int j=i+i;j<=n;j+=i){
sub(G[i],G[j]);
}
ans=((ll)ans + (ll)G[i]*f[i]%MOD*i)%MOD;
}
cout<<ans<<endl;
*/
for(int i=1;i<=MAXN;++i)
coef[i]=f[i]*i;
for(int i=1;i<=MAXN;++i){
for(int j=i+i;j<=MAXN;j+=i){
sub(coef[j],coef[i]);
}
}
ans[1]=F[1]=1;
for(int i=1;i<=MAXN;++i){
for(int j=i;j<=MAXN;j+=i){
factors[j].pb(i);
}
}
for(int i=2;i<=MAXN;++i){
ans[i]=ans[i-1];
for(int j=0;j<SZ(factors[i]);++j){
int g=factors[i][j];
sub(ans[i],(ll)F[g]*coef[g]%MOD);
F[g]=(ll)pwX_of_sum[i/g]*pwX_of_pwK[g]%MOD;
add(ans[i],(ll)F[g]*coef[g]%MOD);
}
}
while(T--){
cin>>n;
cout<<ans[n]<<endl;
}
return 0;
}
