PN筛入门
PN 筛
Powerful Number 筛
oi-wiki 上说是杜教筛的 expand，但是这玩意 只能做积性函数前缀和
若未加说明，以下涉及的 \(f,g,h\) 均为积性函数。\(p\) 均为质数
Powerful number
定义一个数 \(N=\prod_{i=1}^mp_i^{c_i}\) 为 powerful number 当且仅当 \(\forall i,c_i\ge 2\)。特别的，在本文中我们认为 \(1\) 也是powerful number。
性质：
-
\(\forall N\in \mathbb{PN},\exists a,b\in \mathbb{Z},N=a^2b^3\)
证明容易，对于每个 \(c_i\bmod 2=0\)，加入 \(a\)，否则将 \(p_i^{3}\) 加入 \(b\)，余下 \(p_i^{c_i-3}\) 加入 \(a\) 即可
-
\(O(\sum_{i=1}^n[i\in\mathbb{PN}])=O(\sqrt n)\)
Proof:
\(O(\sum_{i=1}^n[i\in\mathbb{PN}])=O(\sum_{a=1}^{\sqrt n}({\frac{n}{a^2})}^{\frac{1}{3}})=O(n^{\frac{1}{3}} * n^{\frac{1}{2}·\frac{1}{3}})=O(\sqrt n)\)
-
\(1\sim n\) 的 $PN $ 每个质因子都不大于 \(\sqrt n\)
算法内容
考虑求出 \(S(n)=\sum_{i=1}^nf(i)\)。
我们考虑构造一个积性函数 \(g\)。（这里与杜教筛的区别在于构造难易程度上）
满足以下要求：
- \(g\) 易求前缀和
- \(\forall p,g(p)=f(p)\)
再构造 \(h * g=f\)。根据定义式有：
函数 \(h\) 有以下性质：
-
\(h(p)=0\)，这是因为 \(f(p)=g(1)h(p)+g(p)h(1),h(1)=g(1)=1,g(p)=h(p)\implies h(p)=0\)
-
\(\forall i\notin\mathbb{PN},h(i)=0\)。根据 \(h\) 为积性函数可知。
-
\(h(p^c)=f(p^c)-\sum_{i=0}^{c-1}h(p^i)g(p^{c-i})\)
这告诉我们：\(1\sim n\) 中的 \(h(p)\) 仅有 \(O(\sqrt n)\) 个有效值
而且只需求出所有 \(h(p_i^{c_i})\) 即可求出所有 \(h(n)\) 。
现在我们设 \(G(n)=\sum_{i=1}^ng(i)\)，可以快速计算 \(G(n)\) 。
则有：
这告诉我们，若可以快速求得 \(G\)，则我们可以 \(O(\sqrt n)\) 筛出所有 \(\mathbb{PN}\)（暴力枚举 \(a,b\) 然后去重即可）。然后快速求出所有有效 \(h\) 值，枚举 \(h\) 即可求出前缀和。
\(h\) 的求值可以考虑定义法，当然更好的是考虑推出通项公式。
整个PN的流程可以梳理为：
-
设出 \(g,g(p^c)\) 必须容易求出（快速计算）
-
逆设出 \(h\)，求出所有 \(\mathbb{PN}\)
-
求 \(G\)
-
求 \(h\)
-
数论分块求出 \(S\)
其中第二步、第五步是 \(O(\sqrt n)\)，问题在于第三步和第四步的求算。
如果可以直接运用公式法得出最好，如果不行，第三步可以考虑杜教筛，第四步考虑根据定义式计算。
考虑第四步的复杂度：\(\le \sqrt n\) 的质因子根据素数分布定理仅有 \(O(\frac{\sqrt n}{\ln n})\) 个，而每个的指数最多有 \(\log n\) 个，因此暴力计算的复杂度是 \(O(\frac{\sqrt n}{\ln n}\log^2 n)\) 的，实际表现更为优异
至于求出所有 \(h\) 有效值，则可以考虑直接从大到小将所有的质因子“拼上去”，类似于DFS。
甚至说如果单次询问不需要记忆化所有 \(PN\) 值，硬算即可。只需要求出所有的 \(h(p^c)\)。
杜教筛与 PN 筛的联系在于求 \(G\)，事实上如果可以构造出杜教筛所要求的 \(g\)，就可以直接利用杜教筛求解而不需要 \(PN\) 筛。
但若构造不出来， \(g\) 是一个简化过程，也许 \(g\) 会更加容易构造出杜教筛所需。
况且一个很优秀的地方是这里面用到的 \(G\) 都是 \(G(\lfloor\frac{n}{i}\rfloor)\)。这意味着杜教筛的总复杂度不变。
其复杂度就是 \(O(\frac{\sqrt n}{\ln n}\log^2 n+n^{\frac{2}{3}})\)。
\(h\) 的构造与求解是一个机械的过程，难点在于设计 \(g\) 与求解 \(G\)
例题：模版-Min_25筛
定义积性函数 \(f\)，且 \(f(p^k)=p^k(p^k-1)\)，求 \(S(n)=\sum_{i=1}^nf(i)\)。
显然构造 \(g(n)=n·\varphi(n)\)，可以在质数处"拟合"。
现在考虑求出 \(G(n)\)。
求解 \(G(n)\)
看到 \(\varphi\) 可以想到设计 \(g'(n)=n\)。
所以 \((g * g')(n)=\sum_{d|n}g(d)\frac{n}{d}=\sum_{d|n}n\varphi(d)=n^2\)。
带入杜教筛的式子，有：
\[\begin{aligned} g'(1)G(n)=G(n)&=\sum_{i=1}^n(g * g')(i)-\sum_{i=2}g'(i)G(\lfloor\frac{n}{i}\rfloor)\\ &=\sum_{i=1}^ni^2-\sum_{i=2}^ni·G(\lfloor\frac{n}{i}\rfloor)\\ \end{aligned} \]
问题至此解决，下附一份板子。
#include<bits/stdc++.h>
using namespace std;
#define N 1050500
#define int long long
const int it=1e6+7,p=1e9+7;
int v[N],pri[N],tot,phi[N],n,m,s[N],inv2,inv6,inv4,g[N];
int power(int a,int b){
int ans=1;
while(b){
if(b&1)ans=ans*a%p;
a=a*a%p;b>>=1;
}
return ans;
}
unordered_map<int,int>sit,h;
void init(){
phi[1]=1;
for(int i=2;i<it;i++){
if(!v[i]){
pri[++tot]=i;phi[i]=i-1;
}
for(int j=1;j<=tot&&i*pri[j]<it;++j){
v[pri[j]*i]=1;
if(i%pri[j]==0){
phi[pri[j]*i]=pri[j]*phi[i];break;
}
phi[pri[j]*i]=phi[pri[j]]*phi[i];
}
}
for(int i=1;i<it;++i)s[i]=(s[i-1]+i%p*phi[i]%p)%p,g[i]=i*phi[i]%p;
}
int ych(int n){
n%=p;
return n*(n+1)%p*inv2%p;
}
int lfh(int n){
n%=p;
return n*n%p*(n+1)%p*(n+1)%p*inv4%p;
}
int pfh(int n){
n%=p;
return n*(n+1)%p*(n+n+1)%p*inv6%p;
}
int calc(int n){
if(n<it)return s[n];
if(sit[n]!=0)return sit[n];
// cout<<"calc: "<<n<<"\n";
int res=pfh(n);int lst=1,cur=0;
for(int l=2,r;l<=n;l=r+1){
r=min(n,n/(n/l));cur=ych(r);
res=(res+p-(cur-lst)%p*calc(n/l)%p)%p;lst=cur;
}res=(res%p+p)%p;
return sit[n]=res;
}
//以上是杜教筛求 G
int res=0;
void sol(int now,int val,int pos){//现在以及算到PN:now，处理到第pos个质数，目前h值为val
res+=val*calc(n/now)%p;res%=p;
// cout<<now<<' '<<val<<" "<<pos<<"\n";
//必须写除法，否则越界死掉
for(int i=pos;i<=tot;i++){
if(now>n/pri[i]/pri[i])return ;
int x=now*pri[i],mul=pri[i];
for(int j=2;j;++j){
if(x>n/pri[i])break;
mul*=pri[i];x*=pri[i];
if(!h[mul]){
int w=mul%p*(mul%p-1)%p;//f(n)
int nowmul=1;
for(int c=0;c<j;++c){
if(c!=1){
int t=mul/nowmul;
w=(w+p-h[nowmul]*((t-t/pri[i])%p)%p*(t%p)%p)%p;
}
nowmul=nowmul*pri[i];
}
h[mul]=w;
}
sol(now*mul,val*h[mul]%p,i+1);
}
}
}
signed main(){
cin>>n;inv2=power(2,p-2);inv4=power(4,p-2);inv6=power(6,p-2);
init();h[1]=1;
sol(1,1,1);
cout<<(res%p+p)%p<<"\n";
}
简单的函数
这告诉我们求G也可以综合利用已知信息
给定函数：
求 \(S(n)=\sum_{i=1}^nf(i)\)
考虑PN筛构造。
注意到 \(f(p)=\varphi(p),p>2,f(2)=3\)。
所以可以构造 \(g(n)=\varphi(n),n\bmod 2=1,g(n)=3\varphi(n),n\bmod 2=0\)。
对于一个积性函数，强行钦定含某个质因子的数多一个带权系数，仍然是积性函数
因此 \(h\) 容易求出。
现在考虑求解 \(G(n)\)。
求解 \(S_{\varphi}(n)\) 是基本功，根据杜教筛可以知道 \(S_{\varphi}(n)=\frac{n(n+1)}{2}-\sum_{i=2}^nS(\lfloor\frac{n}{i}\rfloor)\)
考虑 \(S_{\varphi(2)}(n)\)
所以可以递推求解 \(S_{\varphi(2)}\)。
而且因为求解 \(S_{\varphi}(n)\) 的时候会求解 \(S_{\varphi(i)}(\lfloor\frac{n}{2}\rfloor)\)。所以 \(S_{\varphi(2)}\) 直接记搜的复杂度与杜教筛求解 \(S_{\varphi}\) 的复杂度相同，都是 \(O(n^{\frac{2}{3}})\)
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int it=3e6+7,it2=it>>1;
const int p=1e9+7,inv2=(p+1)>>1;
#define N 3050500
int pi[N],phi[N],n,m,v[N],tot,sph[N],sph2[N];
unordered_map<int,int>s1,s2;
void init(){
phi[1]=1;
for(int i=2;i<it;i++){
if(!v[i]){
pi[++tot]=i;phi[i]=i-1;
}
for(int j=1;j<=tot&&i*pi[j]<it;++j){
v[pi[j]*i]=1;
if(i%pi[j]==0){
phi[pi[j]*i]=pi[j]*phi[i];break;
}
phi[pi[j]*i]=phi[pi[j]]*phi[i];
}
}
for(int i=1;i<it;i++)sph[i]=(sph[i-1]+phi[i])%p;
for(int i=1;i<it2;i++)sph2[i]=(sph2[i-1]+phi[i<<1])%p;
// for(int i=1;i<=20;i++)cout<<sph2[i]<<" ";cout<<"\n";
}
int pfh(int x){
x%=p;
return x*(x+1)%p*inv2%p;
}
int calcS1(int x){
if(x<it)return sph[x];
if(s1[x]!=0)return s1[x];
int res=pfh(x);
// cout<<"! "<<res<<"\n";
for(int l=2,r;l<=x;l=r+1){
r=min(x,x/(x/l));
res-=(r-l+1)%p*calcS1(x/l)%p;
// cout<<(r-l+1)<<" "<<calcS1(x/l)<<"\n";
}
res=(res%p+p)%p;
return s1[x]=res%p;
}
int calcS2(int x){
if(x<it2)return sph2[x];
if(s2[x]!=0)return s2[x];
return s2[x]=(calcS1(x)+calcS2(x>>1))%p;
}
int calcG(int x){
// cout<<"Calc G "<<x<<" "<<calcS1(x)<<" "<<calcS2(x>>1)<<"\n";
return (calcS1(x)+2*calcS2(x>>1));
}
int h[105050][40],ph[105050][40];
int res=0;
void dfs(int now,int val,int pos){
// cout<<now<<" "<<val<<" "<<calcG(n/now)<<"\n";
res+=val*calcG(n/now)%p;res%=p;
for(int i=pos;i<=tot;++i){
int x=pi[i],w=now;
if(now>n/(x*x))return ;
w*=x;
for(int j=2;j;++j){
if(w>n/x)break;
w*=x;
dfs(w,val*h[i][j],i+1);
}
}
}
signed main(){
// ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n;
init();
// for(int i=1;i<=20;i++)cout<<calcS2(i)<<" ";
// cout<<calcS1(6)<<"\n";
for(int i=1;i<=tot;i++){
h[i][0]=1;int x=pi[i],now=1;ph[i][0]=1;
if(x*x>n)break;
for(int j=1;j;++j){
if(now>n/x)break;
// cout<<now<< " "<<x<<" !\n";
int lst=now;
now*=x;ph[i][j]=now-lst;
if(j==1)continue;
h[i][j]=pi[i]^j;
// cout<<" "<<ph[i][j]<<" "<<h[i][j]<<"\n";
for(int k=0;k<j;++k){
int mul=ph[i][j-k];
if(i==1)mul*=3,mul%=p;
// cout<<" "<<x<<" "<<mul<<"\n";
h[i][j]-=h[i][k]*mul;
}
// cout<<" "<<h[i][j]<<"\n";
}
}
// cout<<"output\n";
dfs(1,1,1);res=(res%p+p)%p;
cout<<res<<"\n";
}
