【LOJ#572】Misaka Network 與求和（莫比烏斯反演，杜教篩，min_25篩）

題面

LOJ

\[ans=\sum_{i=1}^n\sum_{j=1}^n f(gcd(i,j))^k \]

其中\(f(x)\)表示\(x\)的次大質因子。

題解

這個數據范圍不是杜教篩就是\(min\_25\)篩了吧。。。
看到次大質因子顯然要\(min\_25\)篩了吧。。。
莫比烏斯反演的部分比較簡單，懶得寫過程了。

\[ans=\sum_{T=1}^n [\frac{n}{T}]^2\sum_{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\n",ans);
	return 0;
}