利用powerful number求积性函数前缀和小结

介绍

powerful number指的是，一个数质因数分解后，每个质因子的指数 $\ge 2$，那么这个数就叫powerful number。

这种数有个很好的性质：$n$ 以内的powerful number数量是 $\sqrt{n}$ 级别的。于是我们往往可以用这个性质将求前缀和的复杂度变成这个级别。

利用这个东西，我们可以很方便地求一些积性函数的前缀和。

正题

考虑一般情况，对于积性函数 $f(x)$，我们构造一个积性函数 $g(x)$，满足：

1. $g$ 的前缀和很好求。
2. 对于任意质数 $p$，有 $g(p)=f(p)$。

条件比较苛刻，所以应用范围不算太广。

令 $h=f/g$，这里除法是狄利克雷卷积的逆运算。那么有 $f(p)=h(1)g(p)+h(p)g(1)$。

别忘记积性函数都满足 $f(1)=g(1)=h(1)=1$，根据 $f(p)=g(p)$，得到 $h(p)=0$。

由于 $f,g$ 都是积性函数，根据狄利克雷卷积的性质，$h$ 也是积性函数，即对于任意 $h(i)$，当 $i$ 不是powerful number时，$h(i)$ 都为 $0$。

有了这个性质，再看回我们的问题：求 $f$ 的前缀和，即：
$$\begin{array}{rl}\sum _{i=1}^{n}f(i)& =\sum _{ij\le n}g(i)h(j)\\ & =\sum _{i=1}^{n}h(i)\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }g(j)\\ & =\sum _{i=1}^{n}h(i)S_g(\lfloor \frac{n}{i}\rfloor )\end{array}$$

由于 $g$ 的前缀和很好求，根据上面得到的性质，我们只需要枚举 $i$ 为所有powerful number，然后求和即可。

至于求 $h$ 的方法，利用定义即可，对于$h(p^k)$，我们这样求：
$$\begin{array}{rl}f(p^k)& =\sum _{i=0}^{k}h(p^i)g(p^{k-i})\\ h(p^k)& =f(p^k)-\sum _{i=0}^{k-1}h(p^i)g(p^{k-i})\end{array}$$

例1

有一个积性函数 $f(x)$，满足 $f(p^k)=p$，求前缀和。

构造 $g=id$ 即可，$id$ 是单位函数，即 $id(x)=x$，前缀和显然很好求，也满足 $f(p)=g(p)=p$。

例2

多组询问（100组），令 $f(x)$ 表示 $x$ 的最大平方因子，求前 $n$ 项之和，$n\le 10^{14}$。

对于质数 $p$，显然有 $f(p)=1$，于是令 $g=I$ 即可。$I$ 是恒等函数，即 $I(x)=1$。

注意到这里是多组询问，如果每次询问都达到 $\sqrt{n}$ 我们是受不了的，观察答案式子：
$$\sum _{i=1}^{n}h(i)\lfloor \frac{n}{i}\rfloor$$

手玩一下求 $h$ 的式子，会发现在这题有个有趣的性质：当 $k$ 为奇数时 $h(p^k)=0$。

也就是说对于任意 $i$，只有当 $i$ 为平方数时 $h(i)$ 才不为 $0$，所以答案式子等于：
$$\sum _{i=1}^{\sqrt{n}}h(i^2)\lfloor \frac{n}{i^2}\rfloor$$

这个东西再套个除法分块就是 $\sqrt[3]{n}$ 的，总时间复杂度 $O(\sqrt{n}+T\sqrt[3]{n})$。

代码如下：

#include <bits/stdc++.h>
using namespace std;
#define maxn 10000010
#define ll long long

int T;ll n;
int prime[maxn],t=0;
bool v[maxn];
const int N=maxn-10;
const ll N2=1ll*N*N;
void SieveInit(){
	for(int i=2;i<=N;i++){
		if(!v[i])prime[++t]=i;
		for(int j=1;j<=t&&i*prime[j]<=N;j++){
			v[i*prime[j]]=true;
			if(i%prime[j]==0)break;
		}
	}
}
int id(ll x){return x<=N?x:x/N+N;}
ll *h[maxn];unsigned ll H[maxn];//这题对2^64取模所以利用自然溢出
void calc_h(){
	for(int i=1;i<=t;i++){//由于只有平方数有值，所以非平方数干脆不存，h[i][j]是i^{2j}的值
		double lim=N2/prime[i]/prime[i];ll now=1;int j=0;
		while(now<=lim)now*=1ll*prime[i]*prime[i],j++;
		h[i]=new long long[j+1];
		ll d=1;for(int k=0;k<=j;k++){
			h[i][k]=d;for(int l=0;l<k;l++)h[i][k]-=h[i][l];
			d*=1ll*prime[i]*prime[i];
		}
	}
}
void dfs(int x,ll now,ll val){//枚举powerful number
	H[(int)sqrt(now)]=val;
	for(int i=x;i<=t;i++){
		double lim=N2/prime[i]/prime[i];
		if(now>lim)break;
		ll Now=now;
		for(int j=2;;j+=2){
			Now*=1ll*prime[i]*prime[i];
			dfs(i+1,Now,val*h[i][j/2]);
			if(Now>lim)break;
		}
	}
}

int main()
{
	SieveInit();calc_h();dfs(1,1,1);
	for(int i=1;i<=N;i++)H[i]+=H[i-1];
	cin>>T;while(T--)
	{
		cin>>n;
		int sqn=sqrt(n);unsigned ll ans=0;
		for(int l=1,r;l<=sqn;l=r+1){
			r=sqrt(n/(n/l/l));
			ans+=(H[r]-H[l-1])*(n/l/l);
		}
		cout<<ans<<"\n";
	}
}