min_25筛小记

引言

现在有一个积性函数 \(f(n)\) ，需要求出 \(f\) 的前缀和，即给定一个数 \(n\) ，求出：

\[\mathrm{ANS} = \sum_{i=1}^n f(i) \]

显然，当 \(n\) 较小的时候，即可以用数组存下所有 \(f\) 的值时，可以使用线性筛，但是当 \(n\) 达到 \(1\times 10^{10}\) 时，线性筛显然不能胜任，需要使用 \(\texttt{min\_25}\) 筛 。

符号约定

• \(\mathbb{P}\) 表示全体素数，\(\mathbb{P}_i\) 表示第 \(i\) 小的素数，比如 \(\mathbb{P}_1=2, \mathbb{P}_3=5\) 。
• \(\mathrm{lpf}(n) = \min_{\mathbb{P}_i | n}\{\mathbb{P}_i\}\) ，表示 \(n\) 的最小素因子，特别的，\(n=1\) 的时候是 \(0\) 。

基本算法

算法分析（瞎吹）

\(\Rightarrow \rm luogu\) 模版题

这里的 \(f(i)\) 满足 \(f\left(\mathbb{P}_i^k\right) = \left(\mathbb{P_i^k}\right)^2 - \mathbb{P}_i^k\) 。

首先考虑将前缀和分成素数和合数两个部分，即：

\[\mathrm{ANS}=\sum_{\mathbb{P}_i \leq n} f(\mathbb{P}_i)+\sum_{\mathbb{P}_i\leq \sqrt{n}\land \mathbb{P}_i^e \leq n} f(\mathbb{P}_i^e)\left\{\sum_{\mathrm{lpf}(j)> \mathbb{P}_i \land j\leq n} f(j)+[e>1]\right\} \]

注意，因为所有合数 \(a (a\not\in \mathbb{P})\)的最小质因子都是小于 \(\sqrt{a}\) 的，因此枚举最小质因子的时候只需要到 \(\sqrt{n}\) 就行了，这里用线性筛先将 \(\sqrt{n}\) 以内的所有素数筛出来，假设最大的素数是 \(\mathbb{P}_m\) 。

接着考虑对于 \(f(i)\) ，我们需要找到一个完全积性函数对其进行拟合，使其在 \(\mathbb{P}_i^k\) 处的值相等。通常我们可以找到一个低次多项式进行拟合，假设次数为 \(k\) ，次数为 \(k\) 的项的系数为 \(a_k\) ，比如对于这一题，选择函数 \(i^2-i\) 即可，此时 \(k=2, a_0=0,a_1=-1,a_2=1\) 。

设 \(g_k(n,x), h_k(n)\) ，这两个函数的定义为：

\[\begin{aligned} & g_k(n,x) = \sum_{i=1}^n [i\in \mathbb{P}\lor\mathrm{lpf}(i)>\mathbb{P}_x] i^k\\ & h_k(n) = \sum_{i=1}^n \mathbb{P}_i^k \end{aligned} \]

这里用递推求出 \(g_k(n,x)\) 。首先容易得到：

\[g_1(n,0)=\frac{n(n+1)}{2}-1, g_2(n,1)=\frac{n(n+1)(2n+1)}{6}-1 \]

这里减一是为了将 \(1\) 去掉，接着转移，容易想到从 \(g_k(*,x-1)\) 转移到 \(g_k(*,x)\) ，这时候需要减去 \(\mathbb{P}_x\) 作为最小质因子的数。

\[g_k(n,x) = \begin{cases} g_k(n,x-1) & \mathbb{P}_x^2 > n\\ g_k(n,x-1) - \mathbb{P}_x^k \left[g_k\left(\frac{n}{\mathbb{P}_x}, x-1\right)-g_k(\mathbb{P}_{x-1}, x-1)\right] & \mathrm{otherwise} \end{cases} \]

注意，这里 \(g_k(\mathbb{P}_{x-1}, x-1)\) 是为了避免减去质数对 \(g_k(n,x)\) 的贡献，这里还可以用 \(h_k(x-1)\) 代替。

接下来考虑计算最终答案。有一个类似 \(\rm dp\) 的东西，设 \(S(n,x)\) 表示 \(1\) 到 \(n\) 中，最小质因子大于 \(\mathbb{P}_x\) 的所有 \(i\) 的 \(f(i)\) 的和，可以想到类似上面式子的转移：

\[S(n,x) = w(n,x) + \sum_{x < i \land \mathbb{P}_i\leq \sqrt{n}\land \mathbb{P}_i^e \leq n} f(\mathbb{P}_i^e)\left\{S\left(\frac{n}{\mathbb{P}_i^e}, x\right)+[e>1]\right\} \]

特别的，令 \(S(n,x)=0 (\mathbb{P}_x>n\lor n=1)\) 。\(w(n,x)\) 表示 \(n\) 以内比 \(\mathbb{P}_x\) 大的素数的 \(f\) 值的和。在这道题中，有：

\[w(n,x)=g_2(n,m)-g_1(n,m)-h_2(x)+h_1(x) \]

这里减去了小于等于 \(\mathbb{P}_x\) 的素数的贡献，提醒一下，这里的 \(m\) 指的是之前筛出的小于等于 \(\sqrt{n}\) 的最大素数的编号。

显然最终答案是 \(\mathrm{ANS} = S(n,0)\) 。因为玄学原因（我也不知道为什么），这个 \(S(n,0)\) 可以直接用递归计算，而且不用记忆化。

接下来考虑具体实现，关键在 \(g\) 的存储，可以发现形式都是 \(\left\lfloor\frac{n}{a}\right\rfloor\) ，同时，因此可以使用整除分块，同时因为小于 \(\sqrt{n}\) ，大于 \(\sqrt{n}\) 的位置的各自不超过 \(\sqrt{n}\) 个，因此可以用两个数组记一下下标。

整除性质 ：

\[\left\lfloor\frac{\left\lfloor\frac{n}{a}\right\rfloor}{b}\right\rfloor=\left\lfloor\frac{n}{ab}\right\rfloor \]

这个性质保证连续的整除，都可以表示整除另一个数。

同时，因为从 \(g_k(*,x-1)\) 转移到 \(g_k(*,x)\) ，因此可以用滚动数组，类似 \(\rm 01\) 背包的倒序枚举之类的进行优化，反正最后只需要使用 \(g_k(*, m)\) 。

代码

typedef Modint<1000000007> mint;
int m, pri[maxn], pri_cnt, ky_cnt, vis[maxn], id1[maxn], id2[maxn];
LL n, ky[maxn];
mint g[maxn][2], h[maxn][2];

inline int getid(LL x) { return x <= m ? id1[x] : id2[n / x]; };

void init() {
    m = sqrt(n)/*抱歉，这里的m和讲解的m不是一个东西，讲解的m=pri_cnt*/; 
  	mint inv2(500000004), inv6(166666668);
    for (int x = 2; x <= m; x++) {
        if (!vis[x]) { 
            pri[++pri_cnt] = x;
            h[pri_cnt][0] = h[pri_cnt - 1][0] + mint(x);
            h[pri_cnt][1] = h[pri_cnt - 1][1] + mint(x) * mint(x);
        }
        for (int j = 1; j <= pri_cnt && 1ll * x * pri[j] <= m; j++) {
            vis[x * pri[j]] = 1;
            if (x % pri[j] == 0) break;
        }
    }
  	//计算 g_k(*,0)
    for (LL i = 1, j; i <= n; i = j + 1) {
        j = n / (n / i), ky[++ky_cnt] = n / j;
        mint t(ky[ky_cnt]);
        g[ky_cnt][0] = t * (t + 1) * inv2 - 1;
        g[ky_cnt][1] = t * (t + 1) * (mint(2) * t + 1) * inv6 - 1;
        if (n / j <= m) id1[n / j] = ky_cnt; else id2[j] = ky_cnt; 
    }
  	//递推计算 g
    for (int i = 1; i <= pri_cnt; i++) {
        LL lim = 1ll * pri[i] * pri[i];
        mint prii(pri[i]);
        for (int j = 1; j <= ky_cnt && ky[j] >= lim; j++) {
            int kyi = getid(ky[j] / pri[i]);
            g[j][0] = g[j][0] - prii * (g[kyi][0] - h[i - 1][0]);
            g[j][1] = g[j][1] - prii * prii * (g[kyi][1] - h[i - 1][1]); 
        }
    }
}

mint S(LL x, LL y) {
    if (pri[y] >= x) return mint(0);
    LL id = getid(x);
    mint ans = (g[id][1] - g[id][0]) - (h[y][1] - h[y][0]);
    for (int i = y + 1; i <= pri_cnt && 1ll * pri[i] * pri[i] <= x; i++) {
        for (LL c = 1, t = pri[i]; t <= x; c++, t *= pri[i]) {
            mint mt(t);
            ans = ans + mt * (mt - 1) * (S(x / t, i) + mint(c != 1));
        }
    }
    return ans;
}

例题

求区间素数个数

\(\Rightarrow \rm LOJ\) 链接

就是 \(g_0(n,m)\) 。（草率）