Min_25 筛

定义

从此种筛法的思想方法来说，其又被称为「Extended Eratosthenes Sieve」。

由于其由 Min_25 发明并最早开始使用，故称「Min_25 筛」。

性质

其可以在或的时间复杂度下解决一类 积性函数 的前缀和问题。

要求：是一个关于可以快速求值的完全积性函数之和（例如多项式）；可以快速求值。

记号

如无特别说明，本节中所有记为的变量的取值集合均为全体质数。
，即为质数时其值为，否则为。
：全体质数中第小的质数（如：）。特别地，令。
，即的最小质因数。特别地，时，其值为。

解释

观察的定义，可以发现答案即为。

考虑如何求出。通过枚举每个的最小质因子及其次数可以得到递推式：

最后一步推导基于这样一个事实：对于满足的，有，故。
其边界值即为。

假设现在已经求出了所有的，那么有两种方式可以求出所有的：

直接按照递推式计算。
从大到小枚举转移，仅当时转移增加值不为零，故按照递推式后缀和优化即可。

现在考虑如何计算。
观察求的过程，容易发现有且仅有这处的点值是有用的。
一般情况下，是一个关于的低次多项式，可以表示为。
那么对于每个，其对的贡献即为。
分开考虑每个的贡献，问题就转变为了：给定，对所有的，求。

Notice：是完全积性函数！

于是设，即埃筛第轮筛完后剩下的数的值之和。
对于一个合数，必定有。设为不大于的最大质数，则，即在埃筛进行轮之后剩下的均为质数。考虑的边界值，显然为。（还记得吗？特别约定了）
对于转移，考虑埃筛的过程，分开讨论每部分的贡献，有：

对于的部分，值不变，即。
对于的部分，被筛掉的数必有质因子，即。
对于第二部分，由于，满足的会被额外减去。这部分应当加回来，即。

则有：

复杂度分析

对于的计算，其第一种方法的时间复杂度被证明为（见 zzt 集训队论文 2.3）；
对于第二种方法，其本质即为洲阁筛的第二部分，在洲阁论文中也有提及（6.5.4），其时间复杂度被证明为。

对于的计算，事实上，其实现与洲阁筛第一部分是相同的。
考虑对于每个，只有在枚举满足的转移时会对时间复杂度产生贡献，则时间复杂度可估计为：

对于空间复杂度，可以发现不论是还是，其均只在处取有效点值，共个，仅记录有效值即可将空间复杂度优化至。

首先，通过一次数论分块可以得到所有的有效值，用一个大小为的数组记录。对于有效值，记为在中的下标，易得：对于所有有效值，。

然后分开考虑小于等于的有效值和大于的有效值：对于小于等于的有效值，用一个数组记录其，即；对于大于的有效值，用一个数组记录，由于过大所以借助记录，即。

这样，就可以使用两个大小为的数组记录所有有效值的并查询。在计算或时，使用有效值的代替有效值作为下标，即可将空间复杂度优化至。

过程

对于的计算，我们实现时一般选择实现难度较低的第一种方法，其在数据规模较小时往往比第二种方法的表现要好；

对于的计算，直接按递推式实现即可。

对于，可以用线性筛预处理出来替代递推式中的。
相应地，递推式中的也可以用此方法预处理。

用 Extended Eratosthenes Sieve 求 积性函数 的前缀和时，应当明确以下几点：

如何快速（一般是线性时间复杂度）筛出前个值；
的多项式表示；
如何快速求出。

明确上述几点之后按顺序实现以下几部分即可：

筛出内的质数与前个值；
对多项式表示中的每一项筛出对应的，合并得到的所有个有用点值；
按照的递推式实现递归，求出。

例题

求莫比乌斯函数的前缀和

求。

易知。则。
直接筛即可得到的所有个所需点值。

求欧拉函数的前缀和

求。

首先易知。
对于的一次项，有；
对于的常数项，有。
筛两次加起来即可得到的所有个所需点值。

「LOJ #6053」简单的函数

给定：

易知。则按照筛的方法筛，对讨论一下即可。
此处给出一种 C++ 实现：

参考代码

/* 「LOJ #6053」简单的函数 */
#include <cmath>
#include <iostream>

constexpr int MAXS = 200000;  // 2sqrt(n)
constexpr int mod = 1000000007;

template <typename x_t, typename y_t>
void inc(x_t &x, const y_t &y) {
  x += y;
  (mod <= x) && (x -= mod);
}

template <typename x_t, typename y_t>
void dec(x_t &x, const y_t &y) {
  x -= y;
  (x < 0) && (x += mod);
}

template <typename x_t, typename y_t>
int sum(const x_t &x, const y_t &y) {
  return x + y < mod ? x + y : (x + y - mod);
}

template <typename x_t, typename y_t>
int sub(const x_t &x, const y_t &y) {
  return x < y ? x - y + mod : (x - y);
}

template <typename _Tp>
int div2(const _Tp &x) {
  return ((x & 1) ? x + mod : x) >> 1;
}

// 以上目的均为防负数和取模
template <typename _Tp>
long long sqrll(const _Tp &x) {  // 平方函数
  return (long long)x * x;
}

int pri[MAXS / 7], lpf[MAXS + 1], spri[MAXS + 1], pcnt;

void sieve(const int &n) {
  for (int i = 2; i <= n; ++i) {
    if (lpf[i] == 0) {  // 记录质数
      lpf[i] = ++pcnt;
      pri[lpf[i]] = i;
      spri[pcnt] = sum(spri[pcnt - 1], i);  // 前缀和
    }
    for (int j = 1, v; j <= lpf[i] && (v = i * pri[j]) <= n; ++j) lpf[v] = j;
  }
}

long long global_n;
int lim;
int le[MAXS + 1],  // x <= \sqrt{n}
    ge[MAXS + 1];  // x > \sqrt{n}
#define idx(v) (v <= lim ? le[v] : ge[global_n / v])

int G[MAXS + 1][2], Fprime[MAXS + 1];
long long lis[MAXS + 1];
int cnt;

void init(const long long &n) {
  for (long long i = 1, j, v; i <= n; i = n / j + 1) {
    j = n / i;
    v = j % mod;
    lis[++cnt] = j;
    (j <= lim ? le[j] : ge[global_n / j]) = cnt;
    G[cnt][0] = sub(v, 1ll);
    G[cnt][1] = div2((long long)(v + 2ll) * (v - 1ll) % mod);
  }
}

void calcFprime() {
  for (int k = 1; k <= pcnt; ++k) {
    const int p = pri[k];
    const long long sqrp = sqrll(p);
    for (int i = 1; lis[i] >= sqrp; ++i) {
      const long long v = lis[i] / p;
      const int id = idx(v);
      dec(G[i][0], sub(G[id][0], k - 1));
      dec(G[i][1], (long long)p * sub(G[id][1], spri[k - 1]) % mod);
    }
  }
  /* F_prime = G_1 - G_0 */
  for (int i = 1; i <= cnt; ++i) Fprime[i] = sub(G[i][1], G[i][0]);
}

int f_p(const int &p, const int &c) {
  /* f(p^{c}) = p xor c */
  return p ^ c;
}

int F(const int &k, const long long &n) {
  if (n < pri[k] || n <= 1) return 0;
  const int id = idx(n);
  long long ans = Fprime[id] - (spri[k - 1] - (k - 1));
  if (k == 1) ans += 2;
  for (int i = k; i <= pcnt && sqrll(pri[i]) <= n; ++i) {
    long long pw = pri[i], pw2 = sqrll(pw);
    for (int c = 1; pw2 <= n; ++c, pw = pw2, pw2 *= pri[i])
      ans +=
          ((long long)f_p(pri[i], c) * F(i + 1, n / pw) + f_p(pri[i], c + 1)) %
          mod;
  }
  return ans % mod;
}

using std::cin;
using std::cout;

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  cin >> global_n;
  lim = sqrt(global_n);  // 上限

  sieve(lim + 1000);  // 预处理
  init(global_n);
  calcFprime();
  cout << (F(1, global_n) + 1ll + mod) % mod << '\n';

  return 0;
}

最后更新: 2024年12月5日
创建日期: 2019年7月24日