积性函数求和问题的一种筛法

command_block · 2020-01-26 18:43:14 · 个人记录

stO zzt Orz

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

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

符号约定 :

• 涉及的积性函数均有 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) ,满足如下性质:

  给出范围 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 的递推式 ：

边界 : 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)

直接使用上述算法即可。

这东西常数较小 (远小于洲阁筛,同样精细实现的话比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;
}