P4213

· · 题解

杜教筛（给定积性函数 f\ast g=h，已知 g,h 块筛，求 f 的块筛）可以不借助 zak 筛手法做到 O(\frac{n^\frac23}{\log n})。

我们考虑人为地让取值变稀疏。有结论：对于任意 0<\alpha<1，n 以内 \operatorname{lpf}(n)\geq n^\alpha 的数的个数是 O(\frac n{\log n}) 的，或者说是 n 以内 n^\alpha-rough 数的个数，证明可以在此处找到。我们构造积性函数 f'(n)=[\operatorname{lpf}(n)\geq n^\alpha]，对 g',h' 同理，则仍有 f'\ast g'=h'。

若 g,h 是完全积性函数，我们可以用 Min_25 筛的第一部分筛除 n^\alpha 以内的质因数得到 g',h'，这是 \sqrt n\pi(n^\alpha)=O(\frac{n^{\frac12+\alpha}}{\log n}) 的。我们也可以用 Min_25 筛的第二部分从 f' 恢复到 f（注意这两步我们会在 f(i)[i\in\mathbb{P}\vee\operatorname{lpf}(i)>n^\alpha],f(i)[\operatorname{lpf}(i)>n^\alpha] 间切换，Min_25 筛时要保留 n^\alpha 以内质数）。取 \alpha=\frac16，此部分 O(\frac{n^\frac23}{\log n})。

对 g',h' 做杜教筛，考虑使用狄利克雷双曲线法。对 f\ast g=h，有 H(n)=\sum_{i=1}^{\lfloor\sqrt n\rfloor}f(i)G(\lfloor\frac ni\rfloor)+\sum_{i=1}^{\lfloor\sqrt n\rfloor}g(i)F(\lfloor\frac ni\rfloor)-F(\lfloor\sqrt n\rfloor)G(\lfloor\sqrt n\rfloor)，移项得到递归式 F(n)=H(n)-\sum_{i=1}^{\lfloor\sqrt n\rfloor}f(i)G(\lfloor\frac ni\rfloor)-\sum_{i=2}^{\lfloor\sqrt n\rfloor}g(i)F(\lfloor\frac ni\rfloor)+F(\lfloor\sqrt n\rfloor)G(\lfloor\sqrt n\rfloor)。因为 f',g' 的密度是 O(\log^{-1}n)，求块筛单点是 O(\frac{\sqrt n}{\log n}) 的。

线性筛 B 个 F，块筛卷积部分是 O(B+\frac{n}{\sqrt B\log n})=O((\frac n{\log n})^\frac23) 的。实际上我们有手段避免线性筛到 B。令 B=n^\frac23，考虑搜索 B 以内的 n^\alpha-rough 数的质因数分解。这样需要筛出 O(n^\frac 23) 的质数。我们发现要搜到的合数含有的质因数不超过 n^{\frac23-\frac16}=\sqrt n，考虑构造一个完全积性函数 f_p(p^k)=f(p)^k，将 f' 先设为 f_p 筛除 n^\alpha 内质数的结果，然后用上面的 DFS 修改贡献。这样就不需要搜质数了。

剩下的问题是 g,h 可能不是完全积性的。但是根据 PN 筛道理，我们可以在 O(\sqrt n\log n) 的时间内在 f(p)=g(p) 的两个积性函数间转换，因此将 f,g,h 切换为 f(p^k)\gets f(p)^k，用上述做法，不是瓶颈。复杂度 O(\frac{n^\frac23}{\log n})。

#include<bits/stdc++.h>
using namespace std;
int t,bn,bn1,bn2,num,prime[46345],cnt,pi[1664515],mu[1664515];
long long n,r[92686],g[92686],f[92686];
bool vis[46345];
int mu_init(int n,int prime[],int mu[],bool vis[],int cnt=0){
  mu[1]=1;
  for(int i=2;i<=n;i++)vis[i]=1;
  for(int i=2;i<=n;i++){
    if(vis[i])prime[++cnt]=i,mu[i]=-1;
    for(int j=1,v,p;i*prime[j]<=n;j++){
      p=prime[j],v=i*p,vis[v]=0;
      if(i%p==0){
        mu[v]=0;
        break;
      }
      else mu[v]=-mu[i];
    }
  }
  return cnt;
}
int id(long long x){
  return x<=bn?x:num-n/x+1;
}
void dfs(int dep,int now,int v){
  if(now>bn||!vis[now])f[id(now)]+=v+1;
  while(dep<=cnt&&prime[dep]<=min(bn1/now,bn)){
    int p=prime[dep],nxt=now;
    for(int i=1;1ll*nxt*p<=bn1;i++,nxt*=p)dfs(dep+1,nxt*p,i==1?-v:0);
    dep++;
  }
}
int main(){
  cnt=mu_init(46340,prime,mu,vis);
  for(int i=1;i<=46340;i++)pi[i]=pi[i-1]+vis[i];
  cin>>t;
  while(t--){
    cin>>n,bn=sqrt(n),bn1=pow(n,2.0/3),bn2=pow(n,1.0/6);
    for(long long l=1;l<=n;l=r[num]+1)r[++num]=n/(n/l),g[num]=r[num]-1;
    for(int j=1;j<=cnt&&prime[j]<=bn2;j++){
      int p=prime[j];
      long long t=p*p;
      for(int i=num;r[i]>=t;i--)g[i]-=g[id(r[i]/p)]-g[p-1];
    }
    for(int i=1;i<=num;i++)g[i]-=pi[min(r[i],(long long)bn2)]-1;
    for(int i=1;i<=num;i++)f[i]=0;
    dfs(pi[bn2]+1,1,1);
    for(int i=1;i<=num&&r[i]<=bn1;i++)f[i]+=f[i-1];
    for(int i=1;i<=num&&r[i]<=bn1;i++)f[i]-=g[i];
    vector<int>rough;
    for(int i=1;i<=bn;i++)if(g[i]!=g[i-1])rough.push_back(i);
    for(int i=1;i<=num;i++){
      if(r[i]<=bn1)continue;
      int bn=sqrt(r[i]);
      f[i]=f[bn]*g[bn]+1;
      for(int j=0;j<rough.size()&&rough[j]<=bn;j++)f[i]-=mu[rough[j]]*g[id(r[i]/rough[j])];
      for(int j=1;j<rough.size()&&rough[j]<=bn;j++)f[i]-=f[id(r[i]/rough[j])];
    }
    for(int i=1;i<=num;i++)f[i]-=pi[min(r[i],(long long)bn2)]+1;
    for(int j=pi[bn2];j>=1;j--){
      int p=prime[j];
      long long t=1ll*p*p;
      for(int i=num;r[i]>=t;i--)f[i]-=f[id(r[i]/p)]+j;
    }
    for(int i=1;i<=num;i++)f[i]++;
    long long phi=0;
    for(int i=1;i<=num;i++)phi+=(f[i]-f[i-1])*(n/r[i])*(n/r[i]+1)/2;
    cout<<phi<<' '<<f[num]<<'\n',num=0;
  }
  return 0;
}

关于 PN 筛求块筛：

定义 Powerful Number（PN）为所有质因数的次数都大于一的数。可以证明 PN 的个数是 O(\sqrt n) 级别的：显然 PN 都可以表示为形如 a^2b^3 的数，个数为 \sum_{i=1}^{\sqrt n}(\frac n{i^2})^\frac 13=O(\sqrt n)。

假如要求 f 的前缀和，构造一个容易求和的积性函数 g，满足对于任意质数 p，满足 f(p)=g(p)。

存在一个函数 h 使得 f=g\ast h，那么 f(p)=g(1)h(p)+g(p)h(1)=g(p)+h(p)=f(p)+h(p)。因此 h(p)=0。h 为积性函数，所以 h 只在 PN 处有值，因为其他函数都可以分离出一项次数为 1。

\begin{aligned}F(n)&=\sum_{i=1}^nf(i)\\&=\sum_{i=1}^n\sum_{d\mid i}g(\frac id)h(d)\\&=\sum_{d=1}^nh(d)\sum_{d\mid i}g(\frac i d)\\&=\sum_{\substack{d=1\\\text{d is PN}}}^nh(d)G(\lfloor\frac n d\rfloor)\end{aligned}

PN 可以通过搜索求出。现在只需要求 h(d) 的值。可以求所有 h(p^k) 的值，再乘出 PN 的值。若函数特殊，可以手玩。一般可以根据 f=g\ast h 有 f(p^k)=\sum_{i=0}^kg(p^i)h(p^{k-i})，则 h(p^k)=f(p^k)-\sum_{i=1}^kg(p^i)h(p^{k-i})，暴力求逆即可，只需要对根号以下质数这么做，复杂度不超过 O(\sqrt n)。

已知 g 的块筛，我们得以 O(\sqrt n) 求 f 的前缀和的单点。求块筛则可以做到 O(\sqrt n\log n)。

我们要对每个整除位置求 \sum_{ij\leq k}g(i)h(j)。对块筛 F 差分。若 ij\leq\sqrt n，枚举 ij 贡献到对应位置，O(\sqrt n\log n)。

剩下的整除位置可表示为 \lfloor\frac nk\rfloor(k\leq\sqrt n)。若 i>\sqrt n,j\leq\sqrt n，枚举 j,k，能贡献的 i 的范围是 \max(\lfloor\frac n{i(j+1)}\rfloor,\sqrt n)<i\leq \lfloor\frac n{ij}\rfloor，可用块筛表示。jk\leq\sqrt n，复杂度 O(\sqrt n\log n)。i\leq n,j>\sqrt n 同理。

若 i,j\leq\sqrt n，枚举 j 是 PN。若 i\leq\sqrt\frac nj，枚举 i 加到上面。否则 i>\sqrt\frac nj,k\leq\sqrt\frac nj。枚举 j,k，i 的范围是 \max(\lfloor\frac n{j(k+1)}\rfloor,\sqrt\frac nj)<i\leq\min(\lfloor\frac n{jk}\rfloor,\sqrt n)。此部分复杂度 \sum_{i=1}^{n^\frac14}\sqrt\frac n{i^2}=O(\sqrt n\log n)。

void powerful_number(long long n,int sg[],int sf[]){
  vector<pair<long long,int> >pn;
  for(int i=1;i<=cnt;i++){
    int p=prime[i],gp[40];
    long long pk=p;
    hp[i][0]=1;
    for(int j=1;pk<=n;j++,pk*=p){
      gp[j]=calcg(pk,p),hp[i][j]=calcf(p,j);
      for(int k=1;k<=j;k++)hp[i][j]=(hp[i][j]-1ll*gp[k]*hp[i][j-k]%mod+mod)%mod;
    }
  }
  dfs(1,1,1,n,pn),sort(pn.begin(),pn.end());
  for(int i=1,j=0,now=0;i<=num;i++){
    while(j<pn.size()&&pn[j].first<=r[i])now=(now+pn[j++].second)%mod;
    sh[i]=now;
  }
  for(int i=1;i<=bn;i++){
    for(int j=1;i*j<=bn;j++){
      int l=id(max(n/i/(j+1),(long long)bn)),r=id(n/i/j);
      if(sh[i]!=sh[i-1]){
        sf[i*j]=(sf[i*j]+1ll*(sh[i]-sh[i-1]+mod)*(sg[j]-sg[j-1]+mod))%mod;
        sf[num-j+1]=(sf[num-j+1]+1ll*(sh[i]-sh[i-1]+mod)*(sg[r]-sg[l]+mod))%mod;
      }
      sf[num-j+1]=(sf[num-j+1]+1ll*(sg[i]-sg[i-1]+mod)*(sh[r]-sh[l]+mod))%mod;
    }
  }
  for(int i=0;i<pn.size()&&pn[i].first<=bn;i++){
    long long x=pn[i].first,lim=sqrt(n/x);
    for(int j=bn/x+1;j<=lim;j++)sf[id(x*j)]=(sf[id(x*j)]+1ll*pn[i].second*(sg[j]-sg[j-1]+mod))%mod;
    for(int j=max((n+x*bn-1)/(x*bn)-1,1ll);j<=lim;j++){
      sf[num-j+1]=(sf[num-j+1]+1ll*pn[i].second*(sg[min(n/x/j,(long long)bn)]-sg[id(max(n/x/(j+1),(long long)lim))]+mod))%mod;
    }
  }
  for(int i=1;i<=num;i++)sf[i]=(sf[i]+sf[i-1])%mod;
}