洛谷 P6222 - 「P6156 简单题」加强版（莫比乌斯反演）
题意：
\(T\) 组数据，求 \(\sum\limits_{i=1}^n\sum\limits_{j=1}^n(i+j)^k\mu^2(\gcd(i,j))\gcd(i,j)\)
弱化版中 \(T=1\)，\(n \leq 5 \times 10^6\)
强化版中 \(T=10^4\)，\(n \leq 10^7\)
推式子：
令 \(s(x)=\sum\limits_{i=1}^x\sum\limits_{j=1}^x(i+j)^k\)
令 \(f(t)=\sum\limits_{d|t}d\mu^2(d)\mu(\frac{t}{d})\)
预处理 \(s(x)\) 和 \(g(x)=\sum\limits_{i=1}^xf(i)i^k\)，就可以使用整除分块在 \(\sqrt{n}\) 的时间内求出。
接下来我们的问题就是如何求出 \(s(x)\) 和 \(g(x)\)。
首先预处理 \(i^k\) 肯定是需要的。不过一个个快速幂会超时，不过发现 \(i^k\) 是一个积性函数，欧拉筛解决，这部分时间复杂度 \(\pi(n)\log k\)。
对于 \(s(x)\)，直接求肯定不太容易，我们不妨转化为枚举 \(t \in [2,2x]\)，观察 \(t^k\) 被贡献了几次。
例如 \(x=4\)，\(s(4)=1 \times 2^k+2 \times 3^k+3 \times 4^k+4 \times 5^k+3 \times 6^k+2 \times 7^k+1 \times 8^k\)
观察每一项前面的系数，发现了什么。呈阶梯状分部！
预处理 \(w(x)=\sum\limits_{i=1}^xi^x \times (x-i+1)\)，稍微画个图就会发现 \(s(x)=w(2x)-2w(x)\)
接下来是 \(g(x)\)，要求出 \(g(x)\)，肯定要先求出 \(f(x)\)。
由于 \(f(x)\) 是 \(x\mu^2(x)\) 与 \(mu(x)\) 的狄利克雷卷积，而两项都是积性函数，故 \(f(x)\) 也是积性函数。
因此可以用欧拉筛求出 \(f(x)\)。假设 \(x\) 质因数分解里面有一项 \(p^q\)，我们考虑这一项对答案的贡献。
若 \(q=1\)，\(1\times\mu^2(1)\times\mu(p)+p\times\mu^2(p)\times\mu(1)=p-1\)
若 \(q=2\)，\(1\times\mu^2(1)\times\mu(p^2)+p\times\mu^2(p)\times\mu(p)+p^2\times\mu^2(p^2)\times\mu(1)=-p\)
若 \(q\geq 3\)，那么 \(d\) 与 \(\frac{p^q}{d}\) 中必定有一项的幂 \(\geq 2\)，故 \(f(p^q)=0\)。
\(s(x)\) 与 \(g(x)\) 都求出来了，本题也就迎刃而解了。
/*
Contest: -
Problem: P6222
Author: tzc_wk
Time: 2020.9.24
*/
#include <bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define pb push_back
#define fz(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
#define foreach(it,v) for(__typeof(v.begin()) it=v.begin();it!=v.end();it++)
#define all(a) a.begin(),a.end()
#define fill0(a) memset(a,0,sizeof(a))
#define fill1(a) memset(a,-1,sizeof(a))
#define fillbig(a) memset(a,0x3f,sizeof(a))
#define y1 y1010101010101
#define y0 y0101010101010
#define int unsigned int
typedef pair<int,int> pii;
typedef long long ll;
inline int read(){
int x=0,neg=1;char c=getchar();
while(!isdigit(c)){
if(c=='-') neg=-1;
c=getchar();
}
while(isdigit(c)) x=x*10+c-'0',c=getchar();
return x*neg;
}
int T=read(),N=read(),k=read();
inline int qpow(int x,int e){
int ans=1;
while(e){
if(e&1) ans=ans*x;
x=x*x;e>>=1;
} return ans;
}
int pr[20000005],pcnt=0,f[20000005],p[20000005];
bool vis[20000005];
inline void prework(int n){
f[1]=p[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i]){pr[++pcnt]=i;f[i]=i-1;p[i]=qpow(i,k);}
for(int j=1;j<=pcnt&&pr[j]*i<=n;j++){
p[pr[j]*i]=p[i]*p[pr[j]];vis[pr[j]*i]=1;
if(i%pr[j]) f[pr[j]*i]=f[pr[j]]*f[i];
else{
int lft=i/pr[j];
if(lft%pr[j]) f[pr[j]*i]=-pr[j]*f[lft];
else f[pr[j]*i]=0;
break;
}
}
}
fz(i,1,n) f[i]=f[i-1]+f[i]*p[i];
fz(i,1,n) p[i]+=p[i-1];
fz(i,1,n) p[i]+=p[i-1];
}
inline int sum(int x){
return (p[x<<1]-(p[x]<<1));
}
signed main(){
prework(N<<1);
while(T--){
int n=read(),ans=0;
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
ans+=(f[r]-f[l-1])*sum(n/l);
}
cout<<ans<<endl;
}
return 0;
}
