PN筛入门

PN 筛

Powerful Number 筛

oi-wiki 上说是杜教筛的 expand，但是这玩意 只能做积性函数前缀和

若未加说明，以下涉及的 $f,g,h$ 均为积性函数。$p$ 均为质数

Powerful number

定义一个数 $N=\prod _{i=1}^m{p}_i^{c_i}$ 为 powerful number 当且仅当 $\mathrm{\forall }i,c_i\ge 2$。特别的，在本文中我们认为 $1$ 也是powerful number。

性质：

• $\mathrm{\forall }N\in \mathbb{PN},\mathrm{\exists }a,b\in \mathbb{Z},N=a^2{b}^3$

  证明容易，对于每个 $c_{i}mod2=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}}\ast n^{\frac{1}{2}·\frac{1}{3}})=O(\sqrt{n})$​

• $1\sim n$ 的 $PN$ 每个质因子都不大于 $\sqrt{n}$

算法内容

考虑求出 $S(n)=\sum _{i=1}^{n}f(i)$。

我们考虑构造一个积性函数 $g$。（这里与杜教筛的区别在于构造难易程度上）

满足以下要求：

• $g$ 易求前缀和
• $\mathrm{\forall }p,g(p)=f(p)$

再构造 $h\ast g=f$。根据定义式有：

$$\begin{array}{rl}& f(n)=\sum _{d|n}h(d)g(\frac{n}{d})\\ ⟹& h(n)=f(n)-\sum _{d|n,d<n}h(d)g(\frac{n}{d})\end{array}$$

函数 $h$ 有以下性质：

• $h(p)=0$，这是因为 $f(p)=g(1)h(p)+g(p)h(1),h(1)=g(1)=1,g(p)=h(p)⟹h(p)=0$

• $\mathrm{\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}^{n}g(i)$，可以快速计算 $G(n)$ 。

则有：

$$\begin{array}{rl}S(n)& =\sum _{i=1}^{n}f(i)\\ & =\sum _{i=1}^n\sum _{d|i}g(d)h(\frac{i}{d})\\ & =\sum _{i=1}^{n}h(i)\sum _{t=1}^{it\le n}g(t)\\ & =\sum _{i=1}^{n}h(i)G(\lfloor \frac{n}{i}\rfloor )\end{array}$$

这告诉我们，若可以快速求得 $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}^{n}f(i)$。

显然构造 $g(n)=n·\phi (n)$，可以在质数处"拟合"。

现在考虑求出 $G(n)$。

求解 $G(n)$

看到 $\phi$ 可以想到设计 $g^{\prime }(n)=n$。

所以 $(g\ast g^{\prime })(n)=\sum _{d|n}g(d)\frac{n}{d}=\sum _{d|n}n\phi (d)=n^2$。

带入杜教筛的式子，有：

$$\begin{array}{rl}{g}^{\prime }(1)G(n)=G(n)& =\sum _{i=1}^n(g\ast g^{\prime })(i)-\sum _{i=2}{g}^{\prime }(i)G(\lfloor \frac{n}{i}\rfloor )\\ & =\sum _{i=1}^n{i}^2-\sum _{i=2}^{n}i·G(\lfloor \frac{n}{i}\rfloor )\end{array}$$

问题至此解决，下附一份板子。

#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也可以综合利用已知信息

给定函数：

$$\begin{array}{rl}f(1)& =1\\ f(p^c)& =p\oplus c\\ f(ab)& =f(a)f(b),gcd(a,b)=1\end{array}$$

求 $S(n)=\sum _{i=1}^{n}f(i)$

考虑PN筛构造。

注意到 $f(p)=\phi (p),p>2,f(2)=3$。

所以可以构造 $g(n)=\phi (n),nmod2=1,g(n)=3\phi (n),nmod2=0$。

对于一个积性函数，强行钦定含某个质因子的数多一个带权系数，仍然是积性函数

因此 $h$ 容易求出。

现在考虑求解 $G(n)$。

$$\begin{array}{rl}G(n)& =\sum _{i=1}^n\phi (i)+2\sum _{i=1}^{\lfloor \frac{n}{2}\rfloor }\phi (2i)\\ & =S_{\phi }(n)+S_{\phi (2)}(\lfloor \frac{n}{2}\rfloor )\end{array}$$

求解 $S_{\phi }(n)$ 是基本功，根据杜教筛可以知道 $S_{\phi }(n)=\frac{n(n+1)}{2}-\sum _{i=2}^{n}S(\lfloor \frac{n}{i}\rfloor )$

考虑 $S_{\phi (2)}(n)$

$$\begin{array}{rl}{S}_{\phi (2)}& =\sum _{i=1}^n\phi (2i)\\ & =\sum _{i=1}^n\phi (i)(1+[imod2=0])\\ & =\sum _{i=1}^n\phi (i)+\sum _{i=1}^{\lfloor \frac{n}{2}\rfloor }\phi (2i)\\ & =S_{\phi }(n)+S_{\phi (2)}(\lfloor \frac{n}{2}\rfloor )\end{array}$$

所以可以递推求解 $S_{\phi (2)}$。

而且因为求解 $S_{\phi }(n)$ 的时候会求解 $S_{\phi (i)}(\lfloor \frac{n}{2}\rfloor )$。所以 $S_{\phi (2)}$ 直接记搜的复杂度与杜教筛求解 $S_{\phi }$ 的复杂度相同，都是 $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";
}