写这份代码出了一堆bug，被 (color{grey} \_ exttt{nobody}\_) D没了/kk
原式
[egin{aligned}&=sum_{i=1}^{n} sum_{j=1}^{n} (i+j)^k mu(gcd(i,j))^2 gcd(i,j)\&=sum_{d=1}^{n} mu(d)^2 dsum_{i=1}^{n} sum_{j=1}^{n} (i+j)^k[gcd(i,j)==d]\&=sum_{d=1}^{n} mu(d)^2 dsum_{i=1}^{frac{n}{d}} sum_{j=1}^{frac{n}{d}} (id+jd)^k[gcd(i,j)==1]\&=sum_{d=1}^{n} mu(d)^2 d^{k+1}sum_{i=1}^{frac{n}{d}} sum_{j=1}^{frac{n}{d}} (i+j)^k [gcd(i,j)==1]\&=sum_{d=1}^{n} mu(d)^2 d^{k+1} sum_{i=1}^{frac{n}{d}} sum_{j=1}^{frac{n}{d}} (i+j)^k sum_{x|gcd(i,j)}mu(x)\&=sum_{d=1}^{n} mu(d)^2 d^{k+1} sum_{x=1}^{frac{n}{d}}mu(x) sum_{i=1}^{frac{n}{dx}} sum_{j=1}^{frac{n}{dx}} (ix+jx)^k\&=sum_{d=1}^{n} mu(d)^2 d^{k+1} sum_{x=1}^{frac{n}{d}}mu(x)x^k sum_{i=1}^{frac{n}{dx}} sum_{j=1}^{frac{n}{dx}} (i+j)^k\&=sum_{d=1}^{n} mu(d)^2 d^{k+1} sum_{T=1,d|T}^{n}mu(dfrac{T}{d})(dfrac{T}{d})^k sum_{i=1}^{frac{n}{T}} sum_{j=1}^{frac{n}{T}} (i+j)^k\&=sum_{T=1}^{n} sum_{i=1}^{frac{n}{T}} sum_{j=1}^{frac{n}{T}} (i+j)^k sum_{d=1,d|T}^{n} mu(dfrac{T}{d}) (dfrac{T}{d})^k mu(d) d^k dmu(d)\&=sum_{T=1}^{n} T^k sum_{i=1}^{frac{n}{T}} sum_{j=1}^{frac{n}{T}} (i+j)^k sum_{d=1,d|T}^{n} dmu^2(d) mu(dfrac{T}{d})\end{aligned}
]
(T^k) 非常显然的完全积性函数可以线性筛。
考虑如何快速求出 (f(T)=sum_{d=1,d|T}^{n} dmu^2(d) mu(dfrac{T}{d}))
(color{grey} \_ exttt{nobody}\_) 告诉我可以发掘 (mu) 的性质。
那个式子是一堆积性函数狄利克雷卷积起来，还是积性函数。
线性筛当 (p|T) 时考虑进行以下分类讨论进行转移
[f(p^k)=egin{cases}1 (k=0)\1 imes mu^2(1) imes mu(p)+p imes mu(p)^2 imesmu(1)=p-1 (k=1)\1 imes mu^2(1) imes mu(p^2)+p^2 imes mu(p^2)^2 imesmu(1)+p imes mu^2(p) imes mu(p)=-p (k=2)\0 (k>2)end{cases}
]
即
[f(p^k)=egin{cases}1 (k=0)\p-1 (k=1)\-p (k=2)\0 (k>2)end{cases}
]
于是可以线性筛了。
现在式子只剩
[sum_{i=1}^{n} sum_{j=1}^{n} (i+j)^k
]
考虑记一个前缀和 (s_n=sumlimits_{i=1}^{n}i^k)
那么原式 (=sumlimits_{i=n+1}^{2*n}s_i-sumlimits_{i=1}^{n}s_i) ，统计一下 (s_i) 的前缀和就可以 (O(1)) 了。
至此，每一个部分都已经预处理至 (O(1)) ，外层套个整除分块就可以单次 (O(sqrt{n})) 回答询问，于是加强版也能过了。
简单版强样例：
Input
5000000 1000000000000000000
Output
305836999
简单版：
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef double db;
#define pb(x) push_back(x)
#define mkp(x,y) make_pair(x,y)
//#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
//char buf[1<<21],*p1=buf,*p2=buf;
inline int read() {
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)) {if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch))x=x*10+(ch^48),ch=getchar();
return x*f;
}
const int N=10000005;
const int mod = 998244353;
int n;LL k;
int pw[N],S[N],f[N], pri[N], cnt;
bool vis[N];
int qpow(int n,int k,int res=1){
for(;k;k>>=1,n=1ll*n*n%mod)if(k&1)res=1ll*res*n%mod;
return res;
}
void init(const int&N){
pw[1]=1,f[1]=1;
for(int i=2;i<=N;++i){
if (!vis[i]) pri[++cnt] = i, pw[i] = qpow(i, k), f[i] = i - 1;
for (int j = 1; j <= cnt && i * pri[j] <= N; ++ j) {
vis[i * pri[j]] = 1, pw[i * pri[j]] = 1ll * pw[i] * pw[pri[j]] % mod;
if (i % pri[j] == 0) {
if(i / pri[j] % pri[j]) f[i * pri[j]] = -1ll * f[i / pri[j]] * pri[j] % mod;
break;
}
pw[i * pri[j]] =1ll * pw[i] * pw[pri[j]] % mod, f [i * pri[j]] = 1ll * f[i] * f[pri[j]] % mod;
}
}
for(int i=1;i<=N;++i)f[i]=(1ll*pw[i]*f[i]%mod+f[i-1])%mod,pw[i]=(pw[i]+pw[i-1])%mod;
for(int i=1;i<=N;++i)S[i]=(S[i-1]+pw[i])%mod;
}
int query(int n){
int res=0;
for(int l=1,r;l<=n;l=r+1)
r=n/(n/l),res=(res+1ll*(S[(n/l)*2]-S[n/l]*2)*(f[r]-f[l-1])%mod)%mod;
return (res+mod)%mod;
}
signed main(){
n=read(), scanf("%lld", &k), k %= mod - 1,init(n<<1);
printf("%d
",query(n));
}
加强版：
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef double db;
typedef unsigned int unt;
#define pb(x) push_back(x)
#define mkp(x,y) make_pair(x,y)
//#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
//char buf[1<<21],*p1=buf,*p2=buf;
inline unt read() {
unt x=0;char ch=getchar();
while(!isdigit(ch))ch=getchar();
while(isdigit(ch))x=x*10+(ch^48),ch=getchar();
return x;
}
const int N=20000005;
unt T,n,k;
unt pw[N],f[N], pri[1270610], cnt;
bool vis[N];
unt qpow(unt n,unt k,unt res=1){
for(;k;k>>=1,n=n*n)if(k&1)res=res*n;
return res;
}
void init(const unt&N){
pw[1]=1,f[1]=1;
for(unt i=2;i<=N;++i){
if (!vis[i]) pri[++cnt] = i, pw[i] = qpow(i, k), f[i] = i - 1;
for (unt j = 1; j <= cnt && i * pri[j] <= N; ++ j) {
vis[i * pri[j]] = 1, pw[i * pri[j]] = pw[i] * pw[pri[j]];
if (i % pri[j] == 0) {
if(i / pri[j] % pri[j]) f[i * pri[j]] = - f[i / pri[j]] * pri[j];
break;
}
pw[i * pri[j]] = pw[i] * pw[pri[j]], f [i * pri[j]] = f[i] * f[pri[j]];
}
}
for(unt i=1;i<=N;++i)f[i]=pw[i]*f[i]+f[i-1],pw[i]=pw[i]+pw[i-1];
for(unt i=1;i<=N;++i)pw[i]=pw[i-1]+pw[i];
}
unt query(unt n){
unt res=0;
for(unt l=1,r;l<=n;l=r+1)
r=n/(n/l),res=(res+(pw[(n/l)<<1]-(pw[n/l]<<1))*(f[r]-f[l-1]));
return res;
}
signed main(){
T = read(), n=read(), k = read(),init(n<<1);
while(T--)printf("%u
",query(read()));
}
混合码风.jpg