积性函数求和问题的一种筛法 更新于 2026/6/16 13:51:51 作者

command_block

stO zzt Orz

本文给出了一种时间,空间复杂度均为 $O\Big(\left(\frac{N}{\log N}\right)^{2/3}\Big)$ 的积性函数求和算法。

前置知识 ： 杜教筛(+贝尔级数+powerful number)

符号约定 :

• $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筛小记。

• ③

当然，可以考虑利用朴素方法处理一部分 $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)

这是论文里面的经典题。

为了避免混淆，题目中的 $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 甚至比我原来精细实现的杜教筛还快。

#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;
}