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)\) 。(草率)

posted @ 2022-02-19 12:59  juruohjr  阅读(34)  评论(0编辑  收藏  举报