比较厉害的积性函数求和

听 zak 讲的，感觉很厉害。

给定一个积性函数 \(S\)，可以快速计算 \(S(p^k)\)，求 \(\sum\limits_{i=1}^n\sum\limits_{j=1}^mS(ij)\)。

把 \(n,m\) 当作同阶。

我们考虑枚举 \(i,j\) 的 \(\gcd\)。

\(\sum\limits_{g=1}^{\min(n,m)}\sum\limits_{i=1}^{n/g}\sum\limits_{j=1}^{m/g}S(ijg^2)[\gcd(i,j)=1]\)。

我们先考虑如何处理 \(S(ijg^2)\)。

我们记 \(H(i)=\dfrac{S(g^2i)}{S(g^2)}\)，我们知道它是积性的。

考虑 \(a,b\)，满足 \(\gcd(a,b)=1\)。
\(H(a)\times H(b)=\dfrac{S(g^2a)}{S(g^2)}\times \dfrac{S(g^2b)}{S(g^2)}=\dfrac{S(g^2a)S(g^2b)}{S(g^2)^2}\)
又因为对于积性函数，有 \(S(a)\times S(b)=S(\gcd(a,b))\times S(\operatorname{lcm}(a,b))\)。
\(H(a)\times H(b)=\dfrac{S(g^2a)S(g^2b)}{S(g^2)^2}=\dfrac{S(\gcd(g^2a,g^2b))\times S(\operatorname{lcm}(g^2a,g^2b))}{S(g^2)^2}=\dfrac{S(g^2)S(g^2ab)}{S(g^2)^2}=\dfrac{S(g^2ab)}{S(g^2)}=H(ab)\)。

所以我们可以将 \(S(ijg^2)\) 替换成 \(H(ij)S(g^2)\)。

\(\sum\limits_{g=1}^{\min(n,m)}\sum\limits_{i=1}^{n/g}\sum\limits_{j=1}^{m/g}H(ij)S(g^2)[\gcd(i,j)=1]\)。

将 \(S(g^2)\) 提出来：\(\sum\limits_{g=1}^{\min(n,m)}S(g^2)\sum\limits_{i=1}^{n/g}\sum\limits_{j=1}^{m/g}H(ij)[\gcd(i,j)=1]\)。

因为 \(\gcd(i,j)=1\) 才会做出贡献，所有可以认为 \(H(ij)=H(i)H(j)\)。

然后在进行常规的推式子：

\[\begin{aligned} \text{原式}=&\sum\limits_{g=1}^{\min(n,m)}S(g^2)\sum\limits_{i=1}^{n/g}\sum\limits_{j=1}^{m/g}H(i)H(j)[\gcd(i,j)=1]\\ =&\sum\limits_{g=1}^{\min(n,m)}S(g^2)\sum\limits_{i=1}^{n/g}\sum\limits_{j=1}^{m/g}H(i)H(j)\sum\limits_{t|i,t|j}\mu(t)\\ =&\sum\limits_{g=1}^{\min(n,m)}S(g^2)\sum\limits_{t=1}^{\min(n,m)/g}\mu(t)\sum\limits_{i=1}^{n/gt}H(it)\sum\limits_{j=1}^{m/gt}H(jt) \end{aligned} \]

现在的问题就是形如所有的 \(t=1,2\dots k\)，求出 \(\sum\limits_{i=1}^{n/gt}H(it)\)，相当于 \(\sum\limits_{i=1}^{n/g}[t|i]H(i)\)，也就是迪利克雷后缀和的形式，使用迪利克雷后缀和，时间复杂度为 \(O(\dfrac{n}{g}\ln\ln n)\)。

我们现在就需要对 \(i=1,2\dots n/g\)，求出 \(H(i)\) 的点值。

具体的，我们求出所有 \(i=p^k\) 出的点值即可（其中 \(p\) 为质数）。

也就是要计算 \(\dfrac{S(g^2p^k)}{S(g^2)}\) 的值。

我们讨论 \(p\) 和 \(g\) 的关系：

• 如果 \(p|g\)，设 \(p^j|g\) 且 \(p^{j+1}\nmid g\)，则 \(S(g^2p^k)=S(\operatorname{lcm}(g^2,p^{2j+k}))=\dfrac{S(g^2)S(p^{2j+k})}{S(\gcd(g^2,p^{2j+k}))}=\dfrac{S(g^2)S(p^{2j+k})}{S(p^{2j})}\)，有 \(H(p^k)=\dfrac{S(p^{2j+k})}{S(p^{2j})}\)。
• 如果 \(p\nmid g\)，则 \(S(g^2p^k)=S(g^2)\times S(p^k)\)，有 \(H(p^k)=S(p^k)\)。

发现所有 \(H(p^k)\) 都可以快速计算，然后 \(O(n/g)\) 吧所有点值计算出来即可。

时间复杂度为 \(\sum\limits_{i=1}^{\min(n,m)}O(\dfrac{n}{i}\ln\ln n)=O(n\ln n\ln\ln n)\)。

洛谷 P8570 [JRKSJ R6] 牵连的世界

此题中 \(S(x)=\sigma_0(x)\varphi(x)\)。（当让这种做法过不了 \(3\times 10^6\)）。

#include <bits/stdc++.h>
#define Mod 1000000007
using namespace std;
inline void add(int &a,int b){a+=b;if(a>=Mod) a-=Mod;}
inline void del(int &a,int b){a-=b;if(a<0) a+=Mod;}
inline int qpow(int a,int p)
{
	int ret=1;
	for(;p;p>>=1,a=1ll*a*a%Mod)
		if(p&1) ret=1ll*ret*a%Mod;
	return ret;
}
int mu[500010];bool vis[500010];
bool pf[500010];int inv[500010],siz[500010];
int minp[500010],cnt[500010];
int H[500010],H_[500010],gS;
vector<int> prim;
vector<int> pw[500010];

inline int S(int p,int k){return 1ll*(k+1)*qpow(p,k-1)*(p-1)%Mod;}

int n,m,ans,tmp;
int nn,mm;
int main()
{
	cin>>n>>m;
	if(n>m) swap(n,m);

	mu[1]=1;
	for(int i=2;i<=m;i++)
	{
		if(!vis[i])
		{
			prim.push_back(i),mu[i]=-1;
			for(long long j=1;j<=m;j*=i)
				pw[i].push_back((int)(j)),pf[j]=true;
			minp[i]=i,cnt[i]=1;
		}
		for(int j:prim)
		{
			if(i*j>m) break;
			vis[tmp=i*j]=true;
			if(i%j==0)
			{
				mu[tmp]=0,minp[tmp]=j,cnt[tmp]=cnt[i]+1;
				break;
			}
			mu[tmp]=-mu[i],minp[tmp]=j,cnt[tmp]=1;
		}
	}

	for(int g=1;g<=n;g++)
	{
		gS=1;nn=n/g,mm=m/g;
		for(int g_=g;g_!=1;g_/=pw[minp[g_]][cnt[g_]])
			gS=1ll*gS*S(minp[g_],cnt[g_]*2)%Mod,siz[minp[g_]]=cnt[g_],
			inv[minp[g_]]=qpow(S(minp[g_],cnt[g_]*2),Mod-2);

		H[1]=1;
		for(int i=2;i<=mm;i++)
		{
			if(pf[i])
			{
				if(siz[minp[i]]) H[i]=1ll*S(minp[i],cnt[i]+2*siz[minp[i]])*inv[minp[i]]%Mod;
				else H[i]=S(minp[i],cnt[i]);
			}
			else H[i]=1ll*H[i/pw[minp[i]][cnt[i]]]*H[pw[minp[i]][cnt[i]]]%Mod;
		}
		memcpy(H_,H,(nn+1)<<2);

		for(int p:prim)
		{
			if(p>mm) break;
			for(int i=mm/p;i;i--) add(H[i],H[i*p]);
		}

		for(int p:prim)
		{
			if(p>nn) break;
			for(int i=nn/p;i;i--) add(H_[i],H_[i*p]);
		}

		tmp=0;

		for(int t=1;t<=nn;t++)
			if(mu[t])
				(mu[t]==1?add:del)(tmp,1ll*H[t]*H_[t]%Mod);
		add(ans,1ll*gS*tmp%Mod);

		for(int g_=g;g_!=1;g_/=pw[minp[g_]][cnt[g_]])
			siz[minp[g_]]=0;
	}
	cout<<ans;
	return 0;
}