您的位置: 首页  >  IT文章  >  Min_25筛 学习笔记

Min_25筛 学习笔记

分类: IT文章 • 2023-11-18 11:20:48

Min_25筛
一、用途：

在小于线性的时间内(O(frac{n^{frac 34}}{log n}))计算出积性函数(f(x))的前缀和(sum_{i=1}^nf(i))

对积性函数(f(x))的要求：(f(p^c)(pin prime))能够快速求值，且(f(p))可以多项式表示，具体原因在之后会体现

二、一些定义：
  • (P)表示全体素数集合
  • (P_i)表示第(i)个素数，且定义(P_0 = 1)
  • (operatorname{minp}(x))表示(x)的最小素因子
三、计算方法：

首先我们用一个方法来计算(sum_{p=2}^nf(p))，(pin P)

定义一个完全积性函数(F(x))，(F(x))在(xin P)的情况下值和(f(x))相同

完全积性函数即：对于任意正整数(a,b)满足(F(acdot b)=F(a)cdot F(b))

  • 举个例子：(f(x))为莫比乌斯函数，在质数的情况下(f(p)=-1)，那么对于任意数(xge 2)，(F(x) = (-1)^{omega(x)})，(omega(x))为(x)的质因子数量（包括相同的）

那么我们计算(sum_{p=2}^nf(p))，(pin P)就相当于计算(sum_{p=2}^nF(p))，(pin P)

[定义g(m,j)=sum_{i=2}^m F(i)cdot [iin P or operatorname{minp}(i)>P_j] ]

即(g(m,j))等于(x)在(2)到(m)之间所有满足(x)是质数或者(x)是合数且最小质因子大于第(j)个质数

考虑转移，存在两种情况：

  1. (P_j^2>m)，最小质因子是(P_j)的合数的最小值是(P_j^2)，而(P_j^2>m)，所以显然合数部分不产生贡献，那么(g(m,j) = g(m,j-1))

  2. (P_j^2le m)，这时候最小质因子是(P_j)的合数是存在的，所以(g(m,j))应该在(g(m,j-1))的基础上减掉最小质因子是(P_j)的合数(x)的贡献(F(x))，即：

    [g(m,j)=g(m,j-1)-X ]

    考虑如何计算(X)，这里考虑容斥来做，首先把所有(P_j)的倍数的贡献去掉，但是其中有最小质因子小于(P_j)的贡献，这部分的贡献我们需要加回来先我们把质因子(P_j)提出来，由于(F(x))是完全积性函数，所以可以得到：

    [X=F(P_j)cdot(g(lfloor frac m{P_j} floor,j-1)-g(P_{j-1},j-1)) ]

    (g(lfloor frac m{P_j} floor,j-1))表示所有小于等于(lfloor frac m{P_j} floor)的质数和最小质因子大于等于(P_j)的合数的贡献，其中包含了的一些小于(P_j)质数的贡献，所以要去掉。而(g(P_{j-1},j-1))就表示小于(P_j)的质数带来的贡献

那么我们就可以得到转移方程：

[g(m,j)=egin{cases}g(m,j-1) & P_j^2>m \ g(m,j-1)-F(P_j)cdot(g(lfloor frac m{P_j} floor,j-1)-g(P_{j-1},j-1)) & P_j^2le mend{cases} ]

考虑(P_j = m)的情况，这时候合数部分没有贡献，显然(g(m,j)=sum_{i=1}^jF(P_i))

再考虑(j=0)的情况，这时候所有小于等于(m)的数都有贡献，那么(g(m,0)=sum_{i=1}^nF(i))

更新一下：

[g(m,j)=egin{cases}g(m,j-1) & P_j^2>m \ g(m,j-1)-F(P_j)cdot(g(lfloor frac m{P_j} floor,j-1)-sum_{i=1}^{j-1}F(P_i)) & P_j^2le mend{cases} ]

所以最开始提到了计算(f(p))需要快速求值

我们需要计算的所有素数的(f(p))的和可以这样表示：

[sum_{p=2}^nf(p)=sum_{p=2}^nF(p)=g(n,mid Pmid) ]

其中(mid Pmid)为集合(P)的大小，(pin P)

观察上面的转移方程，可以发现(P_j^2>n)的素数是没有用的，所以需要用到的素数的数量是不超过(sqrt n)的，这部分素数我们可以线性筛筛出来

对于给定的(n)，由于(lfloor frac{lfloor frac ab floor }c floor=lfloor frac a{bcdot c} floor),观察递推式可知需要计算的(g(m,mid Pmid))的(m)可以通过整除分块来预处理，不超过(sqrt n)个

由于需要计算的数字比较大，但是数量不多，考虑离散化一下用下标记录

那么我们可以从最小的素数开始一步步处理计算出各个值，首先把每个值赋值为(g(w[i],0))，然后枚举每个素数，从最大的(w[i])开始递推，具体可以看代码理解一下

比如计算(1-n)所有素数的和和所有素数的数量：

// g1[i]表示 1~w[i]之间的所有素数和
// g2[i]表示 1~w[i]之间所有素数的数量
// pre[i]表示 前i个素数的和
sqt = sqrt(n);
m = 0;
for(int l = 1, r; l <= n; l = r + 1){
    r = n / (n / l);
    w[++m] = n / l;
    if(w[m]<=sqt) id1[w[m]] = m;
    else id2[n/w[m]] = m;
    g1[m] = w[m] - 1;
    g2[m] = (2 + w[m]) * (w[m] - 1) / 2;
}
for(int p = 1; prime[p] * prime[p] <= w[1]; p++) for(int i = 1; i <= m and prime[p] * prime[p] <= w[i]; i++){
    int k = w[i] / prime[p] <= sqt ? id1[w[i]/prime[p]] : id2[n/(w[i]/prime[p])];
    g1[i] -= g1[k] - p + 1;
    g2[i] -= prime[p] * (g2[k] - pre[p-1]);
}

好了我们现在得到了所有素数对答案的贡献，现在定义一个新函数：

[S(m,j)=sum_{i=2}^mf(i)cdot[operatorname{minp}(i)ge P_j] ]

我们需要计算的前缀和就等于(S(n,1)+f(1))

接下来我们考虑如何计算(S(m,j))

枚举每个(i)的最小质因子和次数

[egin{aligned}S(m,k)&=(sum_{kle i \ P_i^2le m}sum_{cge 1\P_i^cle m}f(P_i^c)cdot ([c>1]+S(lfloor frac{m}{P_i^c} floor,i+1))) + sum_{kle i\ P_ile m}f(P_i) \ &= (sum_{kle i \ P_i^2le m}sum_{cge 1\P_i^cle m}f(P_i^c)cdot ([c>1]+S(lfloor frac{m}{P_i^c} floor,i+1)) )+ g(m,mid Pmid) - sum_{i=1}^{k-1}f(P_i) \ &= (sum_{kle i \ P_i^2le m}sum_{cge 1\P_i^{c+1}le m}f(P_i^c)cdot S(lfloor frac{m}{P_i^c} floor,i+1)+f(P_i^{c+1})) + g(m,mid Pmid) - sum_{i=1}^{k-1}f(P_i) end{aligned} ]

由于(f(P_i^c))是可以快速计算的，所以复杂度是有保障的，具体证明就不知道了阿巴阿巴阿巴

比如计算(phi)和(mu)的前缀和：

// 计算欧拉函数前缀和
LL SP(int m, int k){
    if(m<=1 or prime[k]>m) return 0;
    int id = m <= sqt ? id1[m] : id2[n / m];
    LL ret = g2[id] - pre[k-1] + k - 1;
    for(int i = k; prime[i] * prime[i] <= m; i++){
        LL p = prime[i], f = prime[i] - 1;
        for(int c = 1; p * prime[i] <= m; c++, p *= prime[i], f *= prime[i])
            ret += f * SP(m/p,i+1) + f * prime[i];
    }
    return ret;
}
// 计算莫比乌斯函数前缀和
LL SM(int m, int k){
    if(m<=1 or prime[k]>m) return 0;
    int id = m <= sqt ? id1[m] : id2[n / m];
    LL ret = -g1[id] + k - 1;
    for(int i = k; prime[i] * prime[i] <= m; i++) ret -= SM(m/prime[i],i+1);
    return ret;
}


// 质数情况的莫比乌斯函数都是-1，所以前缀和就是质数数量取负
// 而质数情况下的欧拉函数都是p-1，所以前缀和就是质数前缀和减去质数数量
// 所以在计算之前做下面这步操作
// for(int i = 1; i <= m; i++) g2[i] -= g1[i];

相关推荐