min_25 推导及例题总结

min_25 筛

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

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

算法思路

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

• 先求出$\sqrt{n}$内的所有质数，以及$f(i)$。

• 求出$\sum _{i=1}^n[i\in prime]f(i)$，$p_i$为第$i$个质因子。

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

  考虑从$g(n,i-1)$推到$g(n,i)$，我们也就是要减去最小质因子等于$p_i$的贡献，
  $$\begin{gathered} g(n,i)=\{\begin{array}{cc} \\ 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 _{j=1}^{i-1}f(p_j)\right)& p_i\times p_i\le n\end{array} \end{gathered}$$

• 考虑如何递归计算积性函数前缀和$\sum _{i=1}^{n}f(i)$。

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

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

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

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

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

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

  当我们枚举的幂次为$k+1$时，我们会漏了$f(prim{e}_i^{k+1})$这一项，所以得在前一步给他加上去。

• 优化一点点？如何递推计算积性函数前缀和$\sum _{i=1}^{n}f(i)$。

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

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

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

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

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

P5325 【模板】Min_25筛

有积性函数$f(x)$，且$f(p^k)=p^k(p^k-1),p\in primes$，求$\sum _{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$。
$$\begin{gathered} \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}}\varphi (i)-1\right) \end{gathered}$$
设$S(n)=\sum _{i=1}^n\varphi (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 _{i=1}^{n}F(i)$，对前项考虑$min\_25$求解。

我们定义$S(n,j)=\sum _{i=1}^n[i\in primeormin\_primeofi\ge prim{e}_j]f(i)$，也就是质数或者，最小质因子大于等于$prim{e}_j$的答案。

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

• $prim{e}_j$是次小质因子。
• $prim{e}_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 _{i=1}^n\sum _{j=1}^{n}sgcd(i,j{)}^k$，$sgcd(i,j)$是$gcd(i,j)$的次大公约数，$sgcd(i,j)=\frac{\gcd (i,j)}{min\_primeof\gcd (i,j)}$，另$sgcd(i,j)=f(gcd(i,j))$。
$$\begin{gathered} \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}}\varphi (i)-1\right) \end{gathered}$$
后项可以用杜教筛求解，考虑前项的$\sum _{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 _{i=1}^{j-1}{p}_j\right)$，的转移

后面减去的不就是最小质因子等于$p_j$且为合数的函数值，那么$g(\frac{n}{p_j}-\sum _{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;

  uint a[N], g1[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了，第一次学的时候只会打打模板，理解的不是很深，希望这次我是真的理解了吧！！！