洛谷P5325 【模板】Min_25筛

链接

https://www.luogu.org/problemnew/show/P5325

学习资料

$LoliconAutomaton$出了一道$min\_25$筛题，于是我就开始学$min\_25$筛
推荐学习资料：
https://www.cnblogs.com/zhoushuyu/p/9187319.html
https://www.cnblogs.com/cjyyb/p/9185093.html

两篇一起看食用效果更佳

Min_25筛

日本人的脑洞永远是我们无法想象的大
这算法完全不懂是怎么想到的，只知道太强了，%就完事了

适用条件

对于一个积性函数$f(x)$，当$x\in prime$时$f(x)$的表达式是多项式，而且$f(p^c)$可以快速计算，那么这玩意就可以用$min\_25$筛来做

第一步，求自变量为质数的$f(x)$之和

现在对于这个多项式函数$f(x)$，我把它每一项拆出来单独做(最后再加起来)，设$F(x)=x^k$，现在我们要对所有满足$x\in prime$的$x$求$F(x)$的和
注意，这里的$F(x)=x^k$是对所有的$x$成立的，也就是说这个函数并不是原来的那个积性函数了，但是在质数处的取值是和原来的函数一样的，既然现在只求质数处的函数值之和，那么也就不会影响最后结果的正确性
设$P_j$为第$j$大的质数（$P_1=2,P_2=3,P_3=5...$）
设$g(n,j)$为把最小质因子$<=P_j$的那些合数去掉之后，剩下的数的$F(x)$之和
严谨的写一下：$$g(n,j)=\sum _{x=2}^{n}F(x)[(x\in pirme)\vee (minp(x)>P_j)]$$
这个东西可以用来递推，最终求出$g(n,num)$就是答案，其中$num$表示不大于$\sqrt{n}$的质数的总数（因为每个合数都会被$\le \sqrt{n}$的质数筛掉）
为了便于理解，我举个例子，假设$F(p)=p$，那么$f(p)$序列就是$f(2)=2,f(3)=3,f(4)=4,f(5)=5$，那么
$g(5,0)=f(2)+f(3)+f(4)+f(5)$
$g(5,1)=f(3)+f(5)$
$g(5,2)=f(5)$
$g(5,3)=0$
可以发现随着$j$的增大，这就是素数筛的过程，不停的用素数把某些项去掉
现在定义完了，来看看这东西怎么算
显然当$j>sqrt(n)$时，$g(n,j)=g(n,j-1)$，因为$P_j$这个素数没有筛掉任何一项
当$j\le sqrt(n)$时，看看它筛掉了哪些？显然是那些最小质因子等于$P_j$的数字，假设这个数除以$P_j$之后得到的商为$x$，那么$P_j\le x$且$x{P}_j\le n$，可以得到$P_j\le x\le \lfloor \frac{n}{P_j}\rfloor$，再根据$F$的完全积性，得到$g(n,j)=g(n,j-1)-F(P_j)(g(\lfloor \frac{n}{P_j}\rfloor ,j-1)-{\sum }_{i=1}^{j-1}F(P_i))$
可以发现这里的递推种$g(x,y)$的第一个自变量$x$都是$\lfloor \frac{n}{k}\rfloor$的形式，所以很幸运我们的$x$的取值只有根号种
总结一下
$$g(n,j)=\{\begin{array}{llc}g(n,j-1)& & P_j>\sqrt{n}\\ g(n,j-1)-F(P_j)(g(\lfloor \frac{n}{P_j}\rfloor ,j-1)-{\sum }_{i=1}^{j-1}F(P_i))& & P_j\le \sqrt{n}\end{array}$$
这玩意具体怎么算呢？
首先我需要一个$F$函数的前缀和
然后观察发现$j$在递推过程中每次减$1$，那既然如此我像背包那样倒过来循环的话，就能少一维空间，一维数组就可以了
大体步骤就是：
1.初始化$g(\lfloor \frac{n}{i}\rfloor ,0)$
2.筛素数并计算$F$的前缀和
3.递推算$g(lfloor\frac{n}{i}\rfloor ,j)$

第二步，算答案

刚才算了质数的，现在算合数的并且合并出答案
设$S(n,j)$表示最小质因子$\ge P_j$且$\le n$的$x$的函数值之和
即$$S(n,j)=\sum _{i=2}^{n}f(i)[minp(i)\ge P_j]$$
（注意这里变成$f$了，刚才一直在讨论$F$，不要搞混）
然后分质数、合数分开算
质数的贡献就是$g(n,num)-{\sum }_{i=1}^{j-1}f(i)$，这很容易理解，就是全部的减去太小的
合数的贡献就是${\sum }_{k=j}^{p{k}^2\le n}{\sum }_{e=1}^{p_k^{e+1}\le n}S(\lfloor \frac{n}{P_k^e}\rfloor ,k+1)f(P_k^e)+f(p_k^{e+1})$，这部分啥意思呢？既然我只算合数的，那我先枚举最小质因数$P_k$，并枚举其幂次$e$，剩下的部分可能是大小不超过$\lfloor \frac{n}{P_k^e}\rfloor$且约数只由大于第$k$个质数的质数组成，也可能包含这个质数的幂(为了不重复计算只需要考虑一次幂)，$p_k^{e+1}\le n$这个不等式的含义是，表面上先保证$f(p_k^{e+1})$合法，其次我可以写出$p_k^e{p}_{k+1}\le n$，当这个不等式不满足的时候，$S(\lfloor \frac{n}{P_k^e}\rfloor ,k+1)$显然为$0$，也就没必要算下去
递推式写出来如下：
$$S(n,j)=g(n,num)-\sum _{i=1}^{j-1}f(i)+\sum _{k=j}^{p{k}^2\le n}\sum _{e=1}^{p_k^{e+1}\le n}S(\lfloor \frac{n}{P_k^e}\rfloor ,k+1)f(P_k^e)+f(p_k^{e+1})$$
这一部分的实现非常直接，就是个递归，连记忆化都没有

时间复杂度

参考：
https://www.luogu.org/blog/user54214/solution-p5325
貌似大家口径不是很统一，说法有两种，要么$O(n^{1-ϵ})$，要么$O(\frac{n^{\frac{3}{4}}}{lo{g}_n})$
反正我是不会分析…

代码

#include <bits/stdc++.h>
#define maxn 200010
#define mod 1000000007ll
#define cl(x) memset(x,0,sizeof(x))
using namespace std;
typedef long long ll;
struct EasyMath
{
    ll prime[maxn];
    bool mark[maxn];
    ll fastpow(ll a, ll b, ll c)
    {
        ll t(a%c), ans(1ll);
        for(;b;b>>=1,t=t*t%c)if(b&1)ans=ans*t%c;
        return ans;
    }
    void get_prime(ll N)
    {
        ll i, j;
        for(i=2;i<=N;i++)mark[i]=false;
        *prime=0;
        for(i=2;i<=N;i++)
        {
            if(!mark[i])prime[++*prime]=i;
            for(j=1;j<=*prime and i*prime[j]<=N;j++)
            {
                mark[i*prime[j]]=true;
                if(i%prime[j]==0)break;
            }
        }
    }
}em;
struct min_25   //f(p^k)=p^k * (p^k - 1)
{
    ll n, sf1[maxn], sf2[maxn], g1[maxn], g2[maxn], id1[maxn], id2[maxn], lis[maxn], _m;
    vector<ll> index;
    inline ll getid(ll x)
    {
        return x<=_m ? id1[x] : id2[n/x];
    }
    void calcg()
    {
        ll i, j;
        //初始化g
        ll _2=em.fastpow(2,mod-2,mod), _6=em.fastpow(6,mod-2,mod);
        for(i=1;i<=n;i=j+1)
        {
            j=n/(n/i);
            lis[++*lis]=n/i;
            if(n/i<=_m)id1[n/i]=*lis;
            else id2[n/(n/i)]=*lis;
            ll t=(n/i)%mod;
            g1[*lis] = ( t*(t+1)%mod*_2 -1 )%mod;
            g2[*lis] = ( t*(t+1)%mod*(2*t+1)%mod*_6 -1 )%mod;
            //g[j]=\sum_{i=2}^j f(j)
        }
        //求f(p_j)的前缀和
        em.get_prime(_m);
        for(i=1;i<=*em.prime;i++)
        {
            auto p=em.prime[i];
            sf1[i]=(sf1[i-1]+p)%mod;
            sf2[i]=(sf2[i-1]+p*p)%mod;
        }
        //求g
        for(j=1;j<=*em.prime;j++)
        {
            auto p=em.prime[j];
            for(ll k=1;k<=*lis;k++)
            {
                if(p*p>lis[k])break;
                ll t = getid(lis[k]/p);
                (g1[k]-=p*(g1[t]-sf1[j-1]))%=mod;
                (g2[k]-=p*p%mod*(g2[t]-sf2[j-1]))%=mod;
            }
        }
    }
    ll S(ll x, ll y)
    {
        if(em.prime[y]>x)return 0;
        auto t=getid(x);
        auto ans=(g2[t]-g1[t])-(sf2[y-1]-sf1[y-1]);
        ans%=mod;
        for(ll i=y;i<=*em.prime and em.prime[i]*em.prime[i]<=x;i++)
        {
            for(auto p1=em.prime[i], p2=p1*p1;p2<=x;p1=p2, p2*=em.prime[i])
            {
                (ans+=S(x/p1,i+1)*(p1*(p1-1)%mod)%mod+p2%mod*(p2%mod-1)%mod)%=mod;
            }
        }
        return ans;
    }
    ll run(ll N)
    {
        _m=sqrt(N);
        cl(g1), cl(g2);
        cl(sf1), cl(sf2);
        n=N;
        calcg();
        return (S(N,1)%mod+1+mod)%mod;
    }
}m25;
int main()
{
    ll N;
    cin>>N;
    cout<<m25.run(N);
    return 0;
}