Skip to content

Min_25筛 ​

  • 复杂度:O(n34log⁡n)O(\frac{n^{\frac{3}{4}}}{\log n})
  • 质数和合数分开考虑,那么合数中不同最小质因子个数为 nlog⁡n\frac{\sqrt n}{\log n}

定义PiP_i表示第ii个质数,S(n,j)=∑i=2nf(i)(min⁡p∣ip>Pj)S(n,j)=\sum_{i=2}^nf(i)(\min_{p|i} p>P_j)

质数:

  • 假设质数前缀和为 q(n)=∑i=1nf(i)(i∈Prime)q(n)=\sum_{i=1}^nf(i) (i \in \text{Prime})
  • S(n,j)S(n,j)在质数上的和就是q(n)−∑i=1jf(Pi)q(n)-\sum_{i=1}^jf(P_i)

合数:

  • 直接考虑这个数比较抽象,我们换个角度
  • 枚举最小质因数和他的出现次数,那么由于是积性函数,剩余的部分从 ≤n\leq n 变成 ≤nPke\leq \frac{n}{P_k^e}
  • 因此合数的贡献为 ∑k>j∑e=1Pke≤nf(Pke)(S(nPke,k)+[e>1])\sum_{k \gt j} \sum_{e=1}^{P_k^e \leq n}f(P_k^e)(S(\frac{n}{P_k^e},k)+[e>1]),之所以是e>1e>1需要算贡献,我的理解是e=1e=1他就是那么f(Pke)f(P_k^e)就是质数贡献了,和合数贡献不符。e>1e>1需要加贡献是因为S(∗,j)S(*,j)部分是没有包含11的,而f(Pke)f(P_k^e)此时是合数贡献,需要加上

总贡献:

  • S(n,j)=q(n)−∑i=1jf(Pi)+∑k>j∑e=1Pke≤nf(Pke)(S(nPke,k)+[e>1])S(n,j)=q(n)-\sum_{i=1}^jf(P_i)+\sum_{k \gt j} \sum_{e=1}^{P_k^e \leq n}f(P_k^e)(S(\frac{n}{P_k^e},k)+[e>1])

处理质数前缀和:

  • 设完全积性函数f′f'在质数上和ff一致
  • 定义 g(n,j)=∑i=2nf′(i)(i∈Prime∣min⁡p∣i>Pj)g(n,j)= \sum_{i=2}^nf'(i)(i \in \text{Prime}|\min_{p|i}>P_j )
  • 考虑如何转移
  • 假设Pj2>nP_j^2>n,那么g(n,j)=g(n,j−1)g(n,j)=g(n,j-1),因为都是质数的贡献
  • 否则,我们需要减去 min⁡p=Pj\min p=P_j的这些贡献
  • g(n,j)=g(n,j−1)−f′(Pj)(g(nPj,j−1)−∑i=1j−1f′(Pi))g(n,j)=g(n,j-1)-f'(P_j)(g(\frac{n}{P_j},j-1)-\sum_{i=1}^{j-1}f'(P_i)),因为g(nPj,j−1)g(\frac{n}{P_j},j-1)会多算一些小于PjP_j的质数,我们需要减去,才能保证这一部分的最小质数是PjP_j
  • q(n)=g(n,∣P∣)q(n)=g(n,|P|)
  • 对于gg的取值全部都是nk\frac{n}{k}的形式,只有2n2 \sqrt n种取值,记录id1id_1为k≤nk \leq \sqrt n,id2id_2为 nk≤n\frac{n}{k} \leq \sqrt n

总结:min25筛能用的前提是质数处的f(p)f(p)是关于pp的多项式,质数次方的位置f(pe)f(p^e)可以快速计算

常见的pxp^x次方求和公式:

x
1f(n)=n(n+1)2f(n)=\frac{n(n+1)}{2}
2f(n)=n(n+1)(2n+1)6f(n)=\frac{n(n+1)(2n+1)}{6}
3f(n)=[n(n+1)2]2f(n)=[\frac{n(n+1)}{2}]^2
4f(n)=n(n+1)(2n+1)(3n2+3n−1)30f(n)=\frac{n(n+1)(2n+1)(3n^2+3n-1)}{30}
5f(n)=n2(n+1)2(2n2+2n−1)12f(n)=\frac{n^2(n+1)^2(2n^2+2n-1)}{12}
6f(n)=n(n+1)(2n+1)(3n4+6n3−3n2−4n+2)42f(n)=\frac{n(n+1)(2n+1)(3n^4+6n^3-3n^2-4n+2)}{42}
7f(n)=n2(n+1)2(3n4+6n3−n2−4n+2)24f(n)=\frac{n^2(n+1)^2(3n^4+6n^3-n^2-4n+2)}{24}
8f(n)=n(n+1)(2n+1)(5n6+15n5+5n4−15n3−10n2+9n+3)90f(n)=\frac{n(n+1)(2n+1)(5n^6+15n^5+5n^4-15n^3-10n^2+9n+3)}{90}
9f(n)=n2(n+1)2(2n6+6n5−5n4−15n3+10n2+12n−6)20f(n)=\frac{n^2(n+1)^2(2n^6+6n^5-5n^4-15n^3+10n^2+12n-6)}{20}
10f(n)=n(n+1)(2n+1)(n2+n−1)(3n6+9n5+2n4−11n3+3n2+10n−5)66f(n)=\frac{n(n+1)(2n+1)(n^2+n-1)(3n^6+9n^5+2n^4-11n^3+3n^2+10n-5)}{66}

实际上只有前三个是常见的

upd板子使用步骤:

  • 假设当前要求的是一个积性函数g(x)g(x)
  • 那么我们需要找到若干个完全积性函数,并且都能O(1)O(1)求和,使得g(x)=g0(x)+g1(x)+...+g(x)=g0(x)+g1(x)+...+,当xx是质数时必须相等
  • 然后就是质数幂处的数能O(1)O(1)求出来
  • 如上面的积性函数为f(x)=x(x−1)f(x)=x(x-1),在质数幂处为f(pk)=pk(pk−1)=x(x−1)f(p^k)=p^k(p^k-1)=x(x-1),可以拆分成g0(x)=x2g0(x)=x^2和 g1(x)=xg1(x)=x 两个完全积性函数
cpp
namespace Min25{
    // sum[0]质数前缀和,sum[1]质数平方前缀和
    // g为g1-g0, g0:f'(k)=k, g1:f'(k)=k^2
    // 将积性函数g拆分,注意g1,g0必须满足为完全积性函数即f(p)=p^k
    int prime[N],cnt,sum[2][N],g[2][N];
    int w[N],tot,id1[N],id2[N],sq,n;
    bool vis[N];
    void sieve(int n=N+1){
        int c=0;
        for(int i=2;i<=n;i++){
            if(!vis[i]){
                vis[i]=1;
                prime[++cnt]=i;
                //======= need to be modified =======//
                sum[0][cnt]=(sum[0][cnt-1]+i)%mod;//f(x)=x的前缀和
                sum[1][cnt]=(sum[1][cnt-1]+i*i%mod)%mod;//f(x)=x^2的前缀和
                //============== END ==============//
            }
            for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
                vis[i*prime[j]]=1;
                if(i%prime[j]==0)break;
            }
        }
    }   
    int qpow(int a,int b){
        int ans=1;
        for(;b;b>>=1){
            if(b&1)ans=ans*a%mod;
            a=a*a%mod;
        }
        return ans;
    }
    int getid(int x){
        if(x<=sq)return id1[x];
        return id2[n/x];
    }
    static int inv6=qpow(6,mod-2),inv2=qpow(2,mod-2);
    //====== need to be modified ======//
    int f0(int x){// 1~x的前缀和 
        x%=mod;
        return x*(x+1)/2%mod;
    }
    int f1(int x){// 1~x的前缀平方和
        x%=mod;
        return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
    }
    //=========== END ==================//
    void init(int x){// 初始化下标的问题,并且预处理出要用到的完全积性函数
        n=x;
        sq=sqrtl(n);
        sieve(sq);
        for(int l=1,r;l<=n;l=r+1){
            r=n/(n/l);
            w[++tot]=n/l;
            //======= need to be modified =======//
            g[0][tot]=f0(w[tot])-1;//f(x)=x^2对应的g(n/l,0)
            g[1][tot]=f1(w[tot])-1;//f(x)=x对应的g(n/l,0)
            //============= END ==================//
            if(n/l<=sq)id1[n/l]=tot;//映射
            else id2[n/(n/l)]=tot;
        }
        //======= need to be modified ======//
        for(int i=1;i<=cnt;i++){
            for(int j=1;j<=tot&&prime[i]*prime[i]<=w[j];j++){
                g[0][j]=(g[0][j]-prime[i]*(g[0][getid(w[j]/prime[i])]-sum[0][i-1]+mod)%mod+mod)%mod;
                g[1][j]=(g[1][j]-prime[i]*prime[i]%mod*(g[1][getid(w[j]/prime[i])]-sum[1][i-1]+mod)%mod+mod)%mod;
            }
        }
        //========== END =================//
    }
    //====== need to be modified ======//
    int F(int x){
        x%=mod;
        return x*(x-1)%mod; 
    }
    //=========== END =====================//
    int S(int x,int y){
        if(prime[y]>=x)return 0;
        int k=getid(x);// 获取下标
        //====== need to be modified ======//
        int res=(g[1][k]-g[0][k]+mod-sum[1][y]+mod+sum[0][y])%mod;//q(n)=g[1][k]-g[0][k],\sum_{i=1}^yf(i)=sum[1][y]-sum[0][y]
        //=========== END ==================//
        for(int i=y+1;i<=cnt&&prime[i]*prime[i]<=x;i++){// 枚举p_i
            for(int e=1,p=prime[i];p<=x;p=p*prime[i],e++){// 枚举p_i^e
                res=(res+F(p)*(S(x/p,i)%mod+(e!=1)))%mod;
            }
        }
        return res;
    }
};

Min25筛modint版 ​

cpp
struct mint {
    // constexpr static int mod = 998244353;
    constexpr static int mod = 1e9+7;
	int x;
	inline mint(int o = 0) { x = o%mod; }
	inline mint & operator = (int o) { return x = o%mod, *this; }
	inline mint & operator += (mint o) { return (x += o.x) >= mod && (x -= mod), *this; }
	inline mint & operator -= (mint o) { return (x -= o.x) < 0 && (x += mod), *this; }
	inline mint & operator *= (mint o) { return x = (long long) x * o.x % mod, *this; }
    inline mint& qpow(int b){
        mint w = *this;
        mint ret(1);
        for(; b; b >>= 1, w *= w) if(b & 1) ret *= w;
        return x = ret.x, *this;
    }
	inline mint & operator /= (mint o) { return *this *= (o .qpow( (mod - 2))); }
	friend inline mint operator + (mint a, mint b) { return a += b; }
	friend inline mint operator - (mint a, mint b) { return a -= b; }
	friend inline mint operator * (mint a, mint b) { return a *= b; }
	friend inline mint operator / (mint a, mint b) { return a /= b; }
    friend ostream& operator<<(ostream& os, const mint& m) {
        return os << m.x;
    }
};
template<typename T>
struct Min25{
    // sum[0]质数前缀和,sum[1]质数平方前缀和
    // g为g1-g0, g0:f'(k)=k, g1:f'(k)=k^2
    // 将积性函数g拆分,注意g1,g0必须满足为完全积性函数即f(p)=p^k
    int prime[N],w[N],cnt,tot,id1[N],id2[N],sq,n;
    T sum[2][N],g[2][N];
    bool vis[N];
    void sieve(int n=N+1){
        int c=0;
        for(int i=2;i<=n;i++){
            if(!vis[i]){
                vis[i]=1;
                prime[++cnt]=i;
                //======= need to be modified =======//
                sum[0][cnt]=sum[0][cnt-1]+i;//f(x)=x的前缀和
                sum[1][cnt]=sum[1][cnt-1]+i*i;//f(x)=x^2的前缀和
                //============== END ==============//
            }
            for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
                vis[i*prime[j]]=1;
                if(i%prime[j]==0)break;
            }
        }
    }   
    int getid(int x){
        if(x<=sq)return id1[x];
        return id2[n/x];
    }
    T inv6=T(6).qpow(mod-2),inv2=T(2).qpow(mod-2);
    //====== need to be modified ======//
    T g0(T x){// g0 
        return x;
    }
    T g1(T x){// g1
        return x*x;
    }
    T f0(T x){// 1~x的前缀和 
        return x*(x+1)*inv2;
    }
    T f1(T x){// 1~x的前缀平方和
        return x*(x+1)*(2*x+1)*inv6;
    }
    //=========== END ==================//
    void init(int x){// 初始化下标的问题,并且预处理出要用到的完全积性函数
        n=x;
        sq=sqrtl(n);
        sieve(sq);
        for(int l=1,r;l<=n;l=r+1){
            r=n/(n/l);
            w[++tot]=n/l;
            //======= need to be modified =======//
            g[0][tot]=f0(w[tot])-T(1);//f(x)=x^2对应的g(n/l,0)
            g[1][tot]=f1(w[tot])-T(1);//f(x)=x对应的g(n/l,0)
            //============= END ==================//
            if(n/l<=sq)id1[n/l]=tot;//映射
            else id2[n/(n/l)]=tot;
        }
        //======= need to be modified ======//
        for(int i=1;i<=cnt;i++){
            for(int j=1;j<=tot&&prime[i]*prime[i]<=w[j];j++){
                g[0][j]-=g0(prime[i])*(g[0][getid(w[j]/prime[i])]-sum[0][i-1]);
                g[1][j]-=g1(prime[i])*(g[1][getid(w[j]/prime[i])]-sum[1][i-1]);
            }
        }
        //========== END =================//
    }
    //====== need to be modified ======//
    T F(T x){
        return x*(x-1); 
    }
    //=========== END =====================//
    T S(int x,int y){
        if(prime[y]>=x)return 0;
        int k=getid(x);// 获取下标
        //====== need to be modified ======//
        T res=g[1][k]-g[0][k]-sum[1][y]+sum[0][y];//q(n)=g[1][k]-g[0][k],\sum_{i=1}^yf(i)=sum[1][y]-sum[0][y]
        //=========== END ==================//
        for(int i=y+1;i<=cnt&&prime[i]*prime[i]<=x;i++){// 枚举p_i
            for(int e=1,p=prime[i];p<=x;p=p*prime[i],e++){// 枚举p_i^e
                res+=F(p)*(S(x/p,i)+(e!=1));
            }
        }
        return res;
    }
};
Min25<mint> sieve;
void solve(){
    int n;
    cin>>n;
    sieve.init(n);
    cout<<sieve.S(n,0)+mint(1)<<"\n";
}

欧拉函数和莫比乌斯函数求和 ​

欧拉函数质数幂次表示:ϕ(pk)=pk(1−1p)=pk−pk−1\phi(p^k)=p^k(1-\frac{1}{p})=p^k-p^{k-1}可以O(1)O(1)算

构造完全积性函数g0(x)=1,g1(x)=x,g(x)=g1(x)−g0(x)=x−1g0(x)=1,g1(x)=x,g(x)=g1(x)-g0(x)=x-1,在质数处相等

然后套板子即可

对于莫比乌斯函数质数次幂表示: μ(pk)=(−1)[k=1]\mu(p^k)=(-1)[k=1]

莫比乌斯本身就是完全积性函数,不用构造,而且求解g函数部分比较特殊

cpp
namespace Min25{
    int prime[N],cnt,sum[2][N],g[2][N];
    int w[N],tot,id1[N],id2[N],sq,n;
    bool vis[N];
    int sum_mu[N],g_mu[N];
    void sieve(int n=N+1){
        for(int i=2;i<=n;i++){
            if(!vis[i]){
                vis[i]=1;
                prime[++cnt]=i;
                //======= need to be modified =======//
                sum[0][cnt]=sum[0][cnt-1]+1;// 质数个数前缀和
                sum[1][cnt]=sum[1][cnt-1]+i;// 质数前缀和
                sum_mu[cnt]=sum_mu[cnt-1]-1;
                //============== END ==============//
            }
            for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
                vis[i*prime[j]]=1;
                if(i%prime[j]==0)break;
            }
        }
    }   
    int getid(int x){
        if(x<=sq)return id1[x];
        return id2[n/x];
    }
    //====== need to be modified ======//
    int f0(int x){ 
        return x;
    }
    int f1(int x){
        return x*(x+1)/2;
    }
    //=========== END ==================//
    void init(int x){// 初始化下标的问题,并且预处理出要用到的完全积性函数
        n=x;
        sq=sqrtl(n);
        sieve(sq);
        for(int l=1,r;l<=n;l=r+1){
            r=n/(n/l);
            w[++tot]=n/l;
            //======= need to be modified =======//
            g[0][tot]=f0(w[tot])-1;
            g[1][tot]=f1(w[tot])-1;
            g_mu[tot]=w[tot]>1?-(w[tot]-1):0;
            //============= END ==================//
            if(n/l<=sq)id1[n/l]=tot;//映射
            else id2[n/(n/l)]=tot;
        }
        //======= need to be modified ======//
        for(int i=1;i<=cnt;i++){
            for(int j=1;j<=tot&&prime[i]*prime[i]<=w[j];j++){
                int k=getid(w[j]/prime[i]);
                g[0][j]-=g[0][k]-sum[0][i-1];
                g[1][j]-=prime[i]*(g[1][k]-sum[1][i-1]);
                g_mu[j]+=sum_mu[i-1]-g_mu[k];
            }
        }
        //========== END =================//
    }
    int S(int x,int y){
        if(prime[y]>=x)return 0;
        int k=getid(x);// 获取下标
        //====== need to be modified ======//
        int res=g[1][k]-g[0][k]-sum[1][y]+sum[0][y];
        //=========== END ==================//
        for(int i=y+1;i<=cnt&&prime[i]*prime[i]<=x;i++){// 枚举p_i
            for(int e=1,p=prime[i];p<=x;p=p*prime[i],e++){// 枚举p_i^e
                int val=p-p/prime[i];// phi(p_i^e)=p_i^e-p_i^(e-1)
                res+=val*(S(x/p,i)+(e!=1));
            }
        }
        return res;
    }
    int D(int x,int y){
        if(prime[y]>=x)return 0;
        int k=getid(x);// 获取下标
        int res=g_mu[k]-sum_mu[y];
        for(int i=y+1;i<=cnt&&prime[i]*prime[i]<=x;i++){// 枚举p_i
            res-=D(x/prime[i],i);
        }
        return res;
    }
    void clear(){
        for(int i=0;i<=cnt;i++)sum[0][i]=sum[1][i]=sum_mu[i]=prime[i]=0;
        for(int i=0;i<=sq;i++)vis[i]=id1[i]=id2[i]=0;
        for(int i=0;i<=tot;i++)w[i]=g[0][i]=g[1][i]=g_mu[i]=0;
        cnt=tot=0;
    }
};