[LOJ 572] Misaka Network 与求和
一、题目
二、解法
直接推柿子吧：
\[\sum_{i=1}^n\sum_{j=1}^nf(\gcd(i,j))^k
\]
\[\sum_{d=1}^nf(d)^k\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}[(i,j)=1]
\]
\[\sum_{d=1}^nf(d)^k\sum_{x=1}^{n/d}\mu(x)\frac{n}{dx}\frac{n}{dx}
\]
\[\sum_{T=1}^n(\frac{n}{T})^2\sum_{d|T}f(d)^k\mu(\frac{T}{d})
\]
上面的式子就是普通的莫比乌斯反演，现在发现有些做不动了，算的时候外面可以套整除分块，但是里面就很难算的。因为 \(n\) 太大了所以不能预处理，我们考虑用筛法来做后面那个部分的前缀和。
设 \(h(i)=\sum_{d|i}f(d)\mu(\frac{i}{d})=f*\mu\) ，如果能求出 \(\sum_{i=1}^n h(i)\) 我们就成功了。考虑杜教筛，令辅助函数 \(g\) 为 \(I\) ，那么 \(g\) 的前缀和是很好求的，\(h*g=f*\mu*I=f\) ，它的前缀和虽然看上去还是很难，但是已经得到简化了。
求 \(f\) 的前缀和可以用 \(\tt Min25\) ，由于是次大的质因数产生贡献，所以我们可以在递归的时候，算每一个比他大的质数和他的贡献，也就是这个质数作为最大的质因子，他作为次大的质因子的贡献。还有一种情况是最大的质因子就是他，但是这要求次数要大于 \(1\) ，递归部分代码可以这样写：
int S(int x,int y)
{
if(x<=p[y]) return 0;//出口
int ans=pk[y]*(g[id(x)]-y)%MOD;//y一定是上一个质因数，这里我们算贡献，g是质数个数
for(int i=y+1;i<=cnt && p[i]*p[i]<=x;i++)//Min25
{
int pr=p[i];
for(int e=1;pr<=x;e++,pr*=p[i])
ans=(ans+S(x/pr,i)+(e!=1)*pk[i])%MOD;//如果次数大于1，那么可以直接贡献
}
return ans;
}
然后我们的杜教筛和 \(S(n,0)\) 都是需要记忆化的，由于我写自然溢出 \(\tt Wa\) 了，所以我只能暴力模了。
#include <cstdio>
#include <iostream>
#include <cmath>
#include <map>
using namespace std;
const int N = 200000;
const int M = N+5;
#define int long long
const int MOD = 4294967296ll;
int read()
{
int x=0,f=1;char c;
while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
return x*f;
}
int n,k,cnt,tot,sqr,p[M],vis[M],pk[M],w[M];
int ans,g[M],p1[M],s1[M],p2[M],s2[M],id1[M],id2[M];
int qkpow(int a,int b)
{
int r=1;
while(b>0)
{
if(b&1) r=r*a%MOD;
a=a*a%MOD;
b>>=1;
}
return r;
}
void init(int n)
{
p[0]=pk[0]=1;
for(int i=2;i<=n;i++)
{
if(!vis[i])
{
p[++cnt]=i;
pk[cnt]=qkpow(i,k);
}
for(int j=1;j<=cnt && i*p[j]<=n;j++)
{
vis[i*p[j]]=1;
if(i%p[j]==0) break;
}
}
}
int id(int x)
{
if(x<=sqr) return id1[x];
return id2[n/x];
}
int S(int x,int y)
{
if(x<=p[y]) return 0;
int ans=pk[y]*(g[id(x)]-y)%MOD;
for(int i=y+1;i<=cnt && p[i]*p[i]<=x;i++)
{
int pr=p[i];
for(int e=1;pr<=x;e++,pr*=p[i])
ans=(ans+S(x/pr,i)+(e!=1)*pk[i])%MOD;
}
return ans;
}
int zy(int n)
{
if(p2[id(n)]) return s2[id(n)];
s2[id(n)]=S(n,0);p2[id(n)]=1;
return s2[id(n)];
}
int get(int n)
{
if(n<=1) return 0;
if(p1[id(n)]) return s1[id(n)];
int ans=zy(n);
for(int l=2,r;l<=n;l=r+1)
{
r=n/(n/l);
ans=(ans-(r-l+1)*get(n/l))%MOD;
}
s1[id(n)]=ans;p1[id(n)]=1;
return ans;
}
signed main()
{
n=read();k=read();
init(N);sqr=sqrt(n);
for(int l=1,r;l<=n;l=r+1)
{
r=n/(n/l);//写错了
w[++tot]=n/l;
g[tot]=w[tot]-1;
if(n/l<=sqr) id1[n/l]=tot;
else id2[n/(n/l)]=tot;
}
for(int i=1;i<=cnt;i++)
for(int j=1;j<=tot && p[i]*p[i]<=w[j];j++)
g[j]-=g[id(w[j]/p[i])]-i+1;
for(int l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
ans=(ans+(n/l)*(n/l)*(get(r)-get(l-1)))%MOD;
}
printf("%lld\n",(ans+MOD)%MOD);
}
