题意
给定 (n,x,k)，求：
[sum_{a_1=1}^{n}{sum_{a_2=1}^{n}{...sum_{a_x=1}^{n}{left( prod_{j=1}^{x}{a_j^k}
ight) f(gcd(a_1,a_2,...,a_x))*gcd(a_1,a_2,...,a_x)}}}
]
其中，
[f(x)=egin{cases}
0 & exists k>1,x\%k^2=0\
1 & others
end{cases}
]
(1≤t≤10^4,1≤k≤10^9,1≤x≤10^9,1≤n≤2×10^5)
分析
推导
[egin{align}
ans & = sum_{d=1}^{n}{d·f(d)sum_{a_1=1}^{n}{sum_{a_2=1}^{n}{...sum_{a_x=1}^{n}{left( prod_{j=1}^{x}{a_j^k}
ight) ·[gcd(a_1,a_2,...,a_x)==d]}}}}\
& = sum_{d=1}^{n}{d^{kx+1}·f(d)sum_{a_1=1}^{lfloor n/d
floor}{sum_{a_2=1}^{lfloor n/d
floor}{...sum_{a_x=1}^{lfloor n/d
floor}{left( prod_{j=1}^{x}{a_j^k}
ight) ·[gcd(a_1,a_2,...,a_x)==1]}}}}\
& = sum_{d=1}^{n}{d^{kx+1}·f(d)sum_{a_1=1}^{lfloor n/d
floor}{sum_{a_2=1}^{lfloor n/d
floor}{...sum_{a_x=1}^{lfloor n/d
floor}{left( prod_{j=1}^{x}{a_j^k}
ight) ·sum_{t|a_1,t|a_2,...,t|a_x}{mu(t)}}}}}\
& = sum_{d=1}^{n}{d^{kx+1}·f(d)sum_{a_1=1}^{lfloor n/d
floor}{a_1^k sum_{a_2=1}^{lfloor n/d
floor}{a_2^k ...sum_{a_x=1}^{lfloor n/d
floor}{a_x^k ·sum_{t|a_1,t|a_2,...,t|a_x}{mu(t)}}}}}\
& = sum_{d=1}^{n}{d^{kx+1}·f(d) left(sum_{i=1}^{lfloor n/d
floor}{i^k}
ight)^x sum_{t|a_1,t|a_2,...,t|a_x}{mu(t)}}\
& = sum_{d=1}^{n}{d^{kx+1}·f(d) sum_{t=1}^{lfloor n/d
floor}{mu(t)} ·left(sum_{i=1}^{lfloor n/dt
floor}{(it)^k}
ight)^x }\
end{align}
]
令 (T=dt)，并改变枚举顺序
[egin{align}
ans & = sum_{T=1}^{n}{left( sum_{i=1}^{lfloor n/T
floor}{i^k}
ight)^x · T^{kx} sum_{d|T}{f(d)mu(frac{T}{d})d}}\
end{align}
]
令 (g(T)=sum_{d|T}^{T}{f(d)mu(frac{T}{d})d})，最终结果为：
[ans= sum_{T=1}^{n}{left( sum_{i=1}^{lfloor n/T
floor}{i^k}
ight)^x · T^{kx}·g(T)}\
]
代码实现
其中，(g,f) 函数可以 (O(nlogn)) 预处理出，(O(n)) 预处理出 (T^{kx}·g(T)) 的前缀和，然后利用分块对每个 (n) 进行 (O(sqrt{n})) 求解，总的复杂度为：(O(nlogn+tsqrt{n}))。一开始想预处理出全部的答案，但是超时了。
代码：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=1e9+7;
const int N=2e5+5;
int t,k,x,n;
int prime[N],mu[N],cnt,f[N];
bool vis[N];
ll g[N],sum[N],pre[N];
void init()
{
memset(vis,false,sizeof(vis));
int maxn=2e5;
cnt=0;
mu[1]=1;
for(int i=2;i<=maxn;i++)
{
if(!vis[i])
{
prime[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt&&i*prime[j]<=maxn;j++)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
else
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<=maxn;i++) f[i]=1;
for(int i=2;i*i<=maxn;i++)//预处理f(k)
{
int t=i*i;
for(int j=t;j<=maxn;j+=t)
f[j]=0;
}
//预处理g函数
for(int i=1;i<=maxn;i++)//枚举d
{
for(int j=i;j<=maxn;j+=i)//T
g[j]=(g[j]+f[i]*mu[j/i]*i%mod+mod)%mod;
}
}
ll power(ll a,ll b)
{
ll res=1;
a%=mod;
while(b)
{
if(b&1) res=res*a%mod;
a=a*a%mod;
b>>=1;
}
return res;
}
void solve()
{
int maxn=2e5;
ll tmp=0;
for(int i=1;i<=maxn;i++)
{
ll t=power(1LL*i,1LL*k);
tmp=(tmp+t)%mod;
sum[i]=power(tmp,1LL*x)%mod;
pre[i]=(pre[i-1]+power(t,1LL*x%(mod-1))*g[i]%mod+mod)%mod;//T^(kx)*g(T)的前缀和
}
}
int main()
{
init();
scanf("%d%d%d",&t,&k,&x);
solve();
while(t--)
{
scanf("%d",&n);
ll ans=0;
for(int l=1,r;l<=n;l=r+1)
{
r=min(n,n/(n/l));
ans=(ans+sum[n/l]*(pre[r]-pre[l-1])%mod+mod)%mod;
}
printf("%lld
",ans);
}
return 0;
}
参考博客：https://blog.csdn.net/weixin_44282912/article/details/107844614