筛法 学习笔记

埃拉托斯特尼筛法

一个质数数的倍数一定为合数，因此可以在枚举到一个质数时对其倍数打标记。如果从小到大枚举到一个数未被标记，说明这个数为质数，因为没有更小的质因数了。

int prime_init(int n,int prime[],bool a[],int cnt=0){
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i]){
      prime[++cnt]=i;
      for(int j=i*2;j<=n;j+=i)a[j]=0;
    }
  }
  return cnt;
}

看似暴力，实际上每个数只会被质因数筛，而 \(\omega\) 的增长是很慢的。复杂度 \(O(n\log\log n)\)。

线性筛

线性筛质数

假设钦定每个数都只被最小的质因数筛。对于每个数 \(i\)，枚举可能的最小质因数 \(j\)，对 \(ij\) 打标记。当 \(j\mid i\)，\(j\) 继续变大会让 \(ij\) 的最小质因数变为 \(i\) 的最小质因数，因此退出。复杂度 \(O(n)\)。

int prime_init(int n,int prime[],bool a[],int cnt=0){
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0;
      if(i%prime[j]==0)break;
    }
  }
  return cnt;
}

线性筛积性函数

可以借助每个数只会被最小质因数筛的性质 \(O(n)\) 计算大部分积性函数 \(1\sim n\) 的值。

以欧拉函数 \(\varphi\) 为例，设筛的数为 \(x\)，\(x\) 被 \(ij\) 筛，\(j\) 为 \(x\) 的最小质因数，分三类讨论：

  1. \(x\) 为质数，直接根据定义得 \(\varphi(x)=x-1\)。

  2. \(i\nmid j\)：\(\gcd(i,j)=1\)，根据积性函数性质可得 \(\varphi(x)=\varphi(i)\varphi(j)\)。

  3. \(i\mid j\)：\(\varphi(ij)=ij\prod p_i^{k_i}=j\varphi(i)\)。

int phi_init(int n,int prime[],int phi[],bool a[],int cnt=0){
  phi[1]=1;
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i,phi[i]=i-1;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0;
      if(i%prime[j]==0){
        phi[i*prime[j]]=phi[i]*prime[j];
        break;
      }
      else phi[i*prime[j]]=phi[i]*(prime[j]-1);
    }
  }
  return cnt;
}

另一个常用的积性函数是莫比乌斯函数 \(\mu\)。

  1. \(x\) 为质数，根据定义得 \(\mu(x)=-1\)。

  2. \(i\nmid j\)：根据积性函数性质可得 \(\mu(x)=\mu(i)\mu(j)=-\mu(i)\)。

  3. \(i\mid j\)：\(i\) 中 \(j\) 的次数大于一，\(\mu(x)=0\)。

int mu_init(int n,int prime[],int mu[],bool a[],int cnt=0){
  mu[1]=1;
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i,mu[i]=-1;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0;
      if(i%prime[j]==0){
        mu[i*prime[j]]=0;
        break;
      }
      else mu[i*prime[j]]=-mu[i];
    }
  }
  return cnt;
}

然而当 \(i\mid j\) 时不是所有函数都能求，有一个通用的方法。额外记录数组 \(low\) 表示最小质因数的最高幂次。当 \(i\mid j\) 时，把 \(ij\) 中所有 \(j\) 单独拿出来，有 \(f(ij)=f(low_{ij})f(\frac i{low_i})\)。

有特殊情况：当 \(x\) 为质数的幂，只能根据定义求解或者从 \(i\) 递推。

#include<bits/stdc++.h>
using namespace std;
int mu_init(int n,int prime[],int low[],int f[],bool a[],int cnt=0){
  low[1]=f[1]=1;
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i,low[i]=i/*f[i]=f(i)*/;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0;
      if(i%prime[j]==0){
        low[i*prime[j]]=low[i]*prime[j];
        if(i==low[i])/*f[i*prime[j]]=f(p^k)*/
        else f[i*prime[j]]=f[i/low[i]]*f[low[i]*prime[j]];
        break;
      }
      else low[i*prime[j]]=prime[j],f[i*prime[j]]=f[i]*f[prime[j]];
    }
  }
  return cnt;
}

杜教筛

杜教筛可以在低于线性的时间内计算一类积性函数的前缀和。

令 \(F(n)=\sum_{i=1}^n f(x)\)。构造两个积性函数 \(g,h\) 使 \(h=f\ast g\)，且 \(g,h\) 可以快速求前缀和。

\[\begin{aligned}\sum_{i=1}^nh(i)&=\sum_{i=1}^n\sum_{d\mid i}f(\frac i d)g(d)\\&=\sum_{d=1}^ng(d)\sum_{d\mid i}f(\frac i d)\\&=\sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor\frac n d\rfloor}f(i)\\&=g(1)F(n)+\sum_{d=2}^ng(d)F(\lfloor\frac n d\rfloor)\\F(n)&=\sum_{i=1}^nh(i)-\sum_{d=2}^ng(d)F(\lfloor\frac n d\rfloor)\end{aligned} \]

因此可以在整除分块里递归。复杂度 \(O(n^{\frac 3 4})\)。

复杂度证明：

\[\begin{aligned}\sum_{i=1}^n\frac{2}{\sqrt i+\sqrt{i+1}}&<\sum_{i=1}^n\frac{1}{\sqrt i}<\sum_{i=1}^n\frac{2}{\sqrt{i-1}+\sqrt i}(\text{放缩})\\2\sum_{i=1}^n\sqrt{i+1}-\sqrt i&<\sum_{i=1}^n\frac{1}{\sqrt i}<2\sum_{i=1}^n\sqrt i-\sqrt{i-1}\\2\sqrt{n+1}-2&<\sum_{i=1}^n\frac{1}{\sqrt i}<2\sqrt n \end{aligned}\]

枚举每个 \(\lfloor\frac n x\rfloor\)：

\[T(n)=O\left(\sum_{i=1}^{\lfloor\sqrt n\rfloor}\sqrt i+\sum_{i=1}^{\lfloor\sqrt n\rfloor}\sqrt{\frac n i}\right)=O\left(\sum_{i=1}^{\lfloor\sqrt n\rfloor}\sqrt{\frac n i}\right)=O\left(\sqrt n\sum_{i=1}^{\lfloor\sqrt n\rfloor}\sqrt{\frac 1 i}\right)=O\left(\sqrt n\sqrt{\sqrt n}\right)=O(n^{\frac 3 4}) \]

一个优化：可以用线性筛预处理一部分的值。前 \(n^{\frac 2 3}\) 的值，可以证明复杂度 \(O(n^{\frac 2 3})\)。

\[T(n)=O\left(B+\sum_{i=1}^{\lfloor\frac n B\rfloor}\sqrt{\frac n i}\right)=O\left(B+\frac{n}{\sqrt B}\right) \]

\(B=n^{\frac 2 3}\) 最优，复杂度 \(O(n^{\frac 2 3})\)。

另外，还需要记忆化以保证复杂度，有两种写法：一种是用 map 或 unordered_map，另一种是观察到递归到的数 \(x\) 都为 \(\lfloor\frac n x\rfloor\) 的形式，可以将大于 \(\sqrt n\) 的 \(x\) 映射到下标 \(\sqrt n+\lfloor\frac n x\rfloor\)。第二种写法有一定局限性，比如 \(n\) 变化时要清空，外面套了高维整除分块时无法使用。

发现这样也顺便求出了所有 \(\lfloor\frac n x\rfloor\) 处的前缀和，我们称所有 \(\lfloor\frac n x\rfloor\) 处的前缀和为块筛。

举个例子：\(\varphi\) 可以构造 \(\operatorname{Id}=\varphi\ast\mathbf{1}\)，\(\mu\) 可以构造 \(\varepsilon=\mu\ast\mathbf{1}\)。

long long getsum1(int x){
  if(x<=bn)return phi[x];
  if(sum1.count(x))return sum1[x];
  long long ans=1ll*x*(x+1)/2;
  for(int l=2,r;l<=x;l=r+1)r=x/(x/l),ans-=(r-l+1)*getsum1(x/l);
  return sum1[x]=ans;
}
long long getsum2(int x){
  if(x<=bn)return mu[x];
  if(sum2.count(x))return sum2[x];
  long long ans=1;
  for(int l=2,r;l<=x;l=r+1)r=x/(x/l),ans-=(r-l+1)*getsum2(x/l);
  return sum2[x]=ans;
}

P3768

\[\begin{aligned}\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)&=\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^nij[\gcd(i,j)=d]\\&=\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac n d\rfloor}ij[\gcd(i,j)=1]\\&=\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac n d\rfloor}ij\sum_{p\mid i,p\mid j}\mu(p)\\&=\sum_{d=1}^nd^3\sum_{p=1}^{\lfloor\frac n d\rfloor}\sum_{p\mid i}\sum_{p\mid j}ij\\&=\sum_{d=1}^nd^3\sum_{p=1}^{\lfloor\frac n d\rfloor}p^2\mu(p)\sum_{i=1}^{\lfloor\frac n{dp}\rfloor}\sum_{j=1}^{\lfloor\frac n{dp}\rfloor}ij\\&=\sum_{d=1}^nd^3\sum_{p=1}^{\lfloor\frac n d\rfloor}p^2\mu(p)\frac{\lfloor\frac n{dp}\rfloor(\lfloor\frac n{dp}\rfloor+1)}{4}\\&=\sum_{T=1}^n\frac{\lfloor\frac n T\rfloor(\lfloor\frac n T\rfloor+1)}{4}\sum_{p\mid T}p^2\mu(p)\left(\frac T p\right)^3\\&=\sum_{T=1}^n\frac{\lfloor\frac n T\rfloor(\lfloor\frac n T\rfloor+1)}{4}\sum_{p\mid T}\mu(p)\frac T p\\&=\sum_{T=1}^n\frac{\lfloor\frac n T\rfloor(\lfloor\frac n T\rfloor+1)}{4}T^2\varphi(T)\end{aligned} \]

问题在于求 \(T^2\varphi(p)\) 的前缀和。考虑杜教筛，由于多了一个 \(\operatorname{Id}_2\)，\(g\) 也要配 \(\operatorname{Id}_2\)，卷出来是 \(\operatorname{Id}_3\)。

整除分块套杜教筛虽然调用了 \(\sqrt n\) 次，但外层的整除分块和杜教筛的整除分块是一样的，所以复杂度还是 \(O(n^{\frac 2 3})\)。

#include<bits/stdc++.h>
using namespace std;
int p,prime[4641593],phi[4641593],sum[2159],inv4,inv6,ans;
bool a[4641593],vis[2159];
long long n,bn;
template<typename T>T qpow(T a,T b,T n,T ans=1){
  for(a%=n;b;b>>=1)b&1&&(ans=1ll*ans*a%n),a=1ll*a*a%n;
  return ans%n;
}
template<typename T>T inv(T a,T b)
{
  return qpow(a,b-2,b);
}
int phi_init(int n,int prime[],int phi[],bool a[],int cnt=0){
  phi[1]=1;
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i,phi[i]=i-1;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0;
      if(i%prime[j]==0){
        phi[i*prime[j]]=phi[i]*prime[j];
        break;
      }
      else phi[i*prime[j]]=phi[i]*(prime[j]-1);
    }
  }
  return cnt;
}
int getsum(long long x,int ans=0){
  if(x<=bn)return phi[x];
  if(vis[n/x])return sum[n/x];
  ans=x%p*x%p*(x+1)%p*(x+1)%p*inv4%p;
  for(long long l=2,r;l<=x;l=r+1)r=x/(x/l),ans=(ans-(r%p*(r+1)%p*(2*r%p+1)%p*inv6%p-(l-1)%p*l%p*(2*l%p-1)%p*inv6%p+p)%p*getsum(x/l)%p+p)%p;
  return vis[n/x]=1,sum[n/x]=ans;
}
int main(){
  cin>>p>>n,inv4=inv(4ll,(long long)p),inv6=inv(6ll,(long long)p),bn=pow(n,2.0/3),phi_init(bn,prime,phi,a);
  for(int i=1;i<=bn;i++)phi[i]=(phi[i-1]+1ll*phi[i]*i%p*i%p)%p;
  for(long long l=1,r;l<=n;l=r+1)r=n/(n/l),ans=(ans+(n/l)%p*(n/l)%p*(n/l+1)%p*(n/l+1)%p*inv4%p*(getsum(r)-getsum(l-1)+p)%p)%p;
  return cout<<ans<<'\n',0;
}

Powerful Number 筛

定义 Powerful Number（PN）为所有质因数的次数都大于一的数。可以证明 PN 的个数是 \(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(d)h(\frac i d)\\&=\sum_{d=1}^nh(d)\sum_{d\mid i}^ng(\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})\)，递推即可。

P5325

构造 \(g(n)=n\varphi(n)\)。这题中需要用杜教筛求 \(g\) 的块筛，复杂度 \(O(n^{\frac 2 3})\)。

#include<bits/stdc++.h>
using namespace std;
const int mod=1000000007;
int bn,prime[4641593],f[4641593],h[325166][40],cnt,ans;
long long n,p[325166][40];
bool temp[4641593];
unordered_map<long long,int>sum;
int phi_init(int n,int prime[],int phi[],bool a[],int cnt=0){
  phi[1]=1;
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i,phi[i]=i-1;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0;
      if(i%prime[j]==0){
        phi[i*prime[j]]=1ll*phi[i]*prime[j]%mod;
        break;
      }
      else phi[i*prime[j]]=1ll*phi[i]*(prime[j]-1+mod)%mod;
    }
  }
  return cnt;
}
int getsum(long long x){
  if(x<=bn)return f[x];
  if(sum.count(x))return sum[x];
  int ans=1ll*x%mod*((x+1)%mod)%mod*((2*x+1)%mod)%mod*166666668%mod;
  for(long long l=2,r;l<=x;l=r+1)r=x/(x/l),ans=(ans-1ll*(l+r)%mod*((r-l+1)%mod)%mod*500000004%mod*getsum(x/l)%mod+mod)%mod;
  return sum[x]=ans;
}
void dfs(int d,long long now,int nh){
  if(d>cnt||now>n/p[d][2]){
    ans=(ans+1ll*nh*getsum(n/now)%mod)%mod;
    return;
  }
  dfs(d+1,now,nh);
  for(int i=2;now<=n/p[d][i];i++)dfs(d+1,now*p[d][i],1ll*nh*h[d][i]%mod);
}
int calcf(long long x){
  return x%=mod,x*(x-1+mod)%mod;
}
int calcg(long long x,int p){
  return x%mod*((x/p*(p-1))%mod)%mod;
}
int main(){
  cin>>n,bn=pow(n,2.0/3),cnt=phi_init(bn,prime,f,temp);
  for(int i=1;i<=bn;i++)f[i]=(1ll*f[i]*i%mod+f[i-1])%mod;
  for(int i=1;i<=cnt;i++){
    p[i][0]=h[i][0]=1,p[i][1]=prime[i];
    for(int j=1;p[i][j]<=n;j++,p[i][j]=p[i][j-1]*prime[i]){
      h[i][j]=calcf(p[i][j]);
      for(int k=1;k<=j;k++)h[i][j]=(h[i][j]-1ll*calcg(p[i][k],prime[i])*h[i][j-k]%mod+mod)%mod;
    }
  }
  return dfs(1,1,1),cout<<ans<<'\n',0;
}

Min_25 筛

概述

记号规定：

\(p_k\)：第 \(k\) 小的质数，\(p_0=1\)。

\(\operatorname{lpf}(n)\)：\(n\) 的最小质因数，\(\operatorname{lpf}(1)=1\)。

Min_25 筛可以求一类积性函数的前缀和，若一个质数为 \(p\)，要求 \(f(p)\) 可以被拆成若干能快速求前缀和的完全积性函数（一般是幂函数，也就是 \(f(p)\) 为低阶多项式），且 \(f(p^k)\) 可以快速求值。

质数前缀统计

先求出函数 \(f\) 在质数处的和，即 \(g(n)=\sum_{i=1}^nf(i)[i\in\mathbb{P}]\)。可以将 \(f\) 拆成若干项，每项都是幂函数，最后加减起来。

我们借鉴埃氏筛的思想，设 \(g(n,j)=\sum_{i=1}^nf(i)[i\in\mathbb{P}\vee\operatorname{lpf}(i)>p_j]\)，也就是所有质数与最小质因数大于 \(p_j\) 的数的取值的和。假设合数处也遵循质数的规律，可以看做经过 \(j\) 轮筛法，剩下的数 \(f\) 的和。那么用所有不大于 \(\sqrt n\) 的质数筛完后的 \(g(n)\) 就是答案。

当进行了新一轮的筛法，所有最小质因数为 \(p_j\) 的合数会被筛去。由于 \(f\) 此时是完全积性函数，可以提一个 \(f(p_j)\)，剩下的应为最小质因数大于 \(p_{j-1}\) 的合数，也就是 \(g\left(\lfloor\frac n{p_j}\rfloor,j-1\right)-g(p_{j-1},j-1)\)。因此有递推式 \(g(n,j)=g(n,j-1)-f(p_j)\left(g\left(\lfloor\frac n{p_j}\rfloor,j-1\right)-g(p_{j-1},j-1)\right)\)。这里只需计算 \(\lfloor\frac n x\rfloor\) 处的点值（后面一项的 \(p_{j-1}\leq\sqrt n\)，所以也能表示），可以一次整除分块把这些位置找出。注意这样的位置最多有 \(2\sqrt n\) 个。

边界条件：\(g(n,0)=\sum_{i=2}^n f(i)\)。因为 \(f(p)\) 是幂函数，有通项公式。

复杂度证明：

\[T(n)=O\left(\sum_{i=1}^{\lfloor\sqrt n\rfloor}\pi\left(\sqrt{\frac n i}\right)\right)=O\left(\sum_{i=1}^{\lfloor\sqrt n\rfloor}\frac{\sqrt{\frac n i}}{\log\sqrt{\frac n i}}\right)=O\left(\frac{1}{\log n}\sum_{i=1}^{\lfloor\sqrt n\rfloor}\sqrt{\frac n i}\right)=O(\frac{n^{\frac 3 4}}{\log n}) \]

空间可以使用类似上面杜教筛的方法存 \(g\)，然后压掉 \(j\) 一维，\(O(\sqrt n)\)。最后的 \(g\) 就是 \(\sum_{i=1}^nf(i)[i\in\mathbb{P}]\) 的块筛。

P5493

模板题。求 \(k\) 次幂的前缀和可以拉插，求出块筛后暴力计算这个式子。

#include<bits/stdc++.h>
using namespace std;
int k,mod,bn,prime[200005],f[200005],cnt,y[15],fac[15],vac[15],pre[15],suf[15],num,g[400005],ans;
long long n,r[400005],t;
bool temp[200005];
template<typename T>T qpow(T a,T b,T n,T ans=1){
  for(a%=n;b;b>>=1)b&1&&(ans=1ll*ans*a%n),a=1ll*a*a%n;
  return ans%n;
}
template<typename T>T inv(T a,T b){
  return qpow(a,b-2,b);
}
int prime_init(int n,int prime[],bool a[],int cnt=0){
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0;
      if(i%prime[j]==0)break;
    }
  }
  return cnt;
}
int calcg(long long n,int k,int ans=0){
  n%=mod,pre[0]=suf[k+3]=1;
  for(int i=1;i<=k+2;i++)pre[i]=1ll*pre[i-1]*(n-i+mod)%mod;
  for(int i=k+2;i>=1;i--)suf[i]=1ll*suf[i+1]*(n-i+mod)%mod;
  for(int i=1;i<=k+2;i++)ans=(ans+(((k+2-i)&1?-1ll:1ll)*y[i]*pre[i-1]%mod*suf[i+1]%mod*vac[i-1]%mod*vac[k+2-i]%mod+mod)%mod)%mod;
  return ans;
}
int id(long long x){
  return x<=bn?x:num-n/x+1;
}
int main(){
  cin>>n>>k>>mod,bn=sqrt(n),cnt=prime_init(bn,prime,temp);
  for(int i=1;i<=cnt;i++){
    f[i]=1;
    for(int j=1;j<=k;j++)f[i]=1ll*f[i]*prime[i]%mod;
  }
  for(int i=1;i<=k+2;i++){
    y[i]=1;
    for(int j=1;j<=k;j++)y[i]=1ll*y[i]*i%mod;
    y[i]=(y[i]+y[i-1])%mod;
  }
  fac[0]=1;
  for(int i=1;i<=k+2;i++)fac[i]=1ll*fac[i-1]*i%mod;
  vac[0]=1,vac[k+2]=inv(fac[k+2],mod);
  for(int i=k+1;i>=1;i--)vac[i]=1ll*vac[i+1]*(i+1)%mod;
  for(long long l=1;l<=n;l=r[num]+1)r[++num]=n/(n/l),g[num]=(calcg(r[num],k)-1+mod)%mod;
  for(int j=1;j<=cnt;j++){
    t=1ll*prime[j]*prime[j];
    for(int i=num;r[i]>=t;i--){
      g[i]=(g[i]-1ll*f[j]*(g[id(r[i]/prime[j])]-g[id(prime[j-1])]+mod)%mod+mod)%mod;
    }
  }
  for(int i=1;1ll*i*i<=n;i++)ans=(ans+1ll*i*i%mod*g[id(n/i)]%mod)%mod;
  return cout<<ans<<'\n',0;
}

可以做到 \(O((\frac{n}{\log n})^{\frac 2 3})\)。

考虑另一种维护 \(g\) 的方法：用一个树状数组，进行一轮筛法时在树状数组上修改。

当 \(p_j\geq\sqrt n\) 时，\(g(n,j)\) 不再变化。这时把 \(g\) 存起来，不再在树状数组上查询。

可以根号分治，设阈值为 \(B_1\)，\(B_1\) 以上的值用 DP，\(B_1\) 以下的用树状数组。

还可以再分，对于小于等于 \(B_2\) 的质数使用不经树状数组优化的 DP。

当 \(B_1=((\frac{n}{\log n})^{\frac 2 3}),B_2=n^{\frac 1 6}\) 时复杂度为 \(O((\frac{n}{\log n})^{\frac 2 3})\)。

#include<bits/stdc++.h>
using namespace std;
int k,mod,bn,b1,b2,prime[200005],f[2423524],cnt,y[15],fac[15],vac[15],pre[15],suf[15],num,g[400005],tr[2423524],ans;
long long n,r[400005],t;
bool vis[2423524];
template<typename T>T qpow(T a,T b,T n,T ans=1){
  for(a%=n;b;b>>=1)b&1&&(ans=1ll*ans*a%n),a=1ll*a*a%n;
  return ans%n;
}
template<typename T>T inv(T a,T b){
  return qpow(a,b-2,b);
}
int prime_init(int n,int prime[],bool a[],int cnt=0){
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0;
      if(i%prime[j]==0)break;
    }
  }
  return cnt;
}
int calcg(long long n,int k,int ans=0){
  n%=mod,pre[0]=suf[k+3]=1;
  for(int i=1;i<=k+2;i++)pre[i]=1ll*pre[i-1]*(n-i+mod)%mod;
  for(int i=k+2;i>=1;i--)suf[i]=1ll*suf[i+1]*(n-i+mod)%mod;
  for(int i=1;i<=k+2;i++)ans=(ans+(((k+2-i)&1?-1ll:1ll)*y[i]*pre[i-1]%mod*suf[i+1]%mod*vac[i-1]%mod*vac[k+2-i]%mod+mod)%mod)%mod;
  return ans;
}
void add(int x,int k){
  for(;x<=b1;x+=(x&-x))tr[x]=(tr[x]+k+mod)%mod;
}
int query(int x,int ans=0){
  for(;x;x-=(x&-x))ans=(ans+tr[x])%mod;
  return ans;
}
int id(long long x){
  return x<=bn?x:num-n/x+1;
}
int main(){
  cin>>n>>k>>mod,bn=sqrt(n),b1=max((int)pow(n/log(n),2.0/3),bn),b2=pow(n,1.0/6),cnt=prime_init(bn,prime,vis);
  for(int i=1;i<=b1;i++){
    f[i]=1;
    for(int j=1;j<=k;j++)f[i]=1ll*f[i]*i%mod;
  }
  for(int i=1;i<=k+2;i++){
    y[i]=1;
    for(int j=1;j<=k;j++)y[i]=1ll*y[i]*i%mod;
    y[i]=(y[i]+y[i-1])%mod;
  }
  fac[0]=1;
  for(int i=1;i<=k+2;i++)fac[i]=1ll*fac[i-1]*i%mod;
  vac[0]=1,vac[k+2]=inv(fac[k+2],mod);
  for(int i=k+1;i>=1;i--)vac[i]=1ll*vac[i+1]*(i+1)%mod;
  for(long long l=1;l<=n;l=r[num]+1)r[++num]=n/(n/l),g[num]=(calcg(r[num],k)-1+mod)%mod;
  memset(vis,0,sizeof(vis));
  for(int i=2;i<=b1;i++)tr[i]=f[i];
  for(int j=1;j<=cnt&&prime[j]<=b2;j++){
    t=1ll*prime[j]*prime[j];
    for(int i=num;r[i]>=t;i--){
      g[i]=(g[i]-1ll*f[prime[j]]*(g[id(r[i]/prime[j])]-g[id(prime[j-1])]+mod)%mod+mod)%mod;
    }
    if(t<=b1)for(int i=t;i<=b1;i+=prime[j])if(!vis[i])vis[i]=1,tr[i]=0;
  }
  for(int i=1;i<=b1;i++)if(i+(i&-i)<=b1)tr[i+(i&-i)]=(tr[i+(i&-i)]+tr[i])%mod;
  for(int j=1,p=1;j<=cnt;j++){
    if(prime[j]<=b2)continue;
    t=1ll*prime[j]*prime[j];
    while(r[p]<b1&&r[p]<t&&p<=num)g[p]=query(r[p]),p++;
    for(int i=num;r[i]>=b1&&r[i]>=t;i--){
      if(r[i]/prime[j]>b1||r[i]/prime[j]<t)g[i]=(g[i]-1ll*f[prime[j]]*(g[id(r[i]/prime[j])]-g[id(prime[j-1])]+mod)%mod+mod)%mod;
      else g[i]=(g[i]-1ll*f[prime[j]]*(query(r[i]/prime[j])-g[id(prime[j-1])]+mod)%mod+mod)%mod;
    }
    if(t<=b1)for(int i=t;i<=b1;i+=prime[j])if(!vis[i])vis[i]=1,add(i,-f[i]);
  }
  for(int i=1;1ll*i*i<=n;i++)ans=(ans+1ll*i*i%mod*g[id(n/i)]%mod)%mod;
  return cout<<ans<<'\n',0;
}

第二部分

法一：

类似地，设 \(S(n,j)=\sum_{i=1}^n f(i)[\operatorname{lpf}(i)>p_j]\)。

考虑如何计算。首先我们知道质数的贡献是 \(g(n)-g(p_j)\)。然后枚举最小质因数的次数并提出来，有：

\[S(n,j)=g(n)-g(p_j)+\sum_{j<k,p_k<\sqrt n}\sum_{1\leq e,p_k^e\leq n}f(p_k^e)\left(S\left(\left\lfloor\frac{n}{p_k^e}\right\rfloor,k\right)+[e\ne 1]\right) \]

\([e\ne 1]\) 是因为 \(S\) 不包含 \(1\)。简化一下，当 \(\left\lfloor\frac{n}{p_k^e}\right\rfloor\leq p_k\)，\(S\left(\left\lfloor\frac{n}{p_k^e}\right\rfloor,k\right)=0\)。

\[S(n,j)=g(n)-g(p_j)+\sum_{j<k,p_k<\sqrt n}\sum_{1\leq e,p_k^{e+1}\leq n}f(p_k^e)S\left(\left\lfloor\frac{n}{p_k^e}\right\rfloor,k\right)+f(p_k^{e+1}) \]

直接递归计算 \(S(n,0)+1\) 即可。边界条件为当 \(p_j\geq n\) 时返回 \(0\)。复杂度为 \(O(n^{1-\varepsilon})\)，但在实际运行中表现良好，与 \(O(\frac{n^{\frac 3 4}}{\log n})\) 相近。缺点是无法求块筛。

法二：

再设 \(S(n,j)=\sum_{i=1}^nf(i)[i\in\mathbb{P}\vee\operatorname{lpf}(i)>p_j]\)，求出正确的合数的贡献。考虑从 \(S(n,j+1)\) 转移到 \(S(n,j)\)。枚举 \(p_{j+1}\) 的次数，特判质数和 \(1\)，有：

\[S(n,j)=S(n,j+1)+\sum_{1\leq e,p_{j+1}^e\leq n}f(p_{j+1}^e)\left(S\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,j+1\right)-g\left(\min\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,p_{j+1}\right)\right)+[e\ne 1]\right) \]

简化一下，当 \(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor\leq p_{j+1}\)，\(S\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,j+1\right)-g\left(\min\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,p_{j+1}\right)\right)\) 实际上已经不产生贡献。变为：

\[S(n,j)=S(n,j+1)+\sum_{1\leq e,p_{j+1}^{e+1}\leq n}f(p_{j+1}^e)\left(S\left(\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor,j+1\right)-g\left(p_{j+1}\right)\right)+f(p_{j+1}^{e+1}) \]

边界条件：\(S(n,+\infty)=g(n)\)。质数的幂的密度和质数相同，复杂度同第一步，运行速度比法一慢。但是求出了块筛。

P5325

将 \(f(p)\) 拆成 \(p^2-p\)。

#include<bits/stdc++.h>
using namespace std;
const int mod=1000000007;
int bn,prime[100005],cnt,num,g[200005],g1[200005],ans;
long long n,r[200005],t;
bool temp[100005];
int prime_init(int n,int prime[],bool a[],int cnt=0){
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0;
      if(i%prime[j]==0)break;
    }
  }
  return cnt;
}
int calcf(long long x){
  return x%=mod,x*(x-1+mod)%mod;
}
int calcg(long long x){
  return x%=mod,x*(x+1)%mod*(2*x+1)%mod*166666668%mod;
}
int calcg1(long long x){
  return x%=mod,x*(x+1)%mod*500000004%mod;
}
int id(long long x){
  return x<=bn?x:num-n/x+1;
}
int S(long long n,int j){
  if(prime[j]>=n)return 0;
  int ans=(g[id(n)]-g[id(prime[j])]+mod)%mod;
  for(int k=j+1;k<=cnt&&1ll*prime[k]*prime[k]<=n;k++){
    long long p=prime[k];
    for(int e=1;p*prime[k]<=n;e++,p*=prime[k]){
      ans=(ans+(1ll*calcf(p)*S(n/p,k)%mod+calcf(p*prime[k]))%mod)%mod;
    }
  }
  return ans;
}
int main(){
  cin>>n,bn=sqrt(n),cnt=prime_init(bn,prime,temp);
  for(long long l=1;l<=n;l=r[num]+1)r[++num]=n/(n/l),g[num]=(calcg(r[num])-1+mod)%mod,g1[num]=(calcg1(r[num])-1+mod)%mod;
  for(int j=1;j<=cnt;j++){
    t=1ll*prime[j]*prime[j];
    for(int i=num;r[i]>=t;i--){
      g[i]=(g[i]-1ll*prime[j]*prime[j]%mod*(g[id(r[i]/prime[j])]-g[id(prime[j-1])]+mod)%mod+mod)%mod;
      g1[i]=(g1[i]-1ll*prime[j]*(g1[id(r[i]/prime[j])]-g1[id(prime[j-1])]+mod)%mod+mod)%mod;
    }
  }
  for(int i=1;i<=num;i++)g[i]=(g[i]-g1[i]+mod)%mod;
  return cout<<(S(n,0)+1)%mod<<'\n',0;
}

#include<bits/stdc++.h>
using namespace std;
const int mod=1000000007;
int bn,prime[100005],cnt,num,g[200005],g1[200005],s[200005],ans;
long long n,r[200005],t;
bool temp[100005];
int prime_init(int n,int prime[],bool a[],int cnt=0){
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0;
      if(i%prime[j]==0)break;
    }
  }
  return cnt;
}
int calcf(long long x){
  return x%=mod,x*(x-1+mod)%mod;
}
int calcg(long long x){
  return x%=mod,x*(x+1)%mod*(2*x+1)%mod*166666668%mod;
}
int calcg1(long long x){
  return x%=mod,x*(x+1)%mod*500000004%mod;
}
int id(long long x){
  return x<=bn?x:num-n/x+1;
}
int main(){
  cin>>n,bn=sqrt(n),cnt=prime_init(bn,prime,temp);
  for(long long l=1;l<=n;l=r[num]+1)r[++num]=n/(n/l),g[num]=(calcg(r[num])-1+mod)%mod,g1[num]=(calcg1(r[num])-1+mod)%mod;
  for(int j=1;j<=cnt;j++){
    t=1ll*prime[j]*prime[j];
    for(int i=num;r[i]>=t;i--){
      g[i]=(g[i]-1ll*prime[j]*prime[j]%mod*(g[id(r[i]/prime[j])]-g[id(prime[j-1])]+mod)%mod+mod)%mod;
      g1[i]=(g1[i]-1ll*prime[j]*(g1[id(r[i]/prime[j])]-g1[id(prime[j-1])]+mod)%mod+mod)%mod;
    }
  }
  for(int i=1;i<=num;i++)s[i]=g[i]=(g[i]-g1[i]+mod)%mod;
  for(int j=cnt;j>=1;j--){
    t=1ll*prime[j]*prime[j];
    for(int i=num;r[i]>=t;i--){
      long long p=prime[j];
      for(int e=1;p*prime[j]<=r[i];e++,p*=prime[j]){
        s[i]=(s[i]+(1ll*calcf(p)*(s[id(r[i]/p)]-g[id(prime[j])]+mod)%mod+calcf(p*prime[j]))%mod)%mod;
      }
    }
  }
  return cout<<(s[num]+1)%mod<<'\n',0;
}

LOJ6235

求质数处的取值只需要用到 Min_25 筛的第一步。这题相当于对 \(\mathbf{1}\) 求质数处函数值的和。

这题有一个简洁的写法：将线性筛与 Min_25 筛合并，当 \(g(n,j)\ne g(n-1,j)\) 时 \(n\) 为质数。

#include<bits/stdc++.h>
using namespace std;
int bn,num,cnt;
long long n,r[632459],g[632459],t;
int id(long long x){
  return x<=bn?x:num-n/x+1;
}
int main(){
  cin>>n,bn=sqrt(n);
  for(long long l=1;l<=n;l=r[num]+1)r[++num]=n/(n/l),g[num]=r[num]-1;
  for(int i=2;i<=bn;i++){
    if(g[i]==g[i-1])continue;
    t=1ll*i*i;
    for(int j=num;r[j]>=t;j--)g[j]-=g[id(r[j]/i)]-cnt;
    cnt++;
  }
  return cout<<g[num]<<'\n',0;
}

LOJ6028

有时我们可以修改 Min_25 筛的 DP 进行一些奇怪的计数。

增加一维余数，则转移方程变为 \(g(n,j,p_jr\bmod m)=g(n,j-1,p_jr\bmod m)-f(p_j)\left(g\left(\lfloor\frac n{p_j}\rfloor,j-1,r\right)-g(p_{j-1},j-1,r)\right)\)。

#include<bits/stdc++.h>
using namespace std;
int m,bn,num;
long long n,r[346415],g[346415][12],cnt[12],t;
int id(long long x){
  return x<=bn?x:num-n/x+1;
}
int main(){
  cin>>n>>m,bn=sqrt(n);
  for(long long l=1;l<=n;l=r[num]+1){
    r[++num]=n/(n/l);
    for(int i=0;i<m;i++)g[num][i]=r[num]/m+(i&&r[num]%m>=i);
    g[num][1]--;
  }
  for(int i=2;i<=bn;i++){
    if(g[i][i%m]==g[i-1][i%m])continue;
    t=1ll*i*i;
    for(int j=num;r[j]>=t;j--)for(int k=0;k<m;k++)g[j][i*k%m]-=g[id(r[j]/i)][k]-cnt[k];
    cnt[i%m]++;
  }
  for(int i=0;i<m;i++)cout<<g[id(n)][i]<<'\n';
  return 0;
}

LOJ6783

观察质数处的取值，发现 \(f(p)=\begin{cases}p+1,\,&p=2\\p-1,\,&p\ne 2\end{cases}\)。第一步可以把 \(f(p)\) 当做 \(p-1\)，最后把 \(1\) 以外的都加 \(2\)。这题需要求出块筛。

#include<bits/stdc++.h>
using namespace std;
const int mod=1000000007;
int bn,prime[316232],cnt,num,g[632454],g1[632454],s[632454],ans;
long long n,r[632454],t;
bool temp[316232];
int prime_init(int n,int prime[],bool a[],int cnt=0){
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0;
      if(i%prime[j]==0)break;
    }
  }
  return cnt;
}
int calcg(long long x){
  return x%=mod,x*(x+1)%mod*500000004%mod;
}
int calcg1(long long x){
  return x%mod;
}
int id(long long x){
  return x<=bn?x:num-n/x+1;
}
int main(){
  cin>>n,bn=sqrt(n),cnt=prime_init(bn,prime,temp);
  for(long long l=1;l<=n;l=r[num]+1)r[++num]=n/(n/l),g[num]=(calcg(r[num])-1+mod)%mod,g1[num]=(calcg1(r[num])-1+mod)%mod;
  for(int j=1;j<=cnt;j++){
    t=1ll*prime[j]*prime[j];
    for(int i=num;r[i]>=t;i--){
      g[i]=(g[i]-1ll*prime[j]*(g[id(r[i]/prime[j])]-g[id(prime[j-1])]+mod)%mod+mod)%mod;
      g1[i]=(g1[i]-(g1[id(r[i]/prime[j])]-g1[id(prime[j-1])]+mod)%mod+mod)%mod;
    }
  }
  for(int i=1;i<=num;i++)s[i]=g[i]=(g[i]-g1[i]+(i==1?0:2)+mod)%mod;
  for(int j=cnt;j>=1;j--){
    t=1ll*prime[j]*prime[j];
    for(int i=num;r[i]>=t;i--){
      long long p=prime[j];
      for(int e=1;p*prime[j]<=r[i];e++,p*=prime[j]){
        s[i]=(s[i]+(1ll*(prime[j]^e)*(s[id(r[i]/p)]-g[id(prime[j])]+mod)%mod+(prime[j]^(e+1)))%mod)%mod;
      }
    }
  }
  sort(s+1,s+num+1),num=unique(s+1,s+num+1)-s-1;
  for(int i=1;i<=num;i++)ans^=s[i]+1;
  return cout<<ans<<'\n',0;
}

UOJ188

要求非严格次大质因数的前缀和，考虑一个类 Min_25 筛 DP，设 \(S(n,j)=\sum_{i=1}^n f(i)[\operatorname{lpf}(i)>p_j]\)。

枚举最小质因子和次数 \(p_k,e\)。当除掉这个最小质因数后是个质数，说明答案为 \(p_k\)，贡献为 \(p_k(\pi(\lfloor\frac{n}{p_k^e}\rfloor)-k+1)\)，可以先用 Min_25 求质数个数的块筛；否则除掉 \(p_k^e\) 不影响答案，贡献为 \(S(\lfloor\frac{n}{p_e^k}\rfloor,k)\)。

假如用递推，方程是一样的，相当与把枚举 \(p_k,e\) 变成了前缀和。

#include<bits/stdc++.h>
using namespace std;
int bn,num,prime[316232],cnt;
long long n1,n2,r[632459],g[632459],s[632459],t;
bool vis[316232];
int prime_init(int n,int prime[],bool a[],int cnt=0){
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i;
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0;
      if(i%prime[j]==0)break;
    }
  }
  return cnt;
}
int id(long long n,long long x){
  return x<=bn?x:num-n/x+1;
}
long long solve(long long n){
  num=0,bn=sqrt(n);
  for(int i=1;i<=bn;i++)vis[i]=0;
  cnt=prime_init(bn,prime,vis);
  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;j++){
    t=1ll*prime[j]*prime[j];
    for(int i=num;r[i]>=t;i--)g[i]-=g[id(n,r[i]/prime[j])]-j+1;
  }
  for(int i=1;i<=num;i++)s[i]=0;
  for(int j=cnt;j>=1;j--){
    t=1ll*prime[j]*prime[j];
    for(int i=num;r[i]>=t;i--){
      long long p=prime[j];
      for(int e=1;p*prime[j]<=r[i];e++,p*=prime[j]){
        s[i]+=s[id(n,r[i]/p)]+prime[j]*(g[id(n,r[i]/p)]-j+1);
      }
    }
  }
  return s[num];
}
int main(){
  return cin>>n1>>n2,cout<<solve(n2)-solve(n1-1)<<'\n',0;
}

[[数学]]

posted @ 2024-03-01 09:35  lgh_2009  阅读(25)  评论(0)    收藏  举报