筛法

Fluorine Wang

2024-05-12

ACM

Math

区间筛

例题:每次给一对 $l,r$ ,求区间内的质数有多少个.
范围是 $l\le r\le10^{12},b-a\le10^6$.

首先要明白一个常见的道理: 区间 $[a,b]$ 内最大的质因数是 $\sqrt b$ .所以我们要先打好1e6内的表,然后利用坐标偏移再打好lr的表即可.

看到 $10^{12}$ 的数据范围要想到区间筛.因为数据范围这么大的话,一般要用到,”一个数的最小质因数不会超过 $\sqrt n$ “这个结论,然后先质数筛,再对每个数进行质因数分解O(logn).然后利用唯一分解定理把约数组合拼上去即可.

Min_25筛

Min_25筛用于解决积性函数前缀和的问题(事实上不是积性函数也能写,那是后话),因为常数小,跑的比杜教筛还快.

Min_25构成主要是一个神仙dp,原理是下面四个公式,代码略,直接抄题解,这不是这篇文章的重点.

构造一个完全积性函数

$F(x)=f(x)$

只需要满足在质数点和 $f(x)$ 一样即可,质数次方不需要满足性质.

$g(n,j)=\sum_{i=1}^nF(i)[i\in p||i的最小质因子大于第j个质数]$ $g(n,j)=\begin{cases}g(n,j-1)&p_j^2>n\\ g(n,j-1)-F(p_j)\cdot\left(g\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right)-\sum_{i=1}^{j-1}F(p_{j})\right)&p_j^2\le n\end{cases}$ $S(n,j)=\sum_{i=1}^nf(i)[i的最小质因子大于第j个质数]$ $S(n,j)=g(n,|P|)-\sum_{i=1}^{j-1}f(p_i)+\sum_{k=j+1}^{p_k^2\le n}\sum^{p_k^e\le n}_{e=1}f(p_k^e)\left(S\left(\left\lfloor\frac n{p_k^e}\right\rfloor,k\right)+[e>1]\right)$

本篇侧重于对Min_25筛中完全积性函数的构造.

数论分块的结果是有 $2\sqrt n$ 个取值,所以对于 $1e10$ 大小的数据要开 $2e5$ 大小的数组,,,对于1e11要开7e5大小的数组,直接1e6准没错,再次强调.

要求

$f$ 是积性函数
$f(p^k)$ 能快速计算

大致流程

先把质数点算出来看是什么,即 $f(p)=xxx$ ,得到一个多项式.
然后把 $f(p)$ 按照多项式的形式拆开,比如 $f(p)=p-1$ 会被拆成 $g_1(p)=p$ 和 $g_2(p)=1$ .
然后对拆出来的g分别求和(这个时候不要管正负号): $sum=\sum_{i=2}^ng(i)$ .
剩下的全是Min_25板子,不过需要一点小改动:

正负号的作用体现在 S(n) 的计算中,也就是 $res=\pm g_i(x)$ 在这里体现正负号和倍数(如果有)
函数拆分的作用体现在 solve(n) 函数中,具体地,只需要修改后面两个for循环即可.
一定要注意在算g函数的时候使用的质数应该是 p[i] 不是 p[j] ,Flu在调欧拉函数因为这里没注意到调了俩小时.

luogu P5325

$f(p^k)=p^k(p^k-1)$

这个东西能拆成两部分,前后两部分都是完全积性的,计算的时候拼一起就行了.

long long sq,n;
const int mod=1000000007;
long long id1[200010],id2[200010],w[200010];
int   vis[200010];//存最小质因数,负的表示质数表中的位置(负的)
long long p[200010],ptop=0;//存质数 
void sieve(int n){//[1,n]
    int tmp;
    for(int i=2;i<=n;++i){
        if(!vis[i]){
            vis[i]=-(++ptop);
            p[ptop]=i;
        }
        for(int j=1;j<=ptop&&i*p[j]<=n;++j){
            vis[i*p[j]]=p[j];
            if(i%p[j]==0){
                break;
            }
        }
    }
}
template<typename T1,typename T2,typename T3>
T1 qp(T1 b,T2 po,T3 p){
    T1 res(1);
    while(po>0){
        if(po&1)
            res=res*b%p;
        b=b*b%p;
        po>>=1;
    }
    return res;
}
long long getid(long long x){
    if(x<=sq)return id1[x];
    return id2[n/x];
}
long long m,inv2,inv6;
long long g1[200010],g2[200010],sum1[200010],sum2[200010];

long long f1(long long x){
    x%=mod;
    return x*(x+1)%mod*inv2%mod;
}
long long f2(long long x){
    x%=mod;
    return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
}
long long S(long long n,long long j){
    if(p[j]>n)return 0;
    i64 res=(g2[getid(n)]-g1[getid(n)]-sum2[j]+sum1[j])%mod;
    for(int k=j+1;k<=ptop&&p[k]*p[k]<=n;++k){
        for(i64 e=1,pk=p[k];pk<=n;++e,pk*=p[k]){
            res=(res+(pk%mod*(pk%mod-1)%mod*(S(n/pk,k)+(e>1))%mod))%mod;
        }
    }
    return (res%mod+mod)%mod;
}
long long solve(long long n){
    sq=sqrt(n),inv2=qp((long long)2,mod-2,mod),inv6=qp((long long)6,mod-2,mod);
    sieve(sq);
    int m=0;
    for(i64 i=1;i<=ptop;++i){
        sum1[i]=(sum1[i-1]+p[i])%mod;
        sum2[i]=(sum2[i-1]+p[i]*p[i])%mod;
    }
    for(i64 l=1,r;l<=n;l=r+1){
        r=(n/(n/l));
        w[++m]=n/l;
        g1[m]=f1(w[m])-1;
        g2[m]=f2(w[m])-1;
        if(w[m]<=sq)id1[w[m]]=m;
        else id2[n/w[m]]=m;
    }
    for(int i=1;i<=ptop;++i){
        for(int j=1;j<=m&&p[i]*p[i]<=w[j];++j){
            g1[j]=(g1[j]-p[i]*(g1[getid(w[j]/p[i])]-sum1[i-1])%mod)%mod;
            g2[j]=(g2[j]-p[i]*p[i]%mod*(g2[getid(w[j]/p[i])]-sum2[i-1])%mod)%mod;
        }
    }
    return (S(n,0)+1)%mod;
}

$\sum_{i=1}^n\mu(i)$

拆: $\mu(p)=-1,g(p)=-1,sum=-n+1$ .代码:

long long sq,n;
const int mod=1000000007;
long long id1[200010],id2[200010],w[200010];
int   vis[1000010];//存最小质因数,负的表示质数表中的位置(负的)
long long  p[100010],ptop=0;//存质数 
void sieve(int n){//[1,n]
    int tmp;
    for(int i=2;i<=n;++i){
        if(!vis[i]){
            vis[i]=-(++ptop);
            p[ptop]=i;
        }
        for(int j=1;j<=ptop&&i*p[j]<=n;++j){
            vis[i*p[j]]=p[j];
            if(i%p[j]==0){
                break;
            }
        }
    }
}
template<typename T1,typename T2,typename T3>
T1 qp(T1 b,T2 po,T3 p){
    T1 res(1);
    while(po>0){
        if(po&1)
            res=res*b%p;
        b=b*b%p;
        po>>=1;
    }
    return res;
}
long long getid(long long x){
    if(x<=sq)return id1[x];
    return id2[n/x];
}

long long g[200010],sum[200010];

long long f(long long x){//快速计算构造的完全积性函数前缀和
    return -x+1;
}

long long S(long long n,long long j){
    if(p[j]>n)return 0;
    i64 res=(g[getid(n)]-sum[j])%mod;
    for(int k=j+1;k<=ptop&&p[k]*p[k]<=n;++k){
        for(i64 e=1,pk=p[k];pk<=n;++e,pk*=p[k]){
            res=(res+(/*f(p^e_k)*/(e>1?0:-1)*(S(n/pk,k)+(e>1))%mod))%mod;
        }
    }
    return (res%mod+mod)%mod;
}

int main(){
    fin>>n;
    sq=sqrt(n);
    sieve(sq);
    int m=0;
    for(i64 i=1;i<=ptop;++i){
        sum[i]=(sum[i-1]-1)%mod;
    }
    for(i64 l=1,r;l<=n;l=r+1){
        r=(n/(n/l));
        w[++m]=n/l;
        g[m]=f(w[m]);
        if(w[m]<=sq)id1[w[m]]=m;
        else id2[n/w[m]]=m;
    }
    for(int i=1;i<=ptop;++i){
        for(int j=1;j<=m&&p[i]*p[i]<=w[j];++j){
            g[j]=(g[j]-/*f(p_j)*/1*(g[getid(w[j]/p[i])]-sum[i-1])%mod)%mod;
        }
    }
    fout<<(S(n,0)+1)%mod;
    return 0;
}

关于为什么F(x)=1:
我们对于常数项构建的函数实际上是F(x)=1,完全积性函数,只是在匹配f(x)的时候我们加上了常数项的偏移(也就是-1).

$\sum_{i=1}^n\varphi(i)$

$f(p)=p-1$ ,拆成两部分: $g1(p)=p,g2(p)=1$ ,然后分别计算即可,提交记录

long long sq,n;
const int mod=998244353;
long long id1[200010],id2[200010],w[200010];
int   vis[1000010];//存最小质因数,负的表示质数表中的位置(负的)
long long  p[100010],ptop=0;//存质数 
void sieve(int n){//[1,n]
    int tmp;
    for(int i=2;i<=n;++i){
        if(!vis[i]){
            vis[i]=-(++ptop);
            p[ptop]=i;
        }
        for(int j=1;j<=ptop&&i*p[j]<=n;++j){
            vis[i*p[j]]=p[j];
            if(i%p[j]==0){
                break;
            }
        }
    }
}
template<typename T1,typename T2,typename T3>
T1 qp(T1 b,T2 po,T3 p){
    T1 res(1);
    while(po>0){
        if(po&1)
            res=res*b%p;
        b=b*b%p;
        po>>=1;
    }
    return res;
}
long long getid(long long x){
    if(x<=sq)return id1[x];
    return id2[n/x];
}


long long g1[200010],g2[200010],sum1[200010],sum2[200010];

long long f1(long long x){
    x%=mod;
    return ((x*(x+1)/2)-1)%mod;
}
long long f2(long long x){//快速计算构造的完全积性函数前缀和
    x%=mod;
    return x-1;
}

long long S(long long n,long long j){
    if(p[j]>n)return 0;
    i64 res=((g1[getid(n)]-sum1[j])-(g2[getid(n)]-sum2[j]))%mod;
    for(int k=j+1;k<=ptop&&p[k]*p[k]<=n;++k){
        for(i64 e=1,pk=p[k];pk<=n;++e,pk*=p[k]){
            res=(res+(/*f(p^e_k)*/(pk-pk/p[k])%mod*(S(n/pk,k)+(e>1))%mod))%mod;
        }
    }
    return (res%mod+mod)%mod;
}


long long solve(long long n){
    sq=sqrt(n);
    sieve(sq);
    int m=0;
    for(i64 i=1;i<=ptop;++i){
        sum1[i]=(sum1[i-1]+p[i])%mod;
        sum2[i]=(sum2[i-1]+1)%mod;
    }
    for(i64 l=1,r;l<=n;l=r+1){
        r=(n/(n/l));
        w[++m]=n/l;
        g1[m]=f1(w[m]);
        g2[m]=f2(w[m]);
        if(w[m]<=sq)id1[w[m]]=m;
        else id2[n/w[m]]=m;
    }
    for(int i=1;i<=ptop;++i){
        for(int j=1;j<=m&&p[i]*p[i]<=w[j];++j){
            g1[j]=(g1[j]-/*f1(p_j)*/p[i]*(g1[getid(w[j]/p[i])]-sum1[i-1])%mod)%mod;
            g2[j]=(g2[j]-/*f2(p_j)*/(1)*(g2[getid(w[j]/p[i])]-sum2[i-1])%mod)%mod;
        }
    }
    return (S(n,0)+1)%mod;
}

int main(){
    cin>>n;
    cout<<solve(n);
    return 0;
}

$\pi(n)=\sum_{i=1}^n[i\in prime]$

这个题略微有点不同,并未让你构造完全积性函数,只是让你寻找g函数的dp本来的意义.

long long sq,n;
// const int mod=1000000007;
long long id1[1000010],id2[1000010],w[1000010];
int   vis[1000010];//存最小质因数,负的表示质数表中的位置(负的)
long long  p[1000010],ptop=0;//存质数 
void sieve(int n){//[1,n]
    for(int i=2;i<=n;++i){
        if(!vis[i]){
            vis[i]=-(++ptop);
            p[ptop]=i;
        }
        for(int j=1;j<=ptop&&i*p[j]<=n;++j){
            vis[i*p[j]]=p[j];
            if(i%p[j]==0){
                break;
            }
        }
    }
}
long long getid(long long x){
    return (x<=sq)?id1[x]:id2[n/x];
}
long long g[1000010],sum[1000010];
long long f(long long x){//快速计算构造的完全积性函数前缀和
    return x-1;
}
long long solve(long long n){
    sq=sqrt(n);
    sieve(sq);
    int m=0;
    for(i64 i=1;i<=ptop;++i){
        sum[i]=(sum[i-1]+1);
    }
    for(i64 l=1,r;l<=n;l=r+1){
        r=(n/(n/l));
        w[++m]=n/l;
        g[m]=f(w[m]);
        if(w[m]<=sq)id1[w[m]]=m;
        else id2[n/w[m]]=m;
    }
    for(int i=1;i<=ptop;++i){
        for(int j=1;j<=m&&p[i]*p[i]<=w[j];++j){
            g[j]=(g[j]-/*f(p_j)*/(g[getid(w[j]/p[i])]-sum[i-1]));
        }
    }
    return g[1];
}