【LOJ#572】Misaka Network 与求和（莫比乌斯反演，杜教筛，min_25筛）
题面
[ans=sum_{i=1}^nsum_{j=1}^n f(gcd(i,j))^k
]
其中(f(x))表示(x)的次大质因子。
题解
这个数据范围不是杜教筛就是(min\_25)筛了吧。。。
看到次大质因子显然要(min\_25)筛了吧。。。
莫比乌斯反演的部分比较简单，懒得写过程了。
[ans=sum_{T=1}^n [frac{n}{T}]^2sum_{d|T}f(d)^k mu(frac{T}{d})
]
后面带个指数好麻烦啊，就假装(f(x)=f(x)^k)吧。。。
显然要求的就是(f)和(mu)狄利克雷卷积的前缀和。。。
令(displaystyle S(n)=sum_{i=1}^n (f*mu)(i))，一脸杜教筛的感觉，类似杜教筛来写式子。
[g(1)S(n)=sum_{i=1}^n (f*mu*g)(i)-sum_{i=2}^n g(i)S([frac{n}{i}])
]
看到(mu)了，直接令(g(x)=1)，((f*u*g)(i)=(f*(u*1))(i)=(f*e)(i)=f(i))。
写出来就是：
[S(n)=sum_{i=1}^nf(i)-sum_{i=2}^n S([frac{n}{i}])
]
然后考虑怎么求(displaystyle sum_{i=1}^n f(i))，一脸(min\_25)筛。
行，本来以为不是(min\_25)筛就是杜教筛，没想到两个一起来。
好了，实现啥的就可以看看代码了。
复杂度因为杜教筛不能提前筛好一部分前缀和，所以似乎是(O(n^{3/4}))？？？
不太会算复杂度，那就当做(O(mbox{跑得过}))了。
#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
using namespace std;
#define ll long long
#define uint unsigned int
#define MAX 100000
int K,blk;uint n;
int w[MAX],id1[MAX],id2[MAX],m;
int pri[MAX],tot;
bool zs[MAX];
uint g[MAX],prik[MAX];
uint fpow(uint a,int b)
{
uint s=1;
while(b){if(b&1)s*=a;a*=a;b>>=1;}
return s;
}
int getid(int x){return (x<=blk)?id1[x]:id2[n/x];}
void pre(int n)
{
for(int i=2;i<=n;++i)
{
if(!zs[i])pri[++tot]=i,prik[tot]=fpow(i,K);;
for(int j=1;j<=tot&&i*pri[j]<=n;++j)
{
zs[i*pri[j]]=true;
if(i%pri[j]==0)break;
}
}
}
uint calc(int x,int y)
{
if(x<=1||pri[y]>x)return 0;
uint ret=(g[getid(x)]-y+1)*prik[y-1];
for(int i=y;i<=tot&&1ll*pri[i]*pri[i]<=x;++i)
{
ll t1=pri[i],t2=1ll*pri[i]*pri[i];
for(int e=1;t2<=x;++e,t1=t2,t2*=pri[i])
ret+=calc(x/t1,i+1)+prik[i];
}
return ret;
}
uint M[MAX];bool vis[MAX];
uint S(int n)
{
if(vis[getid(n)])return M[getid(n)];
uint ret=calc(n,1)+g[getid(n)];
for(int i=2,j;i<=n;i=j+1)
j=n/(n/i),ret-=(j-i+1)*S(n/i);
vis[getid(n)]=true;
return M[getid(n)]=ret;
}
int main()
{
scanf("%u%d",&n,&K);pre(blk=sqrt(n));
for(uint i=1,j;i<=n;i=j+1)
{
j=n/(n/i);w[++m]=n/i;g[m]=w[m]-1;
if(w[m]<=blk)id1[w[m]]=m;
else id2[n/w[m]]=m;
}
for(int j=1;j<=tot&&1ll*pri[j]*pri[j]<=n;++j)
for(int i=1;i<=m&&1ll*pri[j]*pri[j]<=w[i];++i)
g[i]-=g[getid(w[i]/pri[j])]-(j-1);
uint ans=0,lt=0,nw;
for(uint i=1,j;i<=n;i=j+1)
{
j=n/(n/i);nw=S(j);
ans+=(uint)1*(n/i)*(n/i)*(nw-lt);
lt=nw;
}
printf("%u
",ans);
return 0;
}