Min_25筛 学习笔记

Min_25筛

一、用途：

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

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

二、一些定义：

• \(P\)表示全体素数集合
• \(P_i\)表示第\(i\)个素数，且定义\(P_0 = 1\)
• \(\operatorname{minp}(x)\)表示\(x\)的最小素因子

三、计算方法：

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

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

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

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

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

\[定义g(m,j)=\sum_{i=2}^m F(i)\cdot [i\in 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^2\le 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}\rfloor,j-1)-g(P_{j-1},j-1)) \]

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

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

\[g(m,j)=\begin{cases}g(m,j-1) & P_j^2>m \\ g(m,j-1)-F(P_j)\cdot(g(\lfloor \frac m{P_j}\rfloor,j-1)-g(P_{j-1},j-1)) & P_j^2\le m\end{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)=\begin{cases}g(m,j-1) & P_j^2>m \\ g(m,j-1)-F(P_j)\cdot(g(\lfloor \frac m{P_j}\rfloor,j-1)-\sum_{i=1}^{j-1}F(P_i)) & P_j^2\le m\end{cases} \]

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

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

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

其中\(\mid P\mid\)为集合\(P\)的大小，\(p\in P\)

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

对于给定的\(n\)，由于\(\lfloor \frac{\lfloor \frac ab\rfloor }c \rfloor=\lfloor \frac a{b\cdot c}\rfloor\),观察递推式可知需要计算的\(g(m,\mid P\mid)\)的\(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\)的最小质因子和次数

\[\begin{aligned}S(m,k)&=(\sum_{k\le i \\ P_i^2\le m}\sum_{c\ge 1\\P_i^c\le m}f(P_i^c)\cdot ([c>1]+S(\lfloor \frac{m}{P_i^c}\rfloor,i+1))) + \sum_{k\le i\\ P_i\le m}f(P_i) \\ &= (\sum_{k\le i \\ P_i^2\le m}\sum_{c\ge 1\\P_i^c\le m}f(P_i^c)\cdot ([c>1]+S(\lfloor \frac{m}{P_i^c}\rfloor,i+1)) )+ g(m,\mid P\mid) - \sum_{i=1}^{k-1}f(P_i) \\ &= (\sum_{k\le i \\ P_i^2\le m}\sum_{c\ge 1\\P_i^{c+1}\le m}f(P_i^c)\cdot S(\lfloor \frac{m}{P_i^c}\rfloor,i+1)+f(P_i^{c+1})) + g(m,\mid P\mid) - \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];