min_25筛学习笔记

min_25筛是一个能快速求解积性函数前缀和的东西。

要保证 $f(p^k)(p\text{是质数})$ 是个关于 $p^k$ 的多项式。(项数也不要太多)

以下以洛谷的模板为例:($f(p^k)=(p^k)^2-p^k(p\text{是质数})$,求前 $N(N\le 10^{10})$ 项的和 $\bmod 10^9+7$ 的值)

也就是 $\sum\limits^N_{i=1}f(i)$。

我们考虑分成质数和合数计算,同时枚举合数的最小质因子及其次数(因为是积性函数,所以直接提出来)。

答案为 $\sum\limits^N_{p=1}f(p)[p\in prime]+\sum p^e[p\in prime,p^e\le N](\sum\limits_{i=1}^{\lfloor\frac{N}{p^e}\rfloor}f(i)[i\text{的最小质因子}>p])$


首先考虑质数部分:

我们对每个次数不同的项分别计算。假设我们正在计算 $k$ 次方项。

我们不知道为什么但就是设 $g(n,j)$ 为 $\sum\limits^n_{i=1}i^k[i\text{是质数或者}i\text{的最小质因子}>p_j]$。

其中 $p_j$ 是第 $j$ 小的质数。

那么有初始状态 $g(n,0)=\sum\limits^n_{i=2}i^k$(不把 $1$ 算进去)。

转移,发现不满足 $g(n,j)$ 的限制且满足 $g(n,j-1)$ 的限制的数就是最小质因子恰为 $p_j$ 的合数。

首先发现 $p_j^2>n$ 时这样的数不存在,此时 $g(n,j)=g(n,j-1)$。

否则由于 $k$ 次方是完全积性函数,所以可以把 $p_j^k$ 提出来,那么这些数的 $k$ 次方和就是 $p_j^k(g(\lfloor\dfrac{n}{p_j}\rfloor,j)-g(p_{j-1},j-1))$。

为什么前面的 $g$ 第二维是 $j$ 呢?因为我们只提出了一个 $p_j$,剩下可能还有。

为什么要减掉后面的 $g$ 呢?因为后面是 $<p_j$ 的质数的 $k$ 次方和,如果把这一部分算上,那么就相当于整个答案多算了一些最小质因子小于 $p_j$ 的数。所以要减掉。(其实这一块就是质数 $k$ 次方的前缀和,可以预处理。以后记作 $s_k(j-1)$)

最后用 $g(n,j-1)$ 减掉不满足 $g(n,j)$ 的限制且满足 $g(n,j-1)$ 的限制的数即可。

方程出来了:

$$g(n,j)=\begin{cases}g(n,j-1)&p_j^2>n\\g(n,j-1)-p_j^k(g(\lfloor\dfrac{n}{p_j}\rfloor,j)-s_k(j-1))&p_j^2\le n\end{cases}$$

这个可以上滚动数组滚掉 $j$,同时对于 $p_j^2>n$ 的不用动就行了,时间复杂度也降到了……等等,这玩意复杂度明明比暴力还高好吗!

但是我们可以发现,一个 $g(n,j)$ 可以由若干个 $g(\lfloor\dfrac{n}{x}\rfloor,k)$ 表示。

这启示我们把所有 $\lfloor\dfrac{N}{n}\rfloor$ 相同的 $n$ 归为一类一起处理。具体实现,把 $g$ 的第一维下标变为 $\lfloor\dfrac{N}{n}\rfloor$。然而这样还是可能太大,于是可以对 $\le\sqrt{n}$ 和 $>\sqrt{n}$ 的开两种编号。

这一部分时间复杂度约为 $O(\frac{N^{3/4}}{\log N})$。

怎么算的?我们发现类似埃氏筛,一个数 $x$ 被更新过 $\frac{\sqrt{x}}{\log\sqrt{x}}$ 次。

那么复杂度 $O(\sum\limits_{i=1}^{\sqrt{N}}\frac{\sqrt{i}}{\log\sqrt{i}}+\sum\limits_{i=1}^{\sqrt{N}}\frac{\sqrt{\frac{N}{i}}}{\log\sqrt{i}})$。经过积分爆算,发现略小于 $O(\frac{N^{3/4}}{\log N})$。

先放一下模板的代码实现:

min_25筛学习笔记
inline int id(ll x){return x<=sq?id1[x]:id2[n/x];}
void init(){
    sq=sqrt(n);
    FOR(i,2,sq){
        if(!vis[i]){
            pri[++pl]=i;
            s1[pl]=(s1[pl-1]+i)%mod;
            s2[pl]=(s2[pl-1]+(ll)i*i)%mod;    //质数一次项和二次项前缀和 
        }
        FOR(j,1,pl){
            if(i*pri[j]>sq) break;
            vis[i*pri[j]]=true;
            if(i%pri[j]==0) break;
        }
    }
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        w[++tot]=n/l;
        if(n/l<=sq) id1[n/l]=tot;
        else id2[n/(n/l)]=tot;    //设置编号 
        int x=w[tot]%mod;
        g1[tot]=(1ll*x*(x+1)%mod*inv2%mod-1+mod)%mod; 
        g2[tot]=(1ll*x*(x+1)%mod*(2*x+1)%mod*inv6%mod-1+mod)%mod;    //2到x的和 
    }
}
void calc_g(){
    FOR(i,1,pl){    //遍历质数(上文中的j) 
        FOR(j,1,tot){    //遍历上文中的n 
            if((ll)pri[i]*pri[i]>w[j]) break;    //只管p_j^2<=n的 
            g1[j]=(g1[j]-1ll*pri[i]*(g1[id(w[j]/pri[i])]-s1[i-1]+mod)%mod+mod)%mod;
            g2[j]=(g2[j]-1ll*pri[i]*pri[i]%mod*(g2[id(w[j]/pri[i])]-s2[i-1]+mod)%mod+mod)%mod;
            //上面这两行就死磕吧……记得g1,g2和w下标的含义就好 
        }
    }
}
View Code

再看整个式子:

我们吸取了经验知道要绞尽脑汁地设 $S(n,j)$ 表示 $\sum\limits_{i=2}^n f(i)[i\text{的最小质因子}>p_j]$。注意这里是 $f$ 而不是 $k$ 次方,下界是 $2$。

那么要求即为 $S(N,0)+f(1)$。边界 $S(n,j)=0(p_j>n)$。

转移,对答案有贡献的质数必定大于 $p_j$,那么贡献就是 $g(n,pl)-\sum s_k(j)$。($pl$ 是 $\le\sqrt{N}$ 的质数个数)

按照套路再枚举合数的最小质因子和次数,就是 $\sum p^e[p\in prime,p>p_j,p^e\le n](S(\lfloor\dfrac{n}{p^e}\rfloor)+[e\ne 1])$。

(为何后面要有 $[e\ne 1]$ 呢?因为 $e=1$ 时,这个值在前面算质数时已经算过,而如果 $e\ne 1$,那么前面算质数时没算过,而因为 $S$ 下界为 $2$ 后面也不会算到,所以要加上 $f(p^e)$)。

即使不记忆化,这一部分复杂度也是 $O(\frac{N^{3/4}}{\log N})$。(这个复杂度我是真不会分析了)

贴上整题代码:

min_25筛学习笔记
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=1000100,mod=1000000007,inv2=500000004,inv6=166666668;
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define ROF(i,a,b) for(int i=(a);i>=(b);i--)
#define MEM(x,v) memset(x,v,sizeof(x))
int sq,pri[maxn],pl,tot,g1[maxn],g2[maxn],id1[maxn],id2[maxn],s1[maxn],s2[maxn];
bool vis[maxn];
ll n,w[maxn];
inline int add(int a,int b){return a+b<mod?a+b:a+b-mod;}
inline int sub(int a,int b){return a<b?a-b+mod:a-b;}
inline int mul(int a,int b){return (ll)a*b%mod;}
inline int id(ll x){return x<=sq?id1[x]:id2[n/x];}
void init(){
    sq=sqrt(n);
    FOR(i,2,sq){
        if(!vis[i]) pri[++pl]=i,s1[pl]=add(s1[pl-1],i),s2[pl]=add(s2[pl-1],mul(i,i));
        FOR(j,1,pl){
            if(i*pri[j]>sq) break;
            vis[i*pri[j]]=true;
            if(i%pri[j]==0) break;
        }
    }
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        w[++tot]=n/l;
        if(n/l<=sq) id1[n/l]=tot;
        else id2[n/(n/l)]=tot;
        g1[tot]=sub(mul(mul(w[tot]%mod,(w[tot]+1)%mod),inv2),1);
        g2[tot]=sub(mul(mul(mul(w[tot]%mod,(w[tot]+1)%mod),(2*w[tot]+1)%mod),inv6),1);
    }
}
void calc_g(){
    FOR(i,1,pl){
        FOR(j,1,tot){
            if((ll)pri[i]*pri[i]>w[j]) break;
            g1[j]=sub(g1[j],mul(pri[i],sub(g1[id(w[j]/pri[i])],s1[i-1])));
            g2[j]=sub(g2[j],mul(mul(pri[i],pri[i]),sub(g2[id(w[j]/pri[i])],s2[i-1])));
        }
    }
}
int solve(ll nn,int x){    //S(nn,x)
    if(pri[x]>=nn) return 0;
    int ans=sub(sub(g2[id(nn)],s2[x]),sub(g1[id(nn)],s1[x]));    //题目中f为二次项减一次项 
    FOR(i,x+1,pl){    //p
        if((ll)pri[i]*pri[i]>nn) break;
        ll pro=1;
        FOR(j,1,50){    //e
            pro*=pri[i];    //p^e
            if(pro>nn) break;
            ans=add(ans,mul(mul(pro%mod,(pro-1)%mod),add(solve(nn/pro,i),j!=1)));
        }
    }
    return ans;
}
int main(){
    scanf("%lld",&n);
    init();
    calc_g();
    printf("%d\n",add(solve(n,0),1));
}
View Code

 

上一篇:c – 类型惩罚 – 编译器如何决定使用哪种类型?


下一篇:Linux进程描述符中union的用法