Powerful number 筛可以用来求出一类积性函数的前缀和，最快可以达到根号复杂度。
- 以下字母大都代表整数；
- 以下 (p) 代表质数；
- 以下 (*) 代表狄利克雷卷积；
- 以下 (a/b=lfloor{aover b} floor)；
以下任何不够严谨的地方都请忽略。
定义
定义一个 Powerful number 为任意一个质因子的次数都不小于 2 的数，如 (72=2^3 imes 3^2)，(4000=2^5 imes5^3)。Powerful number 有如下性质：
- 一个 Powerful number 一定可以表示为 (a^2b^3) 的形式，如 (72=3^2 imes2^3)，(4000=2^{(2+3)} imes5^3=2^2 imes(2 imes5)^3=2^2 imes10^3)；
- (n) 以内的 Powerful number 个数是 (O(sqrt n)) 的：[int_1^{sqrt n}sqrt[3]{nover x^2}dx=O(sqrt n) ]
算法
假设有一个积性函数 (f)，我们希望快速求出 (f) 的前缀和。我们尝试找到另一个积性函数 (g) 使得 (g(p)=f(p))，且可以较快求出 (g) 的前缀和。
再找到一个积性函数 (h)，使得 (f=g*h)。
显然 (f(p)=g(1)h(p)+g(p)h(1)=h(p)+g(p)=h(p)+f(p))，所以 (h(p)=0)。又因为 (h) 为积性函数，所以所有非 Powerful number 的数的 (h) 值都为 0。
因此枚举 (sqrt n) 个 Powerful number 的 (h) 值并求出对应的 (g) 的前缀和便能求出 (f) 的前缀和。
难点在于如何构造 (g,h)，下面举两个例子。
例
例 1
积性函数 (f(p^q)=p^k)（(k) 为常数）。
构造 (g(n)=n^k)，显然 (f(p)=g(p)=p^k)。
考虑找到 (g^{-1}) 使 (g^{-1}*g=epsilon)，则 (h=f*g^{-1})。
假如 (q>0)，还可以推到：
从 (q=0) 开始手玩，很容易发现 (g^{-1}(p^0)=1)，(g^{-1}(p^1)=-p^k)，(g^{-1}(p^q)=0(q>1))。手玩一下 (f*g^{-1}) 可得 (h=p^k-p^{2k})。
例 2
积性函数 (f(p^q)=poplus q)。（LOJ6053 简单的函数）
显然 (f(p)=egin{cases}p-1quad&(p ot=2)\p+1&(p=2)end{cases})。
构造 (g(n)=egin{cases}varphi(n)quad&(2 mid n)\3varphi(n)quad&(2mid n)end{cases})。
下面主要解决两个问题：(g^{-1}) 与 (g) 的前缀和。
用类似于例 1 的方法推导一番可得：(q) 足够大时，(g*g^{-1}(p^q)=egin{cases} g^{-1}(p^q)-g^{-1}(p^{q-1})+pcdot g*g^{-1}(p^{q-1})quad(p
ot =2)\ g^{-1}(p^q)+g^{-1}(p^{q-1})+pcdot g*g^{-1}(p^{q-1})quad(p=2)end{cases})。
(g^{-1}(p^q)=egin{cases} g^{-1}(p^{q-1})quad&(p
ot =2)\ -g^{-1}(p^{q-1})quad&(p=2)end{cases})。
所以：
- (g^{-1}(1)=1)
- (g^{-1}(p^q)=egin{cases}-p+1quad&(p ot =2)\(-1)^q(p+1)quad&(p=2)end{cases}quad(q>0))
然后推导 (g) 的前缀和：(sumlimits_{i=1}^{n}g(i)=sumlimits_{i=1}^nvarphi(i)+2sumlimits_{i=1}^{n/2}varphi(2i))
设 (S'(n)=sumlimits_{i=1}^{n}varphi(2i))。则当 (2mid n) 时：
假如 (2 mid n)，(S'(n)=S'((n-1)/2)+sumlimits_{i=1}^{n-1}varphi(i)+varphi(2n)=S'(n/2)+sumlimits_{i=1}^nvarphi(i))（完 全 一 致）。
因此在杜教筛出 (O(sqrt n)) 个点的 (varphi) 的前缀和后，再处理出 (O(sqrt n)) 个 (S') 的值。(sumlimits_{i=1}^{n}g(i)=sumlimits_{i=1}^nvarphi(i)+2S'(n/2))
于是我们得到了 (g) 的前缀和，并且可以通过暴力求出 (f*g^{-1}(p^q)) 得到 (h(p^q)) 的值，进一步在搜索时同步求出 (h(x))。
代码：
/**********
Author: WLBKR5
Problem: loj 6053
Name: 简单的函数
Source: uploaded by fjzzq2002
Algorithm: powerful number
Date: 2020/06/02
Statue: accepted
Submission: loj.ac/submission/823642
**********/
#include<bits/stdc++.h>
using namespace std;
#define ll long long
ll getint(){
ll ans=0,f=1;
char c=getchar();
while(c<'0'||c>'9'){
if(c=='-')f=-1;
c=getchar();
}
while(c>='0'&&c<='9'){
ans=ans*10-'0'+c;
c=getchar();
}
return ans*f;
}
const int mod=1e9+7,N=1e7+10;
ll n;int sq;
bool boo[N+10];
int pri[1000100],phi[N+10],cnt=0;
void sieve(int n){
//线性筛
boo[1]=1;
phi[1]=1;
for(int i=2;i<=n;i++){
if(!boo[i]){
pri[cnt++]=i;
phi[i]=i-1;
}
for(int j=0;j<cnt&&pri[j]*i<=n;++j){
boo[pri[j]*i]=1;
if(i%pri[j]==0){
phi[pri[j]*i]=phi[i]*pri[j];
break;
}else{
phi[pri[j]*i]=phi[i]*(pri[j]-1);
}
}
}
for(int i=1;i<=n;i++){
phi[i]=(phi[i]+phi[i-1])%mod;
}
}
int calc(ll x){
//x*(x+1)/2
int ans=x%mod*((x+1ll)%mod)%mod;//此处处理不好可能炸 long long
if(ans&1)ans+=mod;
return ans>>1;
}
ll s[10010];//sum_1^{n/i} g(i)
int s_1[200010],s_2[200010];
ll ans=0;
int h[30100][66],d[30100];
int& s_(ll x){
//S'
return x<=200000?s_1[x]:s_2[n/x];
}
ll sumphi(ll x){
//phi 的前缀和
return x<=N?phi[x]:s[n/x];
}
ll sumg(ll x){
//g 的前缀和
return (sumphi(x)+2*s_(x/2))%mod;
}
void dfs(ll x,ll v,ll p){
//枚举 x
//v=h(x)
//h(x)*sum_1^{n/x}g(j)
ans=(ans+v*1ll*sumg(n/x)%mod)%mod;
for(int i=p;i<cnt;i++){
if(x*1ll*pri[i]*pri[i]>n)break;
ll y=x*pri[i];
for(int j=1;y<=n;++j,y*=pri[i]){
if(j>d[i]){
++d[i];
if(pri[i]==2){
h[i][j]=((pri[i]^j)+(j%2?mod-pri[i]-1ll:pri[i]+1ll))%mod;
for(int k=1;k<j;k++){
h[i][j]=(h[i][j]+
(pri[i]^k)*((j-k)%2?mod-pri[i]-1ll:pri[i]+1ll)%mod)%mod;
}
}else{
h[i][j]=((pri[i]^j)+1ll-pri[i]+mod)%mod;
for(int k=1;k<j;k++){
h[i][j]=(h[i][j]+(pri[i]^k)*(mod-pri[i]+1ll)%mod)%mod;
}
}
}
if(!h[i][j])continue;
dfs(y,v*1ll*h[i][j]%mod,i+1);
}
}
}
signed main(){
n=getint();
sieve(N);//预处理质数、phi 的前缀和
for(int i=1;i<=cnt&&pri[i]*1ll*pri[i]<=n;i++){
h[i][0]=1;
}
int u=1;while(n/u>=N)++u;
for(int i=u;i>=1;--i){
//杜教筛筛出 phi 的前缀和
//phi*1=Id
ll n_=n/i;
s[i]=calc(n_);
for(ll l=2,r=0;l<=n_;l=r+1){
r=n_/(n_/l);
s[i]=(s[i]-(r-l+1ll)%mod*1ll*sumphi(n_/l)%mod+mod)%mod;
}
}
vector<ll>v;
for(ll l=1,r=0;l<=n;l=r+1){
r=n/(n/l);
v.push_back(n/l);
}
for(int i=v.size()-1;i>=0;i--){
//算出 S' 的前缀和
s_(v[i])=(s_(v[i]/2)+sumphi(v[i]))%mod;
}
dfs(1,1,0);
cout<<(ans%mod+mod)%mod;
return 0;
}
（原发布于 csdn）
