Powerful Number 筛

PN 筛。有时候确实比 Min_25 筛快。

首先定义 Powerful Number 为形如 \(a^2b^3\) 的数。一个结论是 \(n\) 以内的 Powerful Number 是 \(O(\sqrt n)\) 的。证明枚举 \(a\) 对 \(b\) 乱积分就可以。

假如说我们现在要求 \(S(n)=\sum_{i=1}^nf(i)\) 。那么我们构造一个易求前缀和的积性函数 \(g\) ，使其在质数处取值和 \(f\) 相同。

然后构造函数 \(h\) ，满足 \(f=g*h\)，\(*\) 为 Dirichlet 卷积。显然 \(h\) 只在 Powerful Number 处可能有值。

那么推式子：

\[\begin{aligned} S(n)&=\sum_{i=1}^nf(i)\\ &=\sum_{i=1}^n\sum_{d|i}h(d)g(\frac id)\\ &=\sum_{d=1}^n\sum_{i=1}^{\lfloor \frac nd \rfloor}h(d)g(i)\\ &=\sum_{d=1}^nh(d)\sum_{i=1}^{\lfloor \frac nd\rfloor}g(i) \end{aligned} \]

那么 \(g\) 的前缀和容易求得，接下来就只需要找到所有 \(n\) 以内的 Powerful Number 并找到 \(h\) 在其中的取值。

找 Powerful Number 通常采用 dfs 搜索各个素数的指数的方式。然后关于求 \(h(p^k)\) ，可以求出它的通项公式，也可以根据 \(f(p^k)=\sum_{i=1}^kg^{p^i}h^{p^{k-i}}\) 递推得到。

复杂度分析：这里有三步。

1. 找 Powerful Number：显然是 \(O(\sqrt n)\)。
2. 求 \(g\) 的前缀和：这个因函数而异。有时候需要采用一些其他的筛法。
3. 求 \(h(p^k)\)：质数个数 \(O(\dfrac n{\log n})\)，每个质数需要 \(O(\log^2n)\) 推出所有 \(h(p^k)\) ，总复杂度 \(O(\sqrt n\log n)\)。

空间复杂度显然是 \(O(\sqrt n)\)。由于它的时间复杂度瓶颈大多数时候在杜教筛，所以导致它在不需要杜教筛的时候有着优秀的复杂度，可以处理 \(10^{13}\) 级别的问题。

下面上个洛谷 Min_25 筛板子的代码。构造 \(g(x)=x\varphi(x)\) 即可。

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cmath>
#define int long long
using namespace std;
const int mod=1000000007;
const int inv2=(mod+1)>>1,inv6=166666668;
int n,sq,ans,phi[5000010],p[5000010],g1[5000010],g2[5000010];
bool v[5000010];
void get(int n){
    phi[1]=1;
    for(int i=2;i<=n;i++){
        if(!v[i])p[++p[0]]=i,phi[i]=i-1;
        for(int j=1;j<=p[0]&&i*p[j]<=n;j++){
            v[i*p[j]]=true;
            if(i%p[j]==0){
                phi[i*p[j]]=phi[i]*p[j];break;
            }
            phi[i*p[j]]=phi[i]*phi[p[j]];
        }
    }
    for(int i=1;i<=n;i++)g1[i]=(g1[i-1]+i*phi[i])%mod;
}
int s(int n){
    n%=mod;
    return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
int getf(int x){
    if(x<=sq)return g1[x];
    if(g2[n/x])return g2[n/x];
    int ans=s(x);
    for(int l=2,r;l<=x;l=r+1){
        r=x/(x/l);
        ans=(ans-(l+r)%mod*(r-l+1)%mod*inv2%mod*getf(x/l)%mod+mod)%mod;
    }
    return g2[n/x]=ans;
}
void dfs(int k,int m,int h){
    if(k>p[0]||m*p[k]>n){
        if(n/m<=sq)ans=(ans+h*g1[n/m])%mod;
        else ans=(ans+h*g2[n/(n/m)])%mod;
        return;
    }
    dfs(k+1,m,h);
    for(int j=p[k]*p[k],e=2;m*j<=n;j*=p[k],e++){
        dfs(k+1,m*j,j%mod*(p[k]-1)%mod*(e-1)%mod*h%mod);
    }
}
signed main(){
    scanf("%lld",&n);
    sq=pow(n,2.0/3);
    get(sq);getf(n);
    dfs(1,1,1);
    printf("%lld\n",ans);
}