积性函数求和问题的一种筛法
command_block
2020-01-26 18:43:14
stO zzt Orz
本文给出了一种时间,空间复杂度均为 $O\Big(\left(\frac{N}{\log N}\right)^{2/3}\Big)$ 的积性函数求和算法。
前置知识 ： [杜教筛(+贝尔级数+powerful number)](https://www.luogu.com.cn/blog/command-block/du-jiao-shai)
**符号约定** :
- $S_F(n)=\sum\limits_{i=1}^nF(i)$ ,即前缀和。
- $p$ 泛指质数。
- 涉及的积性函数均有 $F(1)=1$。
- 对函数 $F$ ,给出**范围** $N$ ,对所有 $m=\lfloor N/i\rfloor$ 求解 $S_F(m)$ ,这种操作称之为**块筛**,有时省略范围。
大 $N$ 一般指题目中给出的 $N$，而小 $n$ 则用作分析。
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
------------
# 1.质数的 $c$ 次幂前缀和（优化）
前半部分见 [Min_25筛小记](https://www.luogu.com.cn/blog/command-block/min25-shai-xiao-ji)。
- ③
当然，可以考虑利用朴素方法处理一部分 $h$ 的前缀，假设我们预处理 $[1,T]$ ,即 $n\leq T$ 的所有 $h(T,\_)$。
每次新加质数的时候，会将 $T$ 以内暴力筛除一遍,需要**树状数组**维护动态前缀和，由于每个数只会被筛除一次，复杂度为 $O(T\log T)$。
递推时次数为 $O\bigg(\sum\limits_{d=1}^{N/T}\frac{\sqrt{N/d}}{logN}\bigg)=O\Big(\frac{N}{logN\sqrt{T}}\Big)$ ,不过注意到树状数组需要一个 $O(\log T)$ ,这部分复杂度为 $O(\frac{N}{\sqrt{T}})$。
总的复杂度就是 $O(\frac{N}{\sqrt{T}}+TlogT)$ ,取 $T=O((\frac{N}{logN})^{2/3})$ 可得复杂度为 $O(N^{2/3}log^{1/3}N)$ 。
- ④
考虑省去递推运算中的部分树状数组查询操作。
注意到 $h(n,k)$ 中如果 $p_k>\sqrt{n}$ 的话函数值不再改变，可以直接使用静态前缀和。
此时, $h(\lfloor n/p_k\rfloor,k-1)$ 仍在变化的条件是 $\sqrt{n/p_k}\geq p_k\rightarrow n\geq p_k^3$。
也就是说,在第 $k$ 次转移中，求 $h(n,k)$ 的查询需要树状数组的条件为 $p_k^3\leq n$。
对于每个 $h(n,\_)$ 有 $O(\frac{\sqrt[3]{n}}{\log n})$ 个树状数组查询操作,复杂度即为$O(\sqrt[3]{n})$。
这部分总复杂度是 $O\bigg(\sum\limits_{d=1}^{N/T}\sqrt[3]{N/d}\bigg)=O(\frac{N}{T^{2/3}})<O(\frac{N}{logN\sqrt{T}})$ ,不是瓶颈。
这样子的总复杂度是 $O\Big(\frac{N}{logN\sqrt{T}}+TlogT\Big)$ 取 $T=O(\frac{N^{2/3}}{log^{4/3}N})$可得复杂度为 $O(\frac{N^{2/3}}{log^{1/3}N})$。
- ⑤
考虑在 $p_k\leq T_2$ 时暴力递推，这部分复杂度为 $O(T_2\sqrt{N}/\log N)$。
此时还需筛去的数满足最小素因子大于 $T_2$ ，且是合数。
笔者在 $T\leq 10^8$ 内验证发现，使用 $O(T/\sqrt{T_2})$ 拟合效果不错。
参考 ： 当 $T_2=50,T=5\times 10^6$ 时，树状数组修改约为 $7\times 10^5$ 次。（对应 $n≈10^{11}$）
取 $T_2=N^{1/6}$ ，则复杂度为 $O\Big(N^{2/3}/\log N+\frac{N}{logN\sqrt{T}}+T\Big)$
令 $T=\left(\frac{N}{\log N}\right)^{2/3}$ 可得复杂度为 $O\Big(\left(\frac{N}{\log N}\right)^{2/3}\Big)$。
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
------------
# 2.积性函数求和
- **问题概述** :
给出一个**积性函数** $F(n)$ ,满足如下性质:
- $F(p)$ 为关于 $p$ 的简单多项式。
- $F(p^k)$ 能够较快速求出。
给出范围 $N$ ,求出 $F$ 的块筛。
根据 $\rm powerful\ number$ 理论，构造数论函数 $G(p)=F(p)$ ，只需要求出 $G$ 的块筛，就可以在 $O(n^{0.6})$ 内得到 $F$ 的块筛。
不妨令 $G(p^k)=0\ (k>1)$。
我们以 $k$ **从大到小**的方式,考虑“**最小素因子**为 $p_k$ 的数”。
设 $S_{n,k}$ 为最小素因子 $\geq k$ ，且素因子次数最高一次的数的集合。
（若素因子次数为 $2$ ，则 $G$ 函数值为 $0$ ，可以不考虑）
仿照前面 $h$ 的求解方法,设 $r(n,k)=\sum\limits_{i\in S_{n,k}}^nG(i)$。
设 $D_{n,k}$ 为最小素因子 $=k$ ，且素因子次数最高一次的数的集合。
按照定义有 $S_{n,k}=S_{n,k+1}+D_{n,k}$。
可以发现 $D_{n,k}=p_kS_{\lfloor n/p_k\rfloor,k+1}$。
则 $S_{n,k}=S_{n,k+1}+p_kS_{\lfloor n/p_k\rfloor,k+1}$
可以据此写出 $r$ 的递推式 ：
$$r(n,k)=r(n,k+1)+G(p_k)*\begin{cases}r(\lfloor n/p_k\rfloor,k+1)&(p_k\leq\sqrt{n})\\1&(\sqrt{n}<p_k\leq n)\\0&(n<p_k)\end{cases}$$
边界 : $s(n,|P|)=1$ ,此时只有 $G(1)=1$ 有贡献。
- **解释** :
当 $p_k\leq\sqrt{n}$ 时,由于不会出现多次质因子，利用积性显然。
当 $p_k>\sqrt{n}$ 时, $\lfloor n/p_k\rfloor<p_k$ ,只有 $G(1)=1$ 可能贡献。
当 $p_k>n$ 时, $\lfloor n/p_k\rfloor=0$ ,无贡献。
------------
现在看来,我们似乎需要枚举所有的质数。
注意到当 $p_k>\sqrt{n}$ 时,只有质数处的 $G(p)$ 会被统计到。
当 $p_{t}$ 为第一个大于 $\sqrt{n}$ 的质数时， $s(n,t)=\sum\limits_{\sqrt{n}<p\leq n}G(p)$
$G(p)$ 能拆成若干个单项式分别求和,那么就转化成了 $2$ 中质数的 $c$ 次幂前缀和问题。
现在我们能直接从 $s(n,t)$ 开始往前递推,只需要考虑 $\sqrt{n}$ 以内的质数。
$s$ 的递推式和 $h$ 是相似的,可以使用相同的解法,这里不再赘述了。
不过注意对于 $s(n,k)$ ,当 $p_k>\sqrt{n}$ 的时候,并不是和 $h$ 一样不变。这可以用时间戳解决。
- **例题** : [SP34096 DIVCNTK - Counting Divisors (general)](https://www.luogu.com.cn/problem/SP34096)
这是论文里面的经典题。
为了避免混淆，题目中的 $K$ 使用大写,为了方便我们令 `K++`。
**题意** : 求 $\sum\limits_{i=1}^nd(i^{K-1})$。
在这道题里面,设 $F(n)=d(n^{K-1})$ ,易得这也是个积性函数。
可得 $F(p)=d(p^{K-1})=K;F(p^k)=d(p^{k{K-1}})=(K-1)k+1$
直接使用上述算法即可。
这东西常数较小 (远小于洲阁筛,同样精细实现的话比Min_25不相上下)。
在 DIVCNTK 拿到了 LG·Rk1 ,原站 SPOJ·Rk5。
把 $k$ 改成 2,3 均能通过 DIVCNT2 + DIVCNT3 ,而且踩了一堆杜教筛。
目前在 DIVCNT3 甚至比我原来精细实现的杜教筛还快。
```cpp
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<cmath>
#define ull unsigned long long
#define lbit(x) ((x)&-(x))
#define Limit 805000
using namespace std;
int tn,lim,lp,p[Limit>>3],c[Limit],e[Limit];
void getsth()
{
for (int i=2,t;i<=lim;i++){
if (!e[i])c[p[++tn]=i]=1;
for (int j=1;j<=tn&&(t=p[j]*i)<=lim;j++){
e[t]=p[j];
if (i%p[j]==0)break;
if(c[i])c[t]=c[i]+1;
}
}p[tn+1]=Limit;
}
ull N,K,tH[35];
void getPN()
{
tH[0]=1;
for (int k=2;k<=33;k++)
tH[k]=(K-1)*k+1-K*tH[k-1];
}
int tot;
ull t[Limit],x[Limit],tag;
inline void upd(int p,ull x){
while(p<=lim){t[p]+=x;p+=lbit(p);}
}
inline ull qry(int p){
ull sum=0;
while(p){sum+=t[p];p^=lbit(p);}
return sum;
}
void getP()
{
for (tot=1;N>(ull)lim*tot;tot++)x[tot]=N/tot-1;
tot--;
for (int i=1;i<=lim;i++)t[i]=lbit(i);
upd(1,-1);
for (int tp=1;tp<=lp;tp++){
int P=p[tp],u=min(N/((ull)P*P),(ull)tot),
uu=min(u,tot/P);
ull sav=qry(P-1),t=N/P;
for (int i=1;i<=uu;i++)
x[i]-=x[P*i]-sav;
if (t<(1ull<<31))
for (int tt=t,i=uu+1;i<=u;i++)
x[i]-=qry(tt/i)-sav;
else
for (int i=uu+1;i<=u;i++)
x[i]-=qry(t/i)-sav;
if ((ull)P*P<=lim)
for (int j=P*P;j<=lim;j+=P)
if (e[j]==P)upd(j,-1);
}
}
void sieve()
{
for (lp=1;(ull)p[lp]*p[lp]<=N;lp++);
getP();
memset(t,0,sizeof(ull)*(lim+5));
for (int i=1;i<=tot;i++)x[i]=K*(x[i]-lp);
for (int i=lp+1;p[i]<=lim;i++)upd(p[i],K);
ull tK[35],tag=0;tK[0]=1;
for (int k=1;k<=33;k++)tK[k]=tK[k-1]*K;
for (int i=lp;i;i--){
int P=p[i],u=min(N/((ull)P*P),(ull)tot),
uu=min(u,tot/P);
for (int i=1;i<=uu;i++)
x[i]+=K*(x[P*i]+tag);
ull t=N/P;
if (t<(1ull<<31))
for (int tt=t,i=uu+1;i<=u;i++)
x[i]+=K*qry(tt/i);
else
for (int i=uu+1;i<=u;i++)
x[i]+=K*qry(t/i);
upd(P,K);tag+=K;
if ((ull)P*P<=lim)
for (int j=P*P;j<=lim;j+=P)
if (e[j]==P&&c[j])upd(j,tK[c[j]]);
}for (int i=1;i<=tot;i++)x[i]+=1+tag;
upd(1,1);lim=sqrtl(N);
for (int i=tot+1;i<=lim;i++)x[i]=qry(N/i);
tot=lim;
}
inline ull S(ull d){
if (d<=tot)return x[d];
return qry(N/d);
}
ull ans;
void dfs(ull n,ull num,int t)
{
ans+=num*S(n);
for (int i=t,k;i<=lp;i++){
if ((ull)p[i]*p[i]>N/n)return ;
k=2;
for (ull x=(ull)p[i]*p[i];x*n<=N;k++,x*=p[i])
dfs(n*x,num*tH[k],i+1);
}
}
#define gL(N) max((int)max(1.75*pow(N/log2(N),0.6666),sqrt(N)+5),100)
void solve()
{
lim=gL(N);if (lim>N)lim=N;
getPN();sieve();
ans=0;dfs(1,1,1);
printf("%llu\n",ans);
}
ull Sn[10050],Sk[10050],mN;
int T;
int main()
{
scanf("%d",&T);
for (int i=1;i<=T;i++){
scanf("%llu%llu",&Sn[i],&Sk[i]);
Sk[i]++;mN=max(mN,Sn[i]);
}lim=gL(mN);if (lim>mN)lim=mN;
getsth();
for (int i=1;i<=T;i++){
N=Sn[i];K=Sk[i];
solve();
}return 0;
}
```