Powerful number 筛略解
Powerful number 筛可以用来求出一类积性函数的前缀和，最快可以达到根号复杂度。
- 以下字母大都代表整数；
- 以下 \(p\) 代表质数；
- 以下 \(*\) 代表狄利克雷卷积；
- 以下 \(a/b=\lfloor{a\over b}\rfloor\)；
以下任何不够严谨的地方都请忽略。
定义
定义一个 Powerful number 为任意一个质因子的次数都不小于 2 的数，如 \(72=2^3\times 3^2\)，\(4000=2^5\times5^3\)。Powerful number 有如下性质：
- 一个 Powerful number 一定可以表示为 \(a^2b^3\) 的形式，如 \(72=3^2\times2^3\)，\(4000=2^{(2+3)}\times5^3=2^2\times(2\times5)^3=2^2\times10^3\)；
- \(n\) 以内的 Powerful number 个数是 \(O(\sqrt n)\) 的：\[\int_1^{\sqrt n}\sqrt[3]{n\over 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)=p\oplus q\)。（LOJ6053 简单的函数）
显然 \(f(p)=\begin{cases}p-1\quad&(p\not=2)\\p+1&(p=2)\end{cases}\)。
构造 \(g(n)=\begin{cases}\varphi(n)\quad&(2\nmid n)\\3\varphi(n)\quad&(2\mid n)\end{cases}\)。
下面主要解决两个问题：\(g^{-1}\) 与 \(g\) 的前缀和。
用类似于例 1 的方法推导一番可得：\(q\) 足够大时，\(g*g^{-1}(p^q)=\begin{cases} g^{-1}(p^q)-g^{-1}(p^{q-1})+p\cdot g*g^{-1}(p^{q-1})\quad(p\not =2)\\ g^{-1}(p^q)+g^{-1}(p^{q-1})+p\cdot g*g^{-1}(p^{q-1})\quad(p=2)\end{cases}\)。
\(g^{-1}(p^q)=\begin{cases} g^{-1}(p^{q-1})\quad&(p\not =2)\\ -g^{-1}(p^{q-1})\quad&(p=2)\end{cases}\)。
所以：
- \(g^{-1}(1)=1\)
- \(g^{-1}(p^q)=\begin{cases}-p+1\quad&(p\not =2)\\(-1)^q(p+1)\quad&(p=2)\end{cases}\quad(q>0)\)
然后推导 \(g\) 的前缀和：\(\sum\limits_{i=1}^{n}g(i)=\sum\limits_{i=1}^n\varphi(i)+2\sum\limits_{i=1}^{n/2}\varphi(2i)\)
设 \(S'(n)=\sum\limits_{i=1}^{n}\varphi(2i)\)。则当 \(2\mid n\) 时：
假如 \(2\nmid n\)，\(S'(n)=S'((n-1)/2)+\sum\limits_{i=1}^{n-1}\varphi(i)+\varphi(2n)=S'(n/2)+\sum\limits_{i=1}^n\varphi(i)\)（完 全 一 致）。
因此在杜教筛出 \(O(\sqrt n)\) 个点的 \(\varphi\) 的前缀和后，再处理出 \(O(\sqrt n)\) 个 \(S'\) 的值。\(\sum\limits_{i=1}^{n}g(i)=\sum\limits_{i=1}^n\varphi(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）
