Min_25 筛算法学习笔记

· · 算法·理论

牛客暑期多校连续两场都考了 Min_25 筛，看来不得不学了。

前言

Min_25 筛是一种计算积性函数前缀和的高效算法。

时间复杂度好像是 O\left ( \frac{n^{\frac{3}{4}}}{\log n} \right ) 的。

模板题（洛谷 P5325）

\color{blue}{\texttt{[Problem Discription]}}

已知某积性函数 f(x) 满足 f(p^{k})=p^{k}(p^{k}-1)，其中 p 为质数，k \in \mathbb{N}_{+}。求 \sum\limits_{i=1}^{n} f(i)。

$\color{blue}{\texttt{[Analysis]}}

首先，我们先化简问题，考虑求幂函数 F(p)=p^{k} 的前缀和。

注，对于幂函数而言，F(x) 不仅是积性函数，更是完全积性函数。想一想，怎么证明。

先把求和式子按质数、合数分类：

\sum\limits_{i=1}^{n} F(i)=\sum\limits_{p \text{ is prime, } p \leq n}F(p)+\sum\limits_{i=1}^{n}[i \text{ is not a prime}]F(i)

然后枚举最小质因子：

\sum\limits_{i=1}^{n} F(i)=\sum\limits_{p \text{ is prime, } p \leq n}F(p)+\sum\limits_{1 \leq p \leq \sqrt{n},1 \leq p^{e}\leq n}F(p^{e})\left ( \sum\limits_{i=1}^{\left \lfloor \frac{n}{p^{e}} \right \rfloor} [\text{minp}(i) \geq p]F(i) \right ) 考虑一个天才的 dp。记 $p_{i}$ 表示从小到大的第 $i$ 个素数， $g(n,j)$ 表示前 $n$ 个自然数中满足本身是素数或者最小质因数**大于** $p_{j}$ 的所有整数 $i$ 的 $k$ 次幂的和，即： $$g(n,j)=\sum\limits_{1 \leq i \leq n}i^{k}\left [ i \in \mathbb{P}\text{rime or minp}(i) > p_{j} \right ]$$ 从 $g(n,j-1)$ 转移到 $g(n,j)$ 时，有一些整数不满足条件，需要被剔除。容易发现，被剔除的整数满足**最小质因数**恰好为 $p_{j}$。 可以提出一个公因子 $p_{j}$，剩下的部分就不能有比 $p_{j}$ 更小的质因子了，即 $g\left ( \left \lfloor \frac{n}{p_{j}} \right \rfloor,j \right )-g(p_{j-1},j-1)$，后面那一项是为了把质数全部删除。 这样就得到了 $g$ 的递推关系式： $$g(n,j)=g(n,j-1)-p_{j}^{k} \left [ g\left ( \left \lfloor \frac{n}{p_{j}} \right \rfloor,j \right )-g(p_{j-1},j-1) \right ]$$ > 这里体现出来了完全积性函数的好处，虽然只提取了一个 $p_{j}$，但是函数值可以直接相乘而不需要考虑是否互质。 注意到 $g(p_{j-1},j-1)$ 其实就是前 $(j-1)$ 个质数的 $k$ 次幂的和，且 $p_{j-1} \leq \sqrt{n}$，直接线性筛即可，对应于代码中的 `pre` 数组。 $1$ 到 $n$ 中所有质数的 $k$ 次幂的和其实就是 $g(n,x)$，其中 $p_{x}$ 是**最后一个**小于等于 $\sqrt{n}$ 的质数。为了方便，把 $g(n,x)$ 简单记为 $g(n)$。 为了求出 $g(n)$，我们需要用到以下性质： $$\left \lfloor \frac{\left \lfloor \frac{n}{a} \right \rfloor}{b} \right \rfloor=\left \lfloor \frac{n}{ab} \right \rfloor$$ 因此，无论做了多少次除法，最终得到的数一定可以写成 $\left \lfloor \frac{n}{x} \right \rfloor$ 的形式。可以除法分块。 为了把这 $O(\sqrt{n})$ 个数保存下来，我们需要开两个下标数组，具体看代码。 _____________ 接下来就是求解原问题了。 还是考虑相似的 dp。答案可以分成质数的函数和和合数的函数和两部分。质数的函数和可以**用多项式直接相加减得到**。考虑合数的函数和。 用 $S(n,j)$ 表示 $$S(n,j)=\sum\limits_{i=1}^{n}f(i)[\text{minp}(i)>p_{j}]$$ 则 $$S(n,j)=\color{blue}{g(n)-\text{pre}(x)}\text{\color{black}{+}}\color{red}{\sum\limits_{p_{k}^{e}\leq n,k>x}f(p_{k}^{x}) \left \{ S \left ( \frac{n}{p_{k}^{e}},k \right ) + [e\not = 1] \right \}}$$ 蓝色项是大于 $p_{j}$ 的质数的贡献，红色项是最小质因数大于 $p_{j}$ 的合数的贡献。 具体参看模板题的代码实现。 $\color{blue}{\text{Code}} 
int prime[N],prcnt;
bool is_prime[N];
ll pre1[N],pre2[N];

void get_prime(int n){
    fill(is_prime+2,is_prime+n+1,true);

    for(int i=2;i<=n;i++){
        if (is_prime[i]){
            prime[++prcnt]=i;
            pre1[prcnt]=(pre1[prcnt-1]+i)%mod;
            pre2[prcnt]=(pre2[prcnt-1]+1ll*i*i%mod)%mod;
        }

        for(int j=1;j<=prcnt;j++){
            if (1ll*i*prime[j]>n) break;
            is_prime[i*prime[j]]=false;
            if (i%prime[j]==0) break;
        }
    }
}

ll g1[N<<1],g2[N<<1],w[N<<1];
int ind1[N],ind2[N];

ll n;int m,wcnt;

const int inv2=500000004;
const int inv3=333333336;

void get_g_and_ind(){
    for(ll l=1;l<=n;){
        ll r=n/(n/l);

        w[++wcnt]=n/l;

        int tmp=(n/l)%mod;
        g1[wcnt]=1ll*tmp*(tmp+1)%mod*inv2%mod-1;
        g2[wcnt]=1ll*tmp*(tmp+1)%mod*(2*tmp+1)%mod*inv2%mod*inv3%mod-1;

        if (n/l<=m) ind1[n/l]=wcnt;
        else ind2[r]=wcnt;

        l=r+1;
    }

    for(int i=1;i<=prcnt;i++)
        for(int j=1;j<=wcnt;j++){
            if (1ll*prime[i]*prime[i]>w[j]) break;

            int k=(w[j]/prime[i]<=m)?ind1[w[j]/prime[i]]:ind2[n/(w[j]/prime[i])];

            g1[j]=(g1[j]-1ll*prime[i]*(g1[k]-pre1[i-1]+mod)%mod+mod)%mod;
            g2[j]=(g2[j]-1ll*prime[i]*prime[i]%mod*(g2[k]-pre2[i-1]+mod)%mod+mod)%mod;
        }//这里 g 是可以滚动数组的，最终保留的信息就是 g(n,x)
}

int S(ll x,int y){
    if (prime[y]>=x) return 0;

    int k=(x<=m)?ind1[x]:ind2[n/x];
    int ans=((g2[k]-g1[k]+mod)%mod-(pre2[y]-pre1[y]+mod)%mod+mod)%mod;

    for(int i=y+1;i<=prcnt;i++){
        if (1ll*prime[i]*prime[i]>x) break;

        ll pe=prime[i];
        for(int e=1;pe<=x;e++,pe*=prime[i]){
            int tmp=pe%mod;
            ans=(ans+1ll*tmp*(tmp-1)%mod*(S(x/pe,i)+(e!=1))%mod)%mod;
        }
    }

    return ans;
}

int main(){
    cin>>n;
    m=sqrt(n);

    get_prime(m);
    get_g_and_ind();

    cout<<(S(n,0)+1)%mod;

    return 0;
}

本质理解

虽然 Min_25 筛正宗的用法是用来求积性函数的，但是有些非积性函数的问题也可以用 Min_25 筛来做。

由上面的代码我们可以知道，我们其实不关心 g(n,j) 的值是多少，我们关心的是 g(n,x)，即 g(n)。而 g(n) 里保存的是所有的质数的信息，而且与合数完全无关。

而在 S 的转移式中，蓝色部分求出来的也是函数在所有质数点的取值的和，只有红色部分包含合数，且红色部分只包含合数。

于是，在某些特定的问题中，我们只关心函数在质数的和，或函数在合数的和（两者等价），我们一样可以用 Min_25 筛解决。

一个经典的问题是求某个区间 [l,r] 中质数的个数，可以用前缀和的方法转化为求 [1,n] 中质数的个数。

定义函数

1 \quad n \text{ is a Prime}\\ 0 \quad n \text{ isn't a Prime} \end{cases}

显然这个函数不是积性函数。

但是我们发现 F(n)=1 是积性函数，那么我们完全可以就求 F(n) 的前缀和，不过是求在质数点的函数的和。也就是只要蓝色部分即可。

$\color{blue}{\texttt{[Problem Discription]}}

求 [l,r] 中有多少个数可以加上自己的一个质因数得到一个完全平方数。

$\color{blue}{\texttt{[Analysis]}}

可以证明，满足条件的数一定可以写成如下形式：

x^{2}-p

其中 p 是 x 的一个质因数，且可以证明不同的 x 生成的一定是不同的整数。

转化为求 [1,n]。\left \lfloor \sqrt{n} \right \rfloor+1 是可以生成合法的解的，但是就一个数，特判就好。

对于 1 \dots \sqrt{n} 中的数，其实就是求质因数的个数的和，即

\sum\limits_{1 \leq p \leq \sqrt{n}} \left \lfloor \frac{n}{p} \right \rfloor

对于合适范围内的 p，直接这么干就好了。对于很大的 p，\left \lfloor \frac{n}{p} \right \rfloor 很小，考虑枚举 \left \lfloor \frac{n}{p} \right \rfloor，得到一个区间，求这个区间的质数个数就好。

\color{blue}{\text{Code}} 
const int N=1e5+100;
#define ll long long

int prime[N],prcnt,pre[N];
bool is_prime[N];

void get_prime(int n){
    fill(is_prime+2,is_prime+n+1,true);

    for(int i=2;i<=n;i++){
        if (is_prime[i]){
            prime[++prcnt]=i;
            pre[prcnt]=prcnt;
        }

        for(int j=1;j<=prcnt;j++){
            if (1ll*prime[j]*i>n) break;
            is_prime[i*prime[j]]=false;
            if (i%prime[j]==0) break;
        }
    }
}

int w[N<<1],ind1[N],ind2[N],wcnt;
ll g[N<<1];int n,m;

void get_g_and_ind(){
    for(int l=1;l<=n;){
        int tmp=n/l,r=n/tmp;

        w[++wcnt]=tmp;
        g[wcnt]=tmp-1;

        if (tmp<=m) ind1[tmp]=wcnt;
        else ind2[r]=wcnt;

        l=r+1;
    } 

    for(int i=1;i<=prcnt;i++)
        for(int j=1;j<=wcnt;j++){
            if (1ll*prime[i]*prime[i]>w[j]) break;

            int k=(w[j]/prime[i]<=m)?ind1[w[j]/prime[i]]:ind2[n/(w[j]/prime[i])];
            g[j]-=g[k]-pre[i-1];
        }
} 

ll S(ll x,int y){
    if (prime[y]>=x) return 0;

    int k=(x<=m)?ind1[x]:ind2[n/x];
    return g[k]-pre[y];
}

void init(){
    memset(pre,0,sizeof(pre));
    memset(g,0,sizeof(g));
    memset(ind1,0,sizeof(ind1));
    memset(ind2,0,sizeof(ind2));
    wcnt=prcnt=0;
}

ll Count(int x){
    n=x;
    m=sqrt(n);

    init();
    get_prime(m);
    get_g_and_ind();

    ll res=0;

    for(int i=1;i<=prcnt;i++)
        res+=x/prime[i];

    for(int i=1;i<=m;i++)
        if (x/i>prime[prcnt])
            res+=S(x/i,0)-prcnt;

    return res;
}

ll calc(ll x){
    int r=sqrt(x);
    ll ans=Count(r);

    int tmp=r+1;
    for(int i=2;1ll*i*i<=tmp;i++)
        if (tmp%i==0){
            if (1ll*(r+1)*(r+1)-i<=x) ans++;
            while (tmp%i==0) tmp/=i;
        }

    if (tmp>1&&1ll*(r+1)*(r+1)-tmp<=x) ans++;

    return ans; 
}

int main(){
    ll l,r;
    cin>>l>>r;
    cout<<calc(r)-calc(l-1);

    return 0;
}