UOJ Logo analysis的博客

博客

浅谈min25筛

2023-07-26 14:19:06 By analysis

更新日志

2023-05-31

博客第一版

2023-06-23

更新了博客中个别地方的表达。

对“例题导入”部分作了补充。

(想加几道例题,但快生地会考了,不得不先把注意力放在生地上 qwq,只能先咕咕咕了)

2023-07-01

生地考完了,期末考完了,放一天假(今天)然后就回学校竞赛。

根据大佬 RP_INT_MAX 的要求,暂时先加一道例题,随便简单补充了下 g 的存储部分的讲解。

(等到以后轻松一点再考虑要不要再加例题)

添加了目录


目录

0x00 前言

0x01 符号定义

0x02 min25 筛的用处

0x10 第一层:求解 i=1nf(i)

0x11 第二层:求解 i=2n[iP]f(i)

0x12 回到第一层:求解 i=1nf(i)

0x20 例题导入

0x21 补充例题

0x30 结论汇总

0x40 参考资料


0x00 前言

由于本蒟蒻的搜索能力较差,很多地方在学习是没有找到明确的推导过程或定义的讲解,所以只好手推或自行定义,如有错误,欢迎指正。

本篇主要是在讲 min25 筛的各条公式如何无脑的用推式子的方式推出,所以会有很多推式子的内容,我会在结论前标注以供公式恐惧症患者快速查看。

0x01 符号定义

下文中意义不明的符号,都可以在这找到。

符号/函数 含义
P 质数集合
pj j 小的质数,特殊的,p0=0
lpfi i 的最小质因子
P 质数集合的子集,包含 [1,n] 内所有质数
|P| P 集合的大小(即 n 内的质数个数)

0x02 min25 筛的用处

min25 筛被用于处理一类积性函数的前缀和。

也即,对于一个函数 f,可以使用 min25 筛的先决条件有:

一、函数 f 为积性函数。

二、若 p 为质数,则 f(pk) 可以写成一个低阶多项式形式(低阶没有明确的定义,但要求你能在代码中敲出来)。

三、若 p 为质数,则 f(p) 可以快速求解(快速求解并没有明确的定义,但通常是 Θ(1) 的)。

0x10 第一层:求解 i=1nf(i)

如果你认真地阅读了 min25 筛的用处,你会发现使用条件几乎都是在限制质数,这是因为在 min25 筛的处理中需要用到大量的质数的性质。

我们可以将上式按照代入的自变量划分成三个部分:质数,合数与 1。

i=1nf(i)=i=2n[iP]f(i)+i=2n[iP]f(i)+f(1)

合数想与质数挂上钩,最好的方法就是分解质因数。

我们可以枚举一个合数的最小质因子及其次数。

由于 f 是积性函数,所以可以把他所有的最小质因子提出(显然,这个最小质因子不会超过 n)。

=f(1)+i=2n[iP]f(i)+pk,1pkannf(pka)(i=2npka[lpfi>pk]f(i)+[a1]f(1))

第二个求和符号枚举 ka

上式中,为了不出现 a=1i=1 的情况(这样子乘起来就是质数了,但这里只能是合数),我们令 i 从 2 开始枚举,如果 a1 再补上 f(1) 的情况。

我们不妨设 S(n,j)=i=2n[lpfi>pj]f(i),那么原式(我们要求的前缀和)就等于 S(n,0)+f(1)(任何数的最小质因子都大于 0)。

根据上面的式子,我们可以依照相同的思路得到 S 的递推式:

一、分质数与合数(这里定义的 S 是从二开始枚举的,所以不需要 1)。

S(n,j)=i=2n[iP&&i>pj]f(i)+i=2n[iP]f(i)

二、合数分解质因数

S(n,j)=i=2n[iP&&i>pj]f(i)+pk,1pkannf(pka)(S(npka,k)+[a1]f(1))

还可以再展开一下(结论一):

S(n,j)=i=2n[iP]f(i)i=2pj[iP]f(i)+pk,1pkannf(pka)(S(npka,k)+[a1]f(1))

至此,我们就得到了答案计算的递归式了。

接下来就需要处理前面的质数部分了。

0x11 第二层:求解 i=2n[iP]f(i)

显然,1 不是质数。

所以原式可以变成:

i=1n[iP]f(i)

为了解决它,我们引入一个函数:

g(n,j)=i=1n[iP||lpfi>pj]ik

那么,i=1n[iP]ik=g(n,|P|)

由于 i=1n[iP]f(i) 中的 i 全是质数,根据限制条件二,我们可以将 f(i) 展开,然后用 ik 的线性组合表示。


可以举个例子(究竟如何线性组合需要依照具体题目进行分析):

f(i)=i2i(iP)

它由 i2i 线性组合(减法)而成。

于是,我们对 i2i 分别求 g,再相减。


接下来考虑递推 g

显然根据定义有:

g(n,j)=g(n,j1)i=1n[lpfi=pj]ik

如果 ilpf=pj 的最小合数,则 i=pj2

如果 i=pj2>n,显然 g(n,j)=g(n,j1)

否则继续讨论。

ik 是完全积性函数,可以直接提出质因子 pj

g(n,j)=g(n,j1)pjki=1n[lpfipj](ipj)k

=g(n,j1)pjki=1n[lpfi>pj1](ipj)k

更改求和指标,使枚举 i 变成枚举 ipj

=g(n,j1)pjki=1npj[lpfipj>pj1](ipjpj)k

只要 lpfi>pj1lpfipj>pj1 一定成立。

=g(n,j1)pjki=1npj[lpfi>pj1]ik

将后面向 g 的定义转化:

=g(n,j1)pjk(i=1npj[iP||lpfi>pj1]iki=1npj[iP]ik)

显然上式是有错的,对于大于 pj1 的质数,会满足条件 lpfi>pj1,不应重复筛去,更正为:

=g(n,j1)pjk(i=1npj[iP||lpfi>pj1]iki=1pj1[iP]ik)

我们为了方便,再加入一个定义:

spi=i=1pi[iP]ik

他经过线性组合,得到 Spi(下文会结合例题详细讲解)。

于是:

=g(n,j1)pjk(g(npj,j1)spj1)

总结一下(结论二):

g(n,j)={g(n,j1)pj2>ng(n,j1)pjk(g(npj,j1)spj1)pj2n}

通过线性组合 g,我们可以得到 G(下文会结合例题详细讲解)。

0x12 回到第一层:求解 i=1nf(i)

我们此前已经推出的结论:

S(n,j)=i=2n[iP]f(i)i=2pj[iP]f(i)+pk,1pkannf(pka)(S(npka,k)+[a1]f(1))

带入上面求出的结果(最终结论):

S(n,j)=G(n,|P|)Spj+pk,1pkannf(pka)(S(npka,k)+[a1]f(1))

自此,我们已经完成了大部分的理论推导,接下来我会结合例题进行讲解。

0x20 例题导入

例题解析

请先读题:传送门

首先,他肯定适用 min25 筛。

我们可以简单的改动一下题中的式子(展开为多项式):

f(pk)=pk(pk1)=(pk)2pk

显然他可以由 ik 线性组合得到(设 Fi(x)=xi):

f(pk)=F2(pk)F1(pk)

于是,我们可以分别求出 F2F1 的结果,相减(即这道题的线性组合策略)得到 f 的结果。

分析完题目,我们开始做:

第一步是要处理出 Spi

Spispi 线性组合得到。

根据定义:

spi=i=1pi[iP]Fk(i)

于是:

Spi=i=1pi[iP]f(i)=i=1pi[iP]F2(i)+F1(i)

能这样拆是因为 i 必定为质数。

容易发现,这里 Spipi 必定在 n 范围内。

我们可以使用线性筛处理:

int pri[N],cnt;
int vis[N];
int sp1[N],sp2[N];//sp1存F_1的sp,sp2存F_2的sp
void init(int maxn)
{
    for(int i=2;i<=maxn;i++)
    {
        if(!vis[i])
        {
            pri[++cnt] = i;
            sp1[cnt] = (sp1[cnt-1] + i)%mod;
            sp2[cnt] = (sp2[cnt-1] + i * i)%mod;
        }
        for(int j=1;j<=cnt;j++)
        {
            if(i * pri[j] > maxn)break;
            vis[i * pri[j]] = 1;
            if(i % pri[j] == 0)break;
        }
    }
}

接下来处理 G

由于我们使用 G,必然是这个形式的:G(x,|P|)

所以第二个参数可以滚动掉。

但是由于 n1010,不可以直接存起来。

我们发现,x 只有可能是这样的 na

只有 n 个。

但是如果直接存,下标还是到 n

所以考虑这样的储存策略:

使用离散化的方法,使用序号 tot 表示一个数。

使用 ind1ind2 建立数字到序号的映射,使用 w 建立序号到数字的映射。

当数字 w[tot]=nx 小于(等于) n 时,x 必然大于(等于) n

所以令 ind1[w[tot]]=tot

因为存在 ind1 中的 w[tot]n,所以 ind1 下标小于 n

反之则 x 必然小于 n,故令 ind2[x]=tot

因为存在 ind2 中的 xn,所以 ind2 下标小于 n

代码如下:

其中 f1 为计算 g1(x,0)=i=1x[iP||lpfi>0]f1(i)=i=1xi=x(x+1)2

f2 计算 g2(x,0)=i=1xi2=x(x+1)(2x+1)6

int f1(int x)
{
    x%=mod;
    return x * (x + 1)/2 % mod;
}
int f2(int x)
{
    x%=mod;
    return x * (x+1) % mod * ((2 * x  + 1)% mod) % mod * inv6 % mod;
}
int w[N],ind1[N],ind2[N],tot,g1[N],g2[N];
for(int l = 1,r;l<=n;l = r + 1) //枚举x可能出现的值(整除分块)
{
    r = n / (n / l);
    //w数组的含义:id(第几个可能的x值)转实际值
    w[++tot] = n/l; //表示第tot个值是n/l
    //计算g(x,0)
    //实际处理中,我们一般不含f(1)(这里f(1)=1),只要最后加回即可
    g1[tot] = f1(w[tot])-1;
    g2[tot] = f2(w[tot])-1;
    //ind1[x]为x的id
    //ind2[x]为n/x的id
    if(w[tot] <= sqr) // 当x<=sqrt(n)时,存于ind1,否则存于ind2
        ind1[w[tot]] = tot;
    else
        ind2[n/w[tot]] = tot;
}

然后,根据公式递推 g

int getid(int x) // x to id
{
    if(x <= sqr)
        return ind1[x];
    else
        return ind2[n/x];
}
for(int i=1;i<=cnt;i++)
{
    for(int j=1;j<=tot && pri[i] * pri[i] <= w[j];j++)
    {
        g1[j] =(g1[j] - pri[i] * (g1[getid(w[j]/pri[i])] - sp1[i-1]) % mod + mod)%mod;
        g2[j] =(g2[j] - pri[i] * pri[i] % mod * (g2[getid(w[j]/pri[i])] - sp2[i-1]) % mod+mod)%mod;
    }
}

最后计算 S

int s(int x,int k)
{
    if(pri[k] > x)
        return 0;
    int ans = ((g2[getid(x)] - g1[getid(x)] + mod)%mod - (sp2[k] - sp1[k] + mod)%mod + mod)%mod;
    for(int i=k+1;i<=cnt && pri[i] * pri[i] <= x;i++)
    {
        for(int a = 1,p = pri[i];p <= x;p *= pri[i],a++)
        {
            ans =(ans + p % mod * ((p% mod-1) % mod * (s(x/p,i) + (a!=1)) % mod))%mod;
        }
    }
    return ans;
}

汇总代码:

#include<bits/stdc++.h>
#define int long long
const int N = (int)2e5+7,mod = (int)1e9+7;
const int inv2 = 500000004,inv6 = 166666668;
using namespace std;
int pri[N],cnt;
int vis[N];
int sp1[N],sp2[N];
void init(int maxn)
{
    for(int i=2;i<=maxn;i++)
    {
        if(!vis[i])
        {
            pri[++cnt] = i;
            sp1[cnt] = (sp1[cnt-1] + i)%mod;
            sp2[cnt] = (sp2[cnt-1] + i * i)%mod;
        }
        for(int j=1;j<=cnt;j++)
        {
            if(i * pri[j] > maxn)break;
            vis[i * pri[j]] = 1;
            if(i % pri[j] == 0)break;
        }
    }
}
int n,sqr;
int w[N],ind1[N],ind2[N],tot,g1[N],g2[N];
int getid(int x)
{
    if(x <= sqr)
        return ind1[x];
    else
        return ind2[n/x];
}
int f1(int x)
{
    x%=mod;
    return x * (x + 1)/2 % mod;
}
int f2(int x)
{
    x%=mod;
    return x * (x+1) % mod * ((2 * x  + 1)% mod) % mod * inv6 % mod;
}
int s(int x,int k)
{
    if(pri[k] > x)
        return 0;
    int ans = ((g2[getid(x)] - g1[getid(x)] + mod)%mod - (sp2[k] - sp1[k] + mod)%mod + mod)%mod;
    for(int i=k+1;i<=cnt && pri[i] * pri[i] <= x;i++)
    {
        for(int a = 1,p = pri[i];p <= x;p *= pri[i],a++)
        {
            ans =(ans + p % mod * ((p% mod-1) % mod * (s(x/p,i) + (a!=1)) % mod))%mod;
        }
    }
    return ans;
}
signed main()
{
    cin>>n;
    sqr = sqrt(n);
    init(sqr);
    for(int l = 1,r;l<=n;l = r + 1)
    {
        r = n / (n / l);
        w[++tot] = n/l;
        g1[tot] = f1(w[tot])-1;
        g2[tot] = f2(w[tot])-1;
        if(w[tot] <= sqr)
            ind1[w[tot]] = tot;
        else
            ind2[n/w[tot]] = tot;
    }
    for(int i=1;i<=cnt;i++)
    {
        for(int j=1;j<=tot && pri[i] * pri[i] <= w[j];j++)
        {
            g1[j] =(g1[j] - pri[i] * (g1[getid(w[j]/pri[i])] - sp1[i-1]) % mod + mod)%mod;
            g2[j] =(g2[j] - pri[i] * pri[i] % mod * (g2[getid(w[j]/pri[i])] - sp2[i-1]) % mod+mod)%mod;
        }
    }
    //f(1)=1
    cout<<(s(n,0)+1)%mod;
    return 0;
}

于是,你就拿到了 100pts

对本例题的汇总

min25 筛的过程中,需要求解的有:

  • Sp 函数

  • G 函数

  • S 函数

对于 Sp 函数,我们可以利用线性筛预处理出来。

查看整个 min25 筛的过程,使用到的 Sp 不超过 |P|

所以预处理出 Sp1 ~ Spn 的值。

对于 G 函数,我们只需要使用 G(x,|P|),所以第二维可以滚动掉。

对于 xn 个可能值,使用类似离散化的手法。

Wxx 这个 id 的真实值,ind1/2 为值对应的 id

然后根据推导出的公式递推。

对于 S 函数,递归求解即可。

0x21 补充例题

Loj6235 区间素数个数

首先,我们按照需求,构造出需要求得的函数。

ans=i=1n[iP]

如果你对 min25 筛的过程足够熟悉,你就会发现

{g(n,j)=i=1n[iP||lpfi>pj]ikgans=i=1n[iP]}

的相似之处。

不妨变换一下:

ans=i=1n[iP]=i=1n[iP]i0=g(n,|P|),k=0

于是我们很快就打出代码:

#include <bits/stdc++.h>
#define ll long long
using namespace std;
//变量名的意义和之前时一样的
ll n, sqr,w[1000005],g[1000005];
int pri[1000005], cnt;
int vis[1000005];
int ind1[1000005], ind2[1000005], tot;

void init(int maxn) {
    //k值为0,sp[x]即为1~p_x的质数个数
    for (int i = 2; i <= maxn; i++) {
        if (!vis[i]) {
            pri[++cnt] = i;
            //sp[cnt] = cnt;
            //没有必要多用一个数组储存sp
        }

        for (int j = 1; j <= cnt; j++) {
            if (i * pri[j] > maxn)
                break;
            vis[i * pri[j]] = 1;
            if (i % pri[j] == 0)
                break;
        }
    }
}

int getid(ll x) {
    if (x <= sqr)
        return ind1[x];
    else
        return ind2[n / x];
}
signed main() {
    scanf("%lld",&n);
    sqr = sqrt(n);
    init(sqr);
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        w[++tot] = n / l;
        g[tot] = w[tot] - 1;
        if (w[tot] <= sqr) {
            ind1[w[tot]] = tot;
        } else {
            ind2[n / w[tot]] = tot;
        }
    }
    for (int j = 1; j <= cnt; j++) {
        for (int i = 1; i <= tot; i++) {
            //k=0
            if ((ll)pri[j] * pri[j] > w[i])break;
            g[i] -= (g[getid(w[i] / pri[j])] - (j-1));
        }
    }
    printf("%lld",g[getid(n)]);
}

100pts

但请记住,在这里:

for (int j = 1; j <= cnt; j++) {
    for (int i = 1; i <= tot; i++) {
        //k=0
        if ((ll)pri[j] * pri[j] > w[i])break;
        g[i] -= (g[getid(w[i] / pri[j])] - (j-1));
    }
}

请记住别写成这样(60pts TLE 的痛):

for (int j = 1; j <= cnt; j++) {
    for (int i = 1; i <= tot; i++) {
        //k=0
        if ((ll)pri[j] * pri[j] <= w[i])
        g[i] -= (g[getid(w[i] / pri[j])] - (j-1));
    }
}

0x30 结论汇总

定义

S(n,j)=i=2n[lpfi>pj]f(i)

g(n,j)=i=1n[iP||lpfi>pj]ik

spi=i=1pi[iP]ik

公式

S(n,j)=G(n,|P|)Spj+pk,1pkannf(pka)(S(npka,k)+[a1]f(1))

g(n,j)={g(n,j1)pj2>ng(n,j1)pjk(g(npj,j1)spj1)pj2n}

0x40 参考资料

评论

暂无评论

发表评论

可以用@mike来提到mike这个用户,mike会被高亮显示。如果你真的想打“@”这个字符,请用“@@”。