0%

Min_25筛

一种亚线性复杂度求积性函数前缀和的筛法。

适用条件:

  • \(f(p)\) 为多项式(除代码外,本文 \(p\) 均指质数)

  • \(f(p^k)\) 便于计算

我们以 LuoguP5325 为例简单介绍。

\(\sum\limits_{i=1}^nf(i)\),其中 \(f\) 是一个积性函数,且 \(f(p^k)=p^k(p^k-1)\)\(n\leqslant10^{10}\)

Min_25筛的主要思想是把要求的前缀和转化为质数和非质数两部分进行计算。

第一部分

对于质数部分,我们定义以下函数: \[ g(n,j)=\sum\limits_{i=2}^n[i\ \text{is a prime or}\ \text{minp}(i)>p_j]i^k \] 其中,\(p_j\) 表示第 \(j\) 个质数,\(\text{minp}(i)\) 表示 \(i\) 的最小质因子。翻译成中文就是对所有质数和最小质因子大于第 \(j\) 个质数的 \(k\) 次幂求和。

这里需要注意,如果我们求的 \(f(p)\) 是一个多项式,就拆成单项式。比如这题就是 \(p^2-p\),那么我们就计算 \(g_1\)\(g_2\) 分别代表一次项的和以及二次项的和,这里先不管系数,会在后面算上。

我们考虑计算 \(g(n,j)\),这是一个比较明显的可以递推的式子,于是我们尝试从 \(g(n,j-1)\) 递推到 \(g(n,j)\)

发现 \(g(n,j-1)\)\(g(n,j)\) 实际上只差了最小质因子为 \(p_j\)(不包括 \(p_j\) 本身)的函数值,所以我们只需要减掉这些多出来的值就可以了,下面先给出式子: \[ g(n,j)=g(n,j-1)-p_j^k\left(g(\frac{n}{p_j},j-1)-g(p_{j-1},j-1)\right) \] 注意一下这里除法为了美观没有写取整。我们发现 \(i^k\) 是一个完全积性函数,所以对于所有最小质因子为 \(p_j\) 的数,我们都可以提一个 \(p_j\) 出来而且不需要管是否互质(也就是剩下的仍有 \(p_j\) 这个因子),然后剩下的数一定小于等于 \(\frac{n}{p_j}\) 且最小质因子仍大于等于 \(p_j\),所以就是 \(g(\frac{n}{p_j},j-1)\)。然后发现所有小于 \(p_j\) 的质数都被多减掉了,所以再用 \(g(p_{j-1},j-1)\),即所有小于等于 \(p_{j-1}\) 的质数的幂次和。

我们注意到小于等于 \(n\) 的合数,它的最小质因子一定小于等于 \(\sqrt n\),即如果 \(p_j>\sqrt n\) 那么 \(g(n,j)=g(n,j-1)\),所以我们只需要枚举小于等于 \(\sqrt n\) 的质数。然后我们发现 \(g(p_{j-1},j-1)\) 求的就是所有小于等于 \(p_{j-1}\) 的质数的幂次和,而且 \(p_{j-1}\) 是小于等于 \(\sqrt n\) 的。所以我们一开始线性筛到 \(\sqrt n\) 求出每个质数的这个东西,我们设它为 \(sp(j-1)\)。然后对于 \(1\)\(n\)\(k\) 次方和,也即是 \(g(n,\text{maxp}(n))\),其中 \(p_{\text{maxp}(n)}\) 是小于等于 \(n\) 的最大质数。

但是这样枚举量还是太大,我们发现一个很重要的式子:\(\lfloor\frac{\lfloor\frac{n}{a}\rfloor}{b}\rfloor=\lfloor\frac{n}{ab}\rfloor\)。所以所有能用到的第一个数值都形如 \(\lfloor\frac{n}{a}\rfloor\),而众所周知这样的数只有 \(O(\sqrt n)\) 个,所以我们只需要求出这样的 \(g\) 就可以了。具体实现因为下标还是很大,用 mapunordered_map 都比较慢。根据那个 \(O(\sqrt n)\) 的证明,我们手动离散化一下:记 id1[x]\(x\) 离散化的下标,id2[x]\(\lfloor\frac{n}{x}\rfloor\) 的下标,可以看代码。

这一部分的复杂度是 \(O(\frac{n^{\frac{3}{4}}}{\log n})\) 的。

第二部分

接下来,我们计算所有数的和,记 \(s(n,j)\)\(\sum\limits_{i=2}^n[\text{minp}(i)>j]f(i)\),即所有最小质因子大于 \(p_j\)\(f\) 的和。

我们依旧尝试递推,首先,我们把所有大于 质数的函数和 \(\sum g(n,\text{maxp}(n))-sp(j)\) 算上(这里带上系数,比如板子题就是 \(g2(n,\text{maxp}(n))-sp2(j)-(g1(n,\text{maxp}(n))-sp1(j))\) ,那么剩下合数部分枚举提出一个最小质因子,先看式子: \[ s(n,j)=\sum(g(n,\text{maxp}(n))-sp(j))+\sum\limits_{c>j,p_c^e\leqslant n}f(p_c^e)\left(s(\frac{n}{p_c^e},c)+[e\not=1]\right) \] 因为这一次并不是完全积性函数了,所以我们要一次把最小质因子提全,这样的话剩下的最小质因子都大于 \(p_c\),然后因为我们没算 \(1\),所以当 \(e\not=1\) 的时候 \(f(p_c^e)\) 没有算进去,当 \(e=1\) 时会在前面质数部分算到。

\(s(n,0)+f(1)\) 即为答案。这里直接递归计算,不需要记忆化。

这一部分的复杂度据说是 \(O(n^{1-\epsilon})\) 的,即大于任何一个 \(O(n^k),k<1\),但是常数很小。

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll const p=1e9+7,inv2=500000004,inv6=166666668;
ll n,pr[100005],cnt,sp1[100005],sp2[100005],tot,bas[300005],g1[300005],
g2[300005],sqr;
int id1[100005],id2[100005];
bool flag[100005];
inline ll mod(ll const &x){return x>=p?x-p:x;}
ll s(ll i,int j)
{
if(pr[j]>=i)return 0;
int t=(i<=sqr)?id1[i]:id2[n/i];
if(sp1[j]<0||sp2[j]<0)exit(-1);
ll res=(g2[t]-sp2[j]-(g1[t]-sp1[j])+p+p)%p;//part1:质数部分
for(int c=j+1;c<=cnt&&pr[c]*pr[c]<=i;c++)//part2 注意这里依旧只需要消去 sqrti 的最小质因子
{
ll now=pr[c];
for(int e=1;now<=i;e++,now*=pr[c])
{
ll tmp=now%p;
res=(res+tmp*(tmp-1+p)%p*(s(i/now,c)+(int)(e!=1)))%p;
}
}
return res;
}
int main()
{
scanf("%lld",&n);
sqr=sqrt(n);
for(ll i=2;i<=sqr;i++)//get prime
{
if(!flag[i])pr[++cnt]=i;
for(int j=1;j<=cnt&&i*pr[j]<=sqr;j++)
{
flag[i*pr[j]]=1;
if(!(i%pr[j]))break;
}
}
for(int i=1;i<=cnt;i++)//计算sp
sp1[i]=mod(sp1[i-1]+pr[i]),
sp2[i]=mod(sp2[i-1]+pr[i]*pr[i]%p);
for(ll l=1,r;l<=n;l=r+1)
{//计算 g(x,0),因为 g(n,j) 总是由 g(xxx,j-1) 推出而且后面只需要 g(n,maxp(n)) 所以滚动掉一维
r=n/(n/l);
bas[++tot]=n/l;//这一次计算 g(n/l,0) 也就是 2 到 n/l 的幂次和
g1[tot]=bas[tot]%p;//因为减少取模所以暂时借用一下这个数
g2[tot]=(g1[tot]*(g1[tot]+1)%p*(2*g1[tot]+1)%p*inv6+p-1)%p;
g1[tot]=mod((g1[tot]+1)*g1[tot]%p*inv2%p+p-1);
if(bas[tot]<=sqr)id1[bas[tot]]=tot;
else id2[l]=tot;//此时一定有 l=r
}
for(int j=1;j<=cnt;j++)//为了和写的式子一样换了一下内外层变量
for(int i=1;i<=tot&&pr[j]*pr[j]<=bas[i];i++)//递推 g(bas[i],j)
{
ll t=bas[i]/pr[j];
t=(t<=sqr)?id1[t]:id2[n/t];//找到 bas[i]/pr[j] 对应的下标
g1[i]=mod(g1[i]-pr[j]*(g1[t]-sp1[j-1]+p)%p+p);
g2[i]=mod(g2[i]-pr[j]*pr[j]%p*(g2[t]-sp2[j-1]+p)%p+p);
}
printf("%lld",mod(s(n,0)+1));
return 0;
}