min_25 筛

一个亚线性筛，复杂度大概是 $O(\frac{n ^{\frac{3}{4}}}{ \log n})$ 。

使用 $min\_25$ 求前缀和，有两个基本特征：① 积性函数，② 满足质数点为多项式。

算法思路

给定 $n \leq 10 ^ {11}$ ，设 $f(n)$ 为一个积性函数，质数点为多项式，求 $\sum\limits_{i = 1} ^{n} f(i)$ 。

先求出 $\sqrt n$ 内的所有质数，以及 $f(i)$ 。
求出 $\sum\limits_{i = 1} ^{n} [i \in prime]f(i)$ ， $p_i$ 为第 $i$ 个质因子。

我们设 $g(n, j) = \sum\limits_{i = 1} ^{n}[i \in prime\ or\ min\_prime\_of\_i > p_j]f(i)$ ，质数或者最小质因子大于 $p_j$ 的 $f$ 值之和。

考虑从 $g(n, i - 1)$ 推到 $g(n, i)$ ，我们也就是要减去最小质因子等于 $p_i$ 的贡献，
$g(n, i) = \left\{\begin{array}{**lr**} g(n, i - 1)& p_i \times p_i > n\\ g(n, i - 1) - f(p_i) \times \left(g(\frac{n}{p_i}, i - 1) - \sum\limits_{j = 1} ^{i - 1}f(p_j) \right) & p_i \times p_i \leq n\\ \end{array}\right.$
考虑如何递归计算积性函数前缀和 $\sum\limits_{i = 1} ^{n} f(i)$ 。

设 $S(n, j) = \sum\limits_{i = 1} ^{n} [min\_prime\_of\_i \geq p_j] f(i)$ ，最小质因子大于等于 $p_j$ 的 $f$ 值之和，所以这个积性函数的前缀和即为 $S(n, 1) + f(1)$ 。

对于 $S(n, j)$ ，中质数的贡献我们是已知的 $S(n, j) = \sum\limits_{i = 1} ^{n}[i \in prime]f(i) - \sum\limits_{i = 1} ^{j - 1} f(p_i)$ ，接下来考虑计算合数对答案的贡献。

我们枚举合数的最小质因子以及最小质因子的次幂，考虑从 $S(n, j + 1)$ 推到 $S(n, j)$ ，也就是在原来的基础上加上，
$\sum_{i = j} ^{prime_i \times prime_i \leq n} \sum_{k = 1} ^{prime_i ^ {k + 1} \leq n} \left(S(\frac{n}{prime_i ^ k}, j + 1)f(prime_i ^ k) + f(prime_i ^ {k + 1})\right)\\$
所以有
$S(n, j) = g(n) - \sum_{i = 1} ^{j - 1} f(p_i) +\sum_{i = j} ^{prime_i \times prime_i \leq n} \sum_{k = 1} ^{prime_i ^ {k + 1} \leq n} \left(S(\frac{n}{prime_i ^ k}, j + 1)f(prime_i ^ k) + f(prime_i ^ {k + 1})\right)\\$
这里的 $g(n)$ 是最后筛得的，也就是里面存的是质数 $f$ 的前缀和。

这里解释一下那后面一长串的东西，为什么要加上 $f(prime_i ^{k + 1})$ ，而且上界为 $prime_i ^{k + 1} \leq n$ ，

由于我们枚举的是最小质因子及幂次，那么一定有 $n > prime_i ^ k$ ，所以有上界为 $prime_i ^{k + 1} \leq n$ ，

至于为什么要加上 $f(prime_i ^ {k + 1})$ ，由于我们定义的函数中 $f(1)$ 在执行 $min\_25$ 的过程中始终为 $0$ ，

当我们枚举的幂次为 $k + 1$ 时，我们会漏了 $f(prime_i ^{k + 1})$ 这一项，所以得在前一步给他加上去。
优化一点点？如何递推计算积性函数前缀和 $\sum\limits_{i = 1} ^{n} f(i)$ 。

不就是把上面的 $S(n, j)$ ，从 $S(n, j + 1)$ 一步一步往下推吗？

对于 $\frac{n}{prime_j ^ k} \geq prime_j, \frac{n}{prime_j ^ k} < prime_{j + 1} ^ 2$ ，这一步，中我们递归会返回的是 $g(n) - \sum\limits_{i = 1} ^{j}f(p_i)$ ，

也就是大于 $prime_j$ 的质数的部分，按照 $S(n, j)$ 的定义，我们在递推的时候，得到的却是 $0$ 。

所以我们修改一下递推时候的 $S(n, j)$ 的定义：设 $S(n, j) = \sum\limits_{i = 1} ^{n}[i \in prime\ or\ min\_prime\_of\_i \geq p_j]f(i)$ ，

我们还是同样的枚举最小质因子，及其次幂，考虑从 $S(n, j + 1)$ 推到 $S(n, j)$ ，也就是要加上
$\sum_{i = j} ^{prime_i \times prime_i \leq n} \sum_{k = 1} ^{prime_i ^{k + 1} \leq n} \left(\left(S(\frac{n}{prime_i ^ k}, j + 1) - \sum_{k = 1} ^{j} f(p_i) \right)f(prime_i ^ k) + f(prime_i ^{k + 1})\right)\\$

P5325 【模板】Min_25筛

有积性函数 $f(x)$ ，且 $f(p ^ k) = p ^ k(p ^ k - 1), p \in primes$ ，求 $\sum\limits_{i = 1} ^{n} f(i)$ 。

#include <bits/stdc++.h>

using namespace std;

namespace min_25 {
  const int N = 1e6 + 10, mod = 1e9 + 7, inv2 = 500000004, inv6 = 166666668;

  int prime[N], id1[N], id2[N], m, T, cnt;

  long long a[N], sum1[N], g1[N], sum2[N], g2[N], s[N], n;

  bool st[N];

  inline int ID(long long x) {
    return x <= T ? id1[x] : id2[n / x];
  }

  inline long long f(long long x) {
    x %= mod;
    return x * (x - 1) % mod;
  }

  inline long long calc1(long long x) {
    x %= mod;
    return (x * (x + 1) % mod * (2 * x + 1) % mod * inv6 % mod - 1 + mod) % mod;
  }

  inline long long calc2(long long x) {
    x %= mod;
    return (x * (x + 1) % mod * inv2 % mod - 1 + mod) % mod;
  }

  void init() {
    T = sqrt(n + 0.5);
    for(int i = 2; i <= T; i++) {
      if(!st[i]) {
        prime[++cnt] = i;
        sum1[cnt] = (sum1[cnt - 1] + 1ll * i * i) % mod;
        sum2[cnt] = (sum2[cnt - 1] + i) % mod;
      }
      for(int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {
        st[i * prime[j]] = 1;
        if(i % prime[j] == 0) {
            break;
        }
      }
    }
    for (long long l = 1, r; l <= n; l = r + 1) {
      r = n / (n / l);
      a[++m] = n / l;
      a[m] <= T ? id1[a[m]] = m : id2[n / a[m]] = m;
      g1[m] = calc1(a[m]);
      g2[m] = calc2(a[m]);
    }
    for (int j = 1; j <= cnt && 1ll * prime[j] * prime[j] <= n; j++) {
      for (int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {
        g1[i] = ((g1[i] - 1ll * prime[j] * prime[j] % mod * (g1[ID(a[i] / prime[j])] - sum1[j - 1]) % mod) % mod + mod) % mod;
        g2[i] = ((g2[i] - 1ll * prime[j] * (g2[ID(a[i] / prime[j])] - sum2[j - 1]) % mod) % mod + mod) % mod;
      }
    }
    for (int i = 1; i <= m; i++) {
      s[i] = (g1[i] - g2[i] + mod) % mod;
    }
    for (int j = cnt; j >= 1; j--) {
      for (int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {
        for (long long cur = prime[j]; cur * prime[j] <= a[i]; cur *= prime[j]) {
          s[i] = (s[i] + f(cur % mod) * (s[ID(a[i] / cur)] - sum1[j] + sum2[j]) % mod + f(cur * prime[j])) % mod;
          s[i] = (s[i] + mod) % mod;
        }
      }
    }
    printf("%lld\n", s[ID(n)] + 1);
  }

}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  scanf("%lld", &min_25::n);
  min_25::init();
  return 0;
}

#6235. 区间素数个数

直接是 $min\_25$ 的 $g$ 函数求解。

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int N = 1e6 + 10;

int prime[N], id1[N], id2[N], m, cnt, T;

bool st[N];

ll a[N], sum[N], g[N], n;

inline int ID(ll x) {
  return x <= T ? id1[x] : id2[n / x];
}

void init() {
  T = sqrt(n + 0.5);
  for (int i = 2; i <= T; i++) {
    if (!st[i]) {
      prime[++cnt] = i;
      sum[cnt] = sum[cnt - 1] + 1;
    }
    for (int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {
      st[i * prime[j]] = 1;
      if (i % prime[j] == 0) {
        break;
      }
    }
  }
  for (ll l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    a[++m] = n / l;
    if (a[m] <= T) {
      id1[a[m]] = m;
    }
    else {
      id2[n / a[m]] = m;
    }
    g[m] = a[m] - 1;
  }
  for (int j = 1; j <= cnt && 1ll * prime[j] * prime[j] <= n; j++) {
    for (int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {
      g[i] -= g[ID(a[i] / prime[j])] - sum[j - 1];
    }
  }
}

ll solve(ll n) {
  return g[ID(n)];
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  scanf("%lld", &n);
  init();
  printf("%lld\n", solve(n));
  return 0;
}

#572. 「LibreOJ Round #11」Misaka Network 与求和

$f(n)$ 为次大质因子， $f(6) = 3, f(4) = 2, f(1) = 0, f(p) = 1$ 。

$\sum_{i = 1} ^{n} \sum_{j = 1} ^{n} f(\gcd(i, j)) ^ k\\ \sum_{d = 1} ^{n} f(d) ^ k \sum_{i = 1} ^{n} \sum_{j = 1} ^{n}[gcd(i, j) = d]\\ \sum_{d = 1} ^{n} f(d) ^ k \sum_{i = 1} ^{\frac{n}{d}} \sum_{j = 1} ^{\frac{n}{d}}[gcd(i, j) = 1]\\ \sum_{d = 1} ^{n} f(d) ^ k \left( 2 \times\sum_{i = 1} ^{\frac{n}{d}} \sum_{j = 1} ^{i}[gcd(i, j) = 1] - 1 \right)\\ \sum_{d = 1} ^{n} f(d) ^ k \left(2 \sum_{i = 1} ^{\frac{n}{d}}\phi(i) - 1\right)\\$

设 $S(n) = \sum\limits_{i = 1} ^{n} \phi(i), F(n) = f(n) ^ k$ ，而有：
$ $\sum_{d = 1} ^{n} F(d) \left(2 \times S(\frac{n}{d}) - 1 \right)\\$ $

对于 $2 \times S(\frac{n}{d}) - 1$ ，可以通过杜教筛求得，考虑如何求前面的 $\sum\limits_{i = 1} ^{n} F(i)$ ，对前项考虑 $min\_25$ 求解。

我们定义 $S(n, j) = \sum\limits_{i = 1} ^{n} [i \in prime\ or\ min\_prime\ of\ i \geq prime_j] f(i)$ ，也就是质数或者，最小质因子大于等于 $prime_j$ 的答案。

每一轮我们枚举最小质因子及其幂次，考虑这一部分合数的贡献，分两种：

$prime_j$ 是次小质因子。
$prime_j$ 不是次小质因子。

#include <bits/stdc++.h>

using namespace std;

typedef unsigned int uint;

uint quick_pow(uint a, int n) {
  uint ans = 1;
  while (n) {
    if (n & 1) {
      ans *= a;
    }
    a *= a;
    n >>= 1;
  }
  return ans;
}

int n, k;

namespace min_25 {
  const int N = 1e6 + 10;

  int prime[N], a[N], id1[N], id2[N], m, cnt, T;

  uint g[N], sum[N], s[N], f[N];

  bool st[N];

  inline int ID(int x) {
    return x <= T ? id1[x] : id2[n / x];
  }

  void init() {
    T = sqrt(n + 0.5);
    for(int i = 2; i <= T; i++) {
      if(!st[i]) {
        prime[++cnt] = i;
        f[cnt] = quick_pow(i, k);
        sum[cnt] = sum[cnt - 1] + 1;
      }
      for(int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {
        st[i * prime[j]] = 1;
        if(i % prime[j] == 0) {
          break;
        }
      }
    }
    for(int l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        a[++m] = n / l;
        a[m] <= T ? id1[a[m]] = m : id2[n / a[m]] = m;
        g[m] = a[m] - 1;
    }
    for(int j = 1; j <= cnt; j++) {
      for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {
        g[i] -= g[ID(a[i] / prime[j])] - sum[j - 1];
      }
    }
    for (int i = 1; i <= m; i++) {
      s[i] = g[i];
    }
    for (int j = cnt; j >= 1; j--) {
      for (int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {
        for (int cur = prime[j]; 1ll * cur * prime[j] <= a[i]; cur *= prime[j]) {
          s[i] += (s[ID(a[i] / cur)] - g[ID(a[i] / cur)]) + (g[ID(a[i] / cur)] - sum[j - 1]) * f[j];
          /*
          合数答案贡献分为两类：
          一、当前枚举的质数不是次小的，
          二、当前枚举的质数是次小的，也就是形如prime[j] ^ k * b, b是一个质数且这个质数大于 prime[j], 

          最后再加上 f(prime[j] ^ {k + 1})。
          */
        }
      }
    }
  }

  uint solve(int x) {
    return x <= 1 ? 0 : s[ID(x)];
  }
}

namespace Djs {
  const int N = 2e6 + 10;

  int prime[N], cnt;

  uint phi[N];

  bool st[N];

  void init() {
    phi[1] = 1;
    for (int i = 2; i < N; i++) {
      if (!st[i]) {
        prime[++cnt] = i;
        phi[i] = i - 1;
      }
      for (int j = 1; j <= cnt && 1ll * i * prime[j] < N; j++) {
        st[i * prime[j]] = 1;
        if (i % prime[j] == 0) {
          phi[i * prime[j]] = phi[i] * prime[j];
          break;
        }
        phi[i * prime[j]] = phi[i] * (prime[j] - 1);
      }
    }
    for (int i = 1; i < N; i++) {
      phi[i] += phi[i - 1];
    }
  }

  unordered_map<int, uint> mp;

  uint Djs_ans(int n) {
    if (n < N) {
      return phi[n];
    }
    if (mp.count(n)) {
      return mp[n];
    }
    uint ans = 1ll * n * (n + 1) / 2;
    for (uint l = 2, r; l <= n; l = r + 1) {
      r = n / (n / l);
      ans -= (r - l + 1) * Djs_ans(n / l);
    }
    return mp[n] = ans;
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  scanf("%d %d", &n, &k);
  min_25::init();
  Djs::init();
  uint ans = 0;
  for (uint l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    ans += (min_25::solve(r) - min_25::solve(l - 1)) * (2 * Djs::Djs_ans(n / l) - 1);
  }
  cout << ans << "\n";
  return 0;
}

1847 奇怪的数学题

给定 $n, k$ ，要我们求 $\sum\limits_{i = 1} ^{n} \sum\limits_{j = 1} ^{n}sgcd(i, j) ^ k$ ， $sgcd(i, j)$ 是 $gcd(i, j)$ 的次大公约数， $sgcd(i, j) = \frac{\gcd(i, j)}{min\_prime\ of\ \gcd(i, j)}$ ，另 $sgcd(i, j) = f(gcd(i, j))$ 。
$ $\sum_{i = 1} ^{n} \sum_{j = 1} ^{n} f(gcd(i, j)) ^ k\\ \sum_{d = 1} ^{n} f(d) ^ k \sum_{i = 1} ^{n} \sum_{j = 1} ^{n} [gcd(i, j) = d]\\ \sum_{d = 1} ^{n} f(d) ^ k \sum_{i = 1} ^{\frac{n}{d}} \sum_{j = 1} ^{\frac{n}{d}}[gcd(i,j) = 1]\\ \sum_{d = 1} ^{n} f(d) ^ k \left(2 \times \sum_{i = 1} ^{\frac{n}{d}} \sum_{j = 1} ^{i} [gcd(i, j) = 1] - 1\right)\\ \sum_{d = 1} ^{n} f(d) ^ k \left(2 \times \sum_{i = 1} ^{\frac{n}{d}} \phi(i) -1 \right)\\$ $后项可以用杜教筛求解，考虑前项的$ \sum\limits_{d = 1} ^{n} f(d) ^ k $，$ min_25 $中有枚举最小质因子的一步，下面定义$ f(n) = n ^ k$。

还记得 $g(n, j) = g(n, j - 1) - f(p_j) \left(g(\frac{n}{p_j}) - \sum\limits_{i = 1} ^{j - 1} p_j \right)$ ，的转移

后面减去的不就是最小质因子等于 $p_j$ 且为合数的函数值，那么 $g(\frac{n}{p_j} - \sum\limits_{i = 1} ^{j - 1}p_j)$ 就是除去 $f(p_j)$ 的，也就是我们需要的 $sgcd$ 了。

#include <bits/stdc++.h>

using namespace std;

typedef unsigned int uint;

const int N = 1e6 + 10;

uint quick_pow(uint a, int n) {
    uint ans = 1;
    while(n) {
        if(n & 1) ans = ans * a;
        a = a * a;
        n >>= 1;
    }
    return ans;
}

namespace Sum {
  uint S[60][60], k;

  uint Sumk(int n) {
    uint ret = 0;
    for(int i = 1; i <= k; ++i) {
      uint s = 1;
      for(int j = 0; j <= i; ++j)
        if((n + 1 - j) % (i + 1)) s *= n + 1 - j;
        else s *= (n + 1 - j) / (i + 1);
      ret += S[k][i] *s;
    }
    return ret;
  }

  void init() {
    S[0][0] = 1;
    for(int i = 1; i <= k; i++)
      for(int j = 1; j <= i; j++)
        S[i][j] = S[i-1][j-1] + j * S[i-1][j];
  }     
}

namespace Min_25 {
  uint prime[N], id1[N], id2[N], m, cnt, k, T;

  u***1[N], sum1[N], g2[N], sum2[N], f[N], ans[N], n;

  bool st[N];

  int ID(int x) {
      return x <= T ? id1[x] : id2[n / x];
  }

  void init() {
    T = sqrt(n + 0.5);
    for(int i = 2; i <= T; i++) {
      if(!st[i]) {
        prime[++cnt] = i;
        sum1[cnt] = sum1[cnt - 1] + 1;
        f[cnt] = quick_pow(i, k);
        sum2[cnt] = sum2[cnt - 1] + f[cnt];
      }
      for(int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {
        st[i * prime[j]] = 1;
        if(i % prime[j] == 0) {
          break;
        }
      }
    }
    for(int l = 1, r; l <= n; l = r + 1) {
      r = n / (n / l);
      a[++m] = n / l;
      if(a[m] <= T) id1[a[m]] = m;
      else id2[n / a[m]] = m;
      g1[m] = a[m] - 1;
      g2[m] = Sum::Sumk(a[m]) - 1;
    }
    for(int j = 1; j <= cnt; j++) {
      for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {
        g1[i] -= (g1[ID(a[i] / prime[j])] - sum1[j - 1]);
        g2[i] -= f[j] * (g2[ID(a[i] / prime[j])] - sum2[j - 1]);
        ans[i] += g2[ID(a[i] / prime[j])] - sum2[j - 1];
      }
    }
    for(int i = 1; i <= m; i++) {
      ans[i] += g1[i];
    }
  }

  uint solve(int x) {
    if(x <= 1) return 0;
    return ans[ID(x)];
  }
}

namespace Djs {
  uint prime[N], phi[N], cnt;

  bool st[N];

  void init() {
    phi[1] = 1;
    for(int i = 2; i < N; i++) {
      if(!st[i]) {
        prime[++cnt] = i;
        phi[i] = i - 1;
      }
      for(int j = 1; j <= cnt && 1ll * i * prime[j] < N; j++) {
        st[i * prime[j]] = 1;
        if(i % prime[j] == 0) {
          phi[i * prime[j]] = phi[i] * prime[j];
          break;
        }
        phi[i * prime[j]] = phi[i] * (prime[j] - 1);
      }
    }
    for(int i = 1; i < N; i++) {
      phi[i] += phi[i - 1];
    }
  }

  unordered_map<int, uint> ans_s;

  uint S(int n) {
    if(n < N) return phi[n];
    if(ans_s.count(n)) return ans_s[n];
    uint ans = 1ll * n * (n + 1) / 2;
    for(uint l = 2, r; l <= n; l = r + 1) {
      r = n / (n / l);
      ans -=(r - l + 1) * S(n / l);
    }
    return ans_s[n] = ans;
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  int n, k;
  cin >> n >> k;
  Sum::k = k;
  Sum::init();
  Djs::init();
  Min_25::n = n, Min_25::k = k;
  Min_25::init();
  uint ans = 0;
  for(uint l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    ans += (Min_25::solve(r) - Min_25::solve(l - 1)) * (2 * Djs::S(n / l) - 1);
  }
  cout << ans << "\n";
  return 0;
}

这是第二次学min_25了，第一次学的时候只会打打模板，理解的不是很深，希望这次我是真的理解了吧！！！

min_25推导及例题总结

min_25 筛

P5325 【模板】Min_25筛

#6235. 区间素数个数

#572. 「LibreOJ Round #11」Misaka Network 与求和

1847 奇怪的数学题