Misaka Network 与求和
一方通行成功接入了 Misaka Network。
现在他要使用超能力，自然计算式被送到了御坂网络进行处理。这次的计算式是这样子的：
[sum_{i=1}^{N}sum_{j=1}^{N}f(gcd(i,j))^k mod 2^{32}
]
其中 (f(x)) 表示 (x) 次大的质因数，重复的质因数计算多次，例如 (f(6)=2,f(4)=2)。规定 (f(1)=0,f(p)=1)，其中 (p) 为质数。
但是妹妹们都不会算这个式子……所以御坂 20001 号找到了你，希望你帮她算一下。
对于所有数据 (1 leq N,k leq 2 imes 10^9)。
题解
http://jklover.hs-blog.cf/2020/04/19/Loj-572-Misaka-Network-与求和/#more
为了方便，记(f(i))就表示次大质因子的(k)次幂。
[egin{aligned}
ext{ans}&=sum_{i=1}^nsum_{j=1}^n f(gcd(i,j)) \
&=sum_{d=1}^n f(d) sum_{i=1}^{lfloorfrac n d
floor}sum_{j=1}^{^{lfloorfrac n d
floor}} [gcd(i,j)=1] \
&=sum_{d=1}^nf(d)sum_{t=1}^{lfloorfrac n d
floor}mu(t)lfloorfrac n{dt}
floor^2\
&=sum_{d=1}^n lfloorfrac{n}{d}
floor^2sum_{t|d}mu(t)f(frac{d}{t})\
&=sum_{d=1}^n lfloorfrac{n}{d}
floor^2cdot (f*mu)(d)
end{aligned}
]
对(d)整除分块，就需要多次询问(f* mu)的前缀和。记(F(n)=sum_{i=1}^n (f* mu)(i))。
用杜教筛的套路，可以得出
[sum_{i=1}^n f(n) = sum_{i=1}^n (f * mu * I) (i)\
= sum_{i=1}^n (F * I) (i) = sum_{i=1}^n F(lfloor frac n i
floor)\
]
[F(n)=sum_{i=1}^n f(i)-sum_{i=2}^n F(lfloor frac n i
floor)
]
计算(F)时整除分块，递归下去计算。
每次还会询问(f)的前缀和，可以把Min_25筛魔改一下，变成
-
(G(n,m))表示(leq n)的合数，次大质因子(geq p_m)的(f)的和。
-
(H(n))表示(leq n)的质数的(f)的和。
这样就比较好处理了。
因为我偷懒没有预处理杜教筛，所以时间复杂度(O(n^{frac 34}))，相当的慢。
IN uint fpow(uint a,int b){
uint ans=1;
for(;b;b>>=1,a=a*a)
if(b&1) ans=ans*a;
return ans;
}
CO int N=2e5+10;
int n,m,prime[N],num; uint kth[N];
int pos[N],idx,ref1[N],ref2[N]; uint H[N];
uint G(int x,int j){
if(x<=1 or prime[j]>x) return 0;
uint ans=0;
for(int i=j;i<=num and prime[i]*prime[i]<=x;++i)
for(int64 pw=prime[i];pw*prime[i]<=x;pw*=prime[i]){
int k=x/pw;
k=k<=m?ref1[k]:ref2[n/k];
ans+=G(x/pw,i+1)+kth[prime[i]]*(H[k]-i+1);
}
return ans;
}
uint _F[N];
uint F(int x){
if(x==0) return 0;
int k=x<=m?ref1[x]:ref2[n/x];
if(_F[k]) return _F[k];
uint ans=G(x,1)+H[k];
for(int l=2,r;l<=x;l=r+1){
r=x/(x/l);
ans-=F(x/l)*(r-l+1);
}
return _F[k]=ans;
}
int main(){
read(n),m=ceil(sqrt(n));
int K=read<int>();
kth[1]=1;
for(int i=2;i<=m;++i){
if(!prime[i]) prime[++num]=i,kth[i]=fpow(i,K);
for(int j=1;j<=num and i*prime[j]<=m;++j){
prime[i*prime[j]]=1;
if(i%prime[j]==0) break;
}
}
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
pos[++idx]=n/l;
n/l<=m?ref1[n/l]=idx:ref2[n/(n/l)]=idx;
}
for(int i=1;i<=idx;++i) H[i]=pos[i]-1;
for(int j=1;j<=num;++j)
for(int i=1;i<=idx and prime[j]*prime[j]<=pos[i];++i){
int k=pos[i]/prime[j];
k=k<=m?ref1[k]:ref2[n/k];
H[i]-=H[k]-(j-1);
}
uint ans=0;
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
ans+=(F(r)-F(l-1))*(n/l)*(n/l);
}
printf("%u
",ans);
return 0;
}