[学习笔记]Powerful number 筛
[学习笔记]Powerful number 筛
一.前言
好耶，原来我也会 PN 。这两天整了一些筛法相关，虽然还没看 Min_25 筛，可能过两天会写吧，现在先把 PN 的笔记写了。yysy 这个东西并没有我想像中的那么难，真的搞懂之后还是挺简单的。
当然我现在所有的筛法相关都没看复杂度，😃，所以这篇是不会有复杂度的。
二.DGF与贝尔级数
这个部分不想看的可以不看，没什么大用。只是正好我一起学了就写在这里。
DGF
又称为迪利克雷生成函数。知道一般的的生成函数（OGF）是将一组形式幂级数配上自变量拿来大力和式，然后转为封闭型式。但是对于一些数论函数我们也可以将它的每个项提出来看作形式幂级数，然后大力生成函数，但是柿子有点不一样了。
然后有一个非常重要的性质：两个积性函数 \(f,g\) 对应的 DGF 的乘积为它们迪利克雷卷积 \(h=f*g\) 的 DGF，即
可以理解为一个生成函数是一个函数的性质。然后正确性自证不难。
同时我们发现一个臭名昭著的函数——黎曼函数：
从而可以用一个简单的方式表示出很多数论函数的封闭型式，我不会推导，所以可以当结论来记，部分推导过程可以看这篇
然后可以通过这个大力看关系。不过一般不常用，可以用另外一个。
贝尔级数
贝尔级数是定义在一个质数 p 上的,
然后与 DGF 相同的是，乘积和卷积一样有对应关系
自证不难。于是可以给出一个表。
这个表里的东西还是很好推的，但是由于不是这篇的重点我也就不推了。
稍加扩展，由于贝尔级数是与某一质数的幂处的函数值相关，再加上一般做题都是求积性函数，可以质因子分解然后由若干的质数幂处的取值积起来。所以贝尔级数应用要相对广泛一些，之后的 PN 和 Min_25 筛都是用的贝尔级数。
三.Powerful Number 筛
问题是给出 \(f\) 在质数幂处的表达式，且为积性函数，求前缀和。
先不说它具体是什么，开始考虑对一个积性函数求前缀和。即要求 \(\sum_{i=1}^nf(i)\) 。可以考虑杜教筛，但第一个是它不一定能筛的出来，第二个是杜教筛还是很慢的。从以前的学习中，我们学会了对一些数论函数快速的利用杜教筛求前缀和。
于是开始考虑做逆向的迪利克雷卷积，即找到一个能够快速利用杜教筛（或者其他的什么）求前缀和的积性函数 \(g\)（比如 \(\mu,\varphi\)) ，求出另一个在积性函数满足 \(h=f/g\),而这里得函数都要求在 1 次处取到 1.
先说明这样做的合理性，第一个是有几个结论:
- 积性函数的逆唯一
- 在 1 取到 1 的积性函数必有逆
- 积性函数的逆还是积性函数
第二个是可以从多项式的角度来理解：两个无穷级数做多项式除法，从低位往高位除，若是只需要答案小的有限项，那么必有一个合理解。
于是大力推式子，
然后发现后一项是求 g 的前缀和。这个由于构造可以较快求出，但是对于每一个 d 都去求就显得过于的蛋疼。于是有一个奇奇怪怪的想法：若是 h 只在一小部分有值，剩下的都是零，不就不用做那么多次了嘛。
比方说描述 h 为：只在偶数处有值，只在 5 的倍数处有值，只在质数处有值等等等等
然后发现这个很难构造。不过不知道是哪个地方的天才，引入了 powerful number 这个理念，问题就解决了。
Powerful Number
定义为所有质因子 \(p_i\) 的指数 \(e_i\) 都大于 1 的数。根据上文，描述 h 为只在 powerful number 处有值，先想在小于 n 的范围中，powerful number 的个数。
容易知道所有的 powerful number 都可以被表示成 \(a^2b^3\) 的形式，证明如下：
考虑对所有质因子 \(p_i\) 按指数 \(e_i\) 分类，可以按照模3来分，若是模 3 余 0 则分入 b，余 1 则将 4 个分入 a 剩下的分入b （因为至少为 4 ），余 2 则分 2 个到 a，剩下放入 b
这样可以知道个数不超过 \(\sum\limits_{i=1}^{\sqrt{n}}(\frac{n}{i^2})^{\frac{1}{3}}\),（相当于枚举 a，b）然后积分得到 \(\sqrt{n}\)。
然后考虑等价意义，由于是积性函数，所以由 \(h(p_i) = 0\) 可以推出只在 \(pn\) 处有值，不然就可以反证违反描述，换句话说，是积性函数且在质数处值为 0 与只 有可能 在 Powerful Number 处有值等价。
构造 h
现在已经知道了当将 h 构造为只在 PN 处有值时，只有 \(\sqrt{n}\) 个地方有值。且只需要 h 为积性函数，在质数处取值为 0 .
积性函数在一开始就已经保证（积性函数相除仍为积性函数），现在要 p 处取值为 0.考虑定义在 p 处的贝尔级数。
最上面两个贝尔级数相除需要得到第三个贝尔级数（根据迪利克雷卷积的性质），然后要求 \(h(p)=0\) 手算一下可以发现只需要 \(f(p)=g(p)\) 即可。
于是现在得到构造 h 的方法，找到一个 \(\forall p,g(p)=f(p)\) 且 \(g(1)=1\) 的函数 g 然后对于每一个质数位置大力反向迪利克雷卷积回去得到 h 在每个质数的幂处的取值。
最后深搜利用积性函数的性质算答案即可。
例题
可能以后会加吧，不知道。
Loj6053
题目给出一个积性函数 \(f\) 满足 \(f(p^k)=p \oplus k,f(1)=1\),求其前缀和。
挺裸的，观察对于奇素数 p \(f(p)=p-1\)，而对于 2 则有点不一样。按照套路直接大力构造函数 \(g=\varphi\) 然后大力跑 PN 即可。关于 2 的处理是这样的，由于在质数处并不相等，会导致对 2 的贝尔级数做除法时，\(h(2)\not=0\) ，所造成的影响是 \(h\) 有值的范围变大了，变成 \(powerful~num \or~2\cdot powerful~number\)，对复杂度的影响不算很大，就可以大力跑。
当然一般的 PN 是要对杜教筛预处理，我不是很会，另一半也是为了便于理解和调试，采用了最基础的 map 写法。
给出代码（内附自认为会比较详尽的注释）
#define fe(i,a,b) for(int i=a;i<=b;++i)
#define int long long
inline int read() {...}
const int SZ = 1e6, mod = 1e9 + 7;
int prime[SZ], phi[SZ];
bool vis[SZ];
inline void g_p(int n) {//素数筛预处理，应该都会
phi[1] = 1;
fe(i, 2, n) {
if (!vis[i])
prime[++prime[0]] = i, phi[i] = i - 1;
for (int j = 1; j <= prime[0] && i * prime[j] <= n; ++j) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) {
phi[i * prime[j]] = prime[j] * phi[i];
break;
} else
phi[i * prime[j]] = phi[prime[j]] * phi[i];
}
}
fe(i, 1, n)phi[i] = (phi[i] + phi[i - 1]) % mod;
}
map<int, int> alf;
int sh(int x) {
int y = x + 1;
if (x & 1)
y >>= 1;
else
x >>= 1;
return (x % mod) * (y % mod) % mod;//计算等差数列，注意这里的取模，非常的恶心，我调了两个小时，一定要分开取模之后再乘
}
inline int get_sum(int x) {//杜教筛模板
if (x < SZ)
return phi[x];
if (alf[x])
return alf[x];
int bet = sh(x);
for (int i = 2, j; i <= x; i = j + 1) {
j = x / (x / i);
bet = (bet - (j - i + 1) * get_sum(x / i) % mod) % mod;
}
return alf[x] = (bet % mod + mod) % mod;
}
int ans, n;
int h[SZ][66];//h就是在每一个质数幂次的取值
void dfs(int x, int y, int pos) {
ans = (ans + y * get_sum(n / x) % mod) % mod;
//计算答案，就是最开始的大力柿子，y是当前计算的对应x的h(x)
if (pos > 1 && x > n / prime[pos] / prime[pos])
return ;
//>1是对于2的特判，而在奇素数的情况下必须添加两个，判断是否超界
fe(w, pos, prime[0]) {
if (w > 1 && x > n / prime[w] / prime[w])
break ;
//判断超界，同理
int bet = x * prime[w];
for (int j = 1; bet <= n; j++, bet *= prime[w]) {//bet代表下一次的 x，j则是指数
if (!h[w][j]) {//没有计算过，严格上来讲有可能计算过为0，最好是开一个标记数组记录
int F = prime[w] ^ j, G = prime[w] - 1;//代入题目条件和构造函数
fe(i, 1, j)F = (F - G % mod * h[w][j - i]) % mod, G *= prime[w];//反向迪利克雷卷积因为只有一个值是未知的其余都是已知的
h[w][j] = F;
}
if (h[w][j])//如果为0后面就不会贡献答案了，不用算
dfs(bet, h[w][j]*y % mod, w + 1);
}
}
}
signed main() {
n = read();
g_p(SZ - 1);
fe(i, 1, prime[0])h[i][0] = 1;
dfs(1, 1, 1);
printf("%lld", (ans % mod + mod) % mod);
return 0;
}
Luogu Min_25 筛模板
我笑了，模板题总是被别的算法过掉。
还是观察条件 \(f(p)=p(p-1)=p\varphi(p)\) 于是构造函数 \(g(i)=Id(i)\cdot\varphi(i)\) 可以知道是积性函数，且在 1 处取值为 1。
直接跑 PN 就行，需要注意的是对 \(Id(i)\cdot\varphi(i)\) 杜教筛略有难度，将它和 \((1\cdot Id(i))\) 相卷，然后运用运算律得到 \((Id\cdot\varphi)*(1\cdot Id)=Id\cdot(\varphi*1)=Id_2\) 然后再结合平方和求和公式，就做出来了。
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<map>
using namespace std;
#define fe(i,a,b) for(int i=a;i<=b;++i)
#define int long long
inline int read(){
char ch=getchar();
int res=0,f=1;
for(;ch<'0'||ch>'9';ch=getchar())if(ch=='-')f=-1;
for(;ch>='0'&&ch<='9';ch=getchar())res=(res<<3)+(res<<1)+(ch-'0');
return res*f;
}
const int SZ=1e6+10,mod=1e9+7;
int prime[SZ],phii[SZ];
bool vis[SZ];
inline void g_p(int n){//筛
phii[1]=1;
fe(i,2,n){
if(!vis[i])prime[++prime[0]]=i,phii[i]=i-1;
for(int j=1;j<=prime[0]&&i*prime[j]<=n;++j){
vis[i*prime[j]]=1;
if(i%prime[j]==0){
phii[i*prime[j]]=phii[i]*prime[j];
break;
}
phii[i*prime[j]]=phii[i]*phii[prime[j]];
}
}
fe(i,1,n)phii[i]=(phii[i-1]+phii[i]*i%mod)%mod;//这里再乘上Id
}
int n;
inline int SH(int x){//平方和
int y=x+1,z=2*x+1;
if(x%2==0)x>>=1;
else if(y%2==0)y>>=1;
else z>>=1;
if(x%3==0)x/=3;
else if(y%3==0)y/=3;
else z/=3;
return (x%mod)*(y%mod)%mod*(z%mod)%mod;//坑点
}
inline int SG(int x){//一次和
int y=x+1;
if(y%2==0)y>>=1;
else x>>=1;
return (x%mod)*(y%mod)%mod;
}
map<int,int> alf;
inline int get_sum(int x){//嘟嘟嘟嘟嘟教筛
if(x<SZ)return phii[x];
else if(alf[x])return alf[x];
int bet=SH(x);
for(int l=2,r;l<=x;l=r+1){
r=x/(x/l);
bet=(bet-(SG(r)-SG(l-1))%mod*get_sum(x/l)%mod)%mod;
}
return alf[x]=bet;
}
int ans,h[SZ][66];
void dfs(int x,int y,int pos){//基本一模一样，只是没有关于2的特判了
ans=(ans+y*get_sum(n/x)%mod)%mod;
if(x>n/prime[pos]/prime[pos])return ;
fe(w,pos,prime[0]){
if(x>n/prime[w]/prime[w])return ;
int bet=x*prime[w];
for(int j=1;bet<=n;j++,bet*=prime[w]){
if(!h[w][j]){
int F=((bet/x)%mod)*((bet/x-1)%mod)%mod,G=(prime[w]%mod)*((prime[w]-1)%mod)%mod;
fe(i,1,j)F=(F-G*h[w][j-i]%mod)%mod,G=G*(prime[w]%mod)%mod*(prime[w]%mod)%mod;
h[w][j]=F;
}
if(h[w][j])dfs(bet,y*h[w][j]%mod,w+1);
}
}
}
signed main(){
g_p(SZ-1);
n=read();
fe(i,1,prime[0])h[i][0]=1;
dfs(1,1,1);
printf("%lld",(ans%mod+mod)%mod);
return 0;
}
