Bỏ qua

Powerful Number Sieve

Định nghĩa

Powerful Number (gọi tắt PN) sieve tương tự Du Jiao sieve, hoặc có thể xem là một mở rộng của Du Jiao sieve, dùng để tính prefix sum của một số hàm tích tính.

Yêu cầu:

  • Tồn tại hàm \(g\) thỏa:
    • \(g\) là hàm tích tính.
    • \(g\) dễ tính prefix sum.
    • Với số nguyên tố \(p\), \(g(p) = f(p)\).

Giả sử cần tính prefix sum của hàm tích tính \(f\): \(F(n) = \sum_{i=1}^{n} f(i)\).

Powerful Number

Định nghĩa: Với số nguyên dương \(n\), phân tích \(n = \prod_{i=1}^{m} p_{i}^{e_{i}}\). \(n\) là PN khi và chỉ khi \(\forall 1 \le i \le m, e_{i} > 1\).

Tính chất 1: Mọi PN đều có thể viết dưới dạng \(a^{2}b^{3}\).

Chứng minh: Nếu \(e_i\) chẵn thì gộp \(p_{i}^{e_{i}}\) vào \(a^{2}\); nếu \(e_i\) lẻ thì gộp \(p_{i}^{3}\) vào \(b^{3}\), sau đó gộp \(p_{i}^{e_{i}-3}\) vào \(a^{2}\).

Tính chất 2: Số PN không vượt quá \(n\) nhiều nhất là \(O(\sqrt{n})\).

Chứng minh: Xét liệt kê \(a\), rồi đếm số \(b\) thỏa điều kiện, số PN xấp xỉ

\[ \int_{1}^{\sqrt{n}} \sqrt[3]{\frac{n}{x^2}} \mathrm{d}x = O(\sqrt{n}) \]

Vậy làm sao tìm tất cả PN \(\le n\)? Dùng sàng tuyến tính tìm các số nguyên tố \(\le \sqrt{n}\), rồi DFS liệt kê các số mũ của các nguyên tố. Do số PN \(\le n\) nhiều nhất \(O(\sqrt{n})\), nên số lần tìm kiếm tối đa \(O(\sqrt{n})\).

PN sieve

Đầu tiên, xây dựng hàm tích tính \(g\) dễ tính prefix sum và thỏa \(g(p)=f(p)\) với mọi nguyên tố \(p\). Đặt \(G(n) = \sum_{i=1}^{n} g(i)\).

Sau đó, định nghĩa \(h = f / g\), trong đó \(/\) là phép chia theo chập Dirichlet. Theo tính chất chập Dirichlet, \(h\) cũng là hàm tích tính, nên \(h(1)=1\). Ta có \(f = g * h\) (với \(*\) là chập Dirichlet).

Với nguyên tố \(p\), \(f(p) = g(1)h(p) + g(p)h(1) = h(p) + g(p) \implies h(p) = 0\). Từ \(h(p)=0\) và tính tích tính của \(h\) suy ra với số không phải PN thì \(h(n)=0\), tức \(h\) chỉ khác \(0\) tại PN.

Từ \(f = g * h\):

\[ \begin{aligned} F(n) &= \sum_{i = 1}^{n} f(i)\\ &= \sum_{i = 1}^{n} \sum_{d|i} h(d) g\left(\frac{i}{d}\right)\\ &= \sum_{d=1}^{n} \sum_{i=1}^{\lfloor \frac{n}{d}\rfloor} h(d) g(i)\\ &= \sum_{d=1}^{n} h(d) \sum_{i=1}^{\lfloor \frac{n}{d}\rfloor} g(i) \\ &= \sum_{d=1}^{n} h(d) G\left(\left\lfloor \frac{n}{d}\right\rfloor\right)\\ &= \sum_{\substack{d=1 \\ d \text{ là PN}}}^{n}h(d) G\left(\left\lfloor \frac{n}{d}\right\rfloor\right) \end{aligned} \]

Tìm tất cả PN trong \(O(\sqrt{n})\), tính các giá trị hữu hiệu của \(h\). Để tính \(h\) chỉ cần tính \(h(p^c)\), sau đó dùng tính tích tính để suy ra \(h\) ở các PN khác. Với mỗi giá trị hữu hiệu \(d\), tính \(h(d)G\left(\left\lfloor \dfrac{n}{d} \right\rfloor\right)\) và cộng lại để thu được \(F(n)\).

Xét cách tính \(h(p^c)\), có hai cách: suy ra công thức chỉ phụ thuộc \(p,c\); hoặc dùng \(f = g * h\) có \(f(p^c) = \sum_{i=0}^c g(p^i)h(p^{c-i})\), suy ra \(h(p^c) = f(p^c) - \sum_{i=1}^{c}g(p^i)h(p^{c-i})\), rồi liệt kê \(p,c\) để tính toàn bộ \(h(p^c)\).

Quy trình

  1. Xây dựng \(g\)
  2. Xây dựng cách tính nhanh \(G\)
  3. Tính \(h(p^c)\)
  4. Duyệt PN, trong quá trình cộng dồn đáp án
  5. Lấy kết quả

Bước 3 có thể dùng công thức trực tiếp, hoặc tiền xử lý bằng liệt kê, hoặc tính tạm khi cần.

Phân tích

Lấy cách 2 để phân tích. Chia thành phần tính \(h(p^c)\) và phần tìm kiếm.

Với phần 1, số nguyên tố \(\le \sqrt{n}\) là \(O\left(\dfrac{\sqrt{n}}{\log n}\right)\). Mỗi nguyên tố có số mũ tối đa \(\log n\). Tính \(h(p^c)\) cần lặp \((c-1)\) lần, nên độ phức tạp \(O\left(\dfrac{\sqrt{n}}{\log n} \cdot \log n \cdot \log n\right) = O(\sqrt{n}\log{n})\). Đây là cận trên lỏng; tùy bài có thể tối ưu thêm.

Phần tìm kiếm: do số PN \(\le n\) là \(O(\sqrt{n})\), nên tối đa \(O(\sqrt{n})\) lượt. Với mỗi PN, tùy cách tính \(G\) mà độ phức tạp khác nhau. Ví dụ nếu \(G\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)\) là \(O(1)\) thì phần này là \(O(\sqrt{n})\).

Đặc biệt, nếu dùng Du Jiao sieve để tính \(G\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)\) thì độ phức tạp phần 2 là độ phức tạp Du Jiao sieve, tức \(O(n^{\frac{2}{3}})\). Vì nếu đã tính \(G(n)\) một lần và dùng sàng tuyến tính cộng với cấu trúc truy cập nhanh (như std::map/std::unordered_map) để lưu các giá trị lớn, thì các giá trị \(G\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)\) đều đã có trong sàng hoặc trong map.

Về bộ nhớ, nút thắt là lưu \(h(p^c)\). Nếu dùng mảng 2 chiều \(a\), \(a_{i,j}\) là \(h(p_i^j)\), thì bộ nhớ là \(O\left(\dfrac{\sqrt{n}}{\log n} \cdot \log n\right) = O(\sqrt{n})\).

Ví dụ

Luogu P5325【模板】Min_25 筛

Đề: Cho hàm tích tính \(f(p^k) = p^k(p^k-1)\), tính \(\sum_{i=1}^{n} f(i)\).

Dễ thấy \(f(p) = p(p-1) = \operatorname{id}(p)\varphi(p)\), xây dựng \(g(n) = \operatorname{id}(n)\varphi(n)\).

Dùng Du Jiao sieve để tính \(G(n)\), theo \((\operatorname{id}\cdot \varphi) * \operatorname{id} = \operatorname{id}_2\) suy ra \(G(n)= \sum_{i=1}^{n} i^2 - \sum_{d=2}^{n} d \cdot G\left(\left\lfloor \dfrac{n}{d} \right\rfloor\right)\).

Sau đó \(h(p^k)\) có thể tính bằng liệt kê, phần này không lặp lại.

Ngoài ra có thể trực tiếp suy ra công thức \(h(p^k)\) chỉ phụ thuộc \(p,k\):

\[ \begin{aligned} & f(p^k) = \sum_{i=0}^{k} g(p^{k-i})h(p^i)\\ \iff & p^k(p^k-1) = \sum_{i=0}^{k} p^{k-i}\varphi(p^{k-i}) h(p^i)\\ \iff & p^k(p^k-1) = \sum_{i=0}^{k} p^{2k-2i-1}(p - 1) h(p^i)\\ \iff & p^k(p^k-1) = h(p^k) + \sum_{i=0}^{k-1} p^{2k-2i-1}(p - 1) h(p^i)\\ \iff & h(p^k) = p^k(p^k-1) - \sum_{i=0}^{k-1} p^{2k-2i-1}(p - 1) h(p^i)\\ \iff & h(p^k) - p^2h(p^{k-1}) = p^{k}(p^k-1)-p^{k+1}(p^{k-1}-1) - p(p-1)h(p^{k-1})\\ \iff & h(p^k) - ph(p^{k-1}) = p^{k+1} - p^k\\ \iff & \frac{h(p^k)}{p^k} - \frac{h(p^{k-1})}{p^{k-1}} = p - 1\\ \end{aligned} \]

Từ \(h(p)=0\), cộng dồn suy ra \(h(p^k) = (k-1)(p-1)p^k\).

Mã tham khảo 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
#include <iostream>
#include <map>
using namespace std;

constexpr int MOD = 1e9 + 7;

template <typename T>
int mint(T x) {
  x %= MOD;
  if (x < 0) x += MOD;
  return x;
}

int add(int x, int y) { return x + y >= MOD ? x + y - MOD : x + y; }

int mul(int x, int y) { return (long long)1 * x * y % MOD; }

int sub(int x, int y) {
  return x < y ? x - y + MOD : x - y;  // 防止负数
}

int qp(int x, int y) {
  int r = 1;
  for (; y; y >>= 1) {
    if (y & 1) r = mul(r, x);
    x = mul(x, x);
  }
  return r;
}

int inv(int x) { return qp(x, MOD - 2); }

namespace PNS {
constexpr int N = 2e6 + 5;
constexpr int M = 35;

long long global_n;

int g[N], sg[N];

int h[N][M];
bool vis_h[N][M];

int ans;

int pcnt, prime[N], phi[N];
bool isp[N];

void sieve(int n) {
  pcnt = 0;
  for (int i = 2; i <= n; ++i) isp[i] = true;  // 判断质数数组
  phi[1] = 1;
  for (int i = 2; i <= n; ++i) {
    if (isp[i]) {
      ++pcnt;
      prime[pcnt] = i;
      phi[i] = i - 1;
    }
    for (int j = 1; j <= pcnt; ++j) {  // 筛去非质数
      long long nxt = (long long)1 * i * prime[j];
      if (nxt > n) break;
      isp[nxt] = false;
      if (i % prime[j] == 0) {  // i是非质数的情况
        phi[nxt] = phi[i] * prime[j];
        break;
      }
      phi[nxt] = phi[i] * phi[prime[j]];
    }
  }

  for (int i = 1; i <= n; ++i) g[i] = mul(i, phi[i]);

  sg[0] = 0;
  for (int i = 1; i <= n; ++i) sg[i] = add(sg[i - 1], g[i]);  // g函数的前缀和
}

int inv2, inv6;

void init() {
  sieve(N - 1);
  for (int i = 1; i <= pcnt; ++i) h[i][0] = 1, h[i][1] = 0;
  for (int i = 1; i <= pcnt; ++i) vis_h[i][0] = vis_h[i][1] = true;
  inv2 = inv(2);
  inv6 = inv(6);
}

int S1(long long n) { return mul(mul(mint(n), mint(n + 1)), inv2); }

int S2(long long n) {
  return mul(mul(mint(n), mul(mint(n + 1), mint(n * 2 + 1))), inv6);
}

map<long long, int> mp_g;

int G(long long n) {
  if (n < N) return sg[n];
  if (mp_g.count(n)) return mp_g[n];

  int ret = S2(n);
  for (long long i = 2, j; i <= n; i = j + 1) {
    j = n / (n / i);
    ret = sub(ret, mul(sub(S1(j), S1(i - 1)), G(n / i)));
  }
  mp_g[n] = ret;
  return ret;
}

void dfs(long long d, int hd, int pid) {
  ans = add(ans, mul(hd, G(global_n / d)));

  for (int i = pid; i <= pcnt; ++i) {
    if (i > 1 && d > global_n / prime[i] / prime[i]) break;  // 剪枝

    int c = 2;
    for (long long x = d * prime[i] * prime[i]; x <= global_n;
         x *= prime[i], ++c) {  // 计算f.g函数
      if (!vis_h[i][c]) {
        int f = qp(prime[i], c);
        f = mul(f, sub(f, 1));
        int g = mul(prime[i], prime[i] - 1);
        int t = mul(prime[i], prime[i]);

        for (int j = 1; j <= c; ++j) {
          f = sub(f, mul(g, h[i][c - j]));
          g = mul(g, t);
        }
        h[i][c] = f;
        vis_h[i][c] = true;
      }

      if (h[i][c]) dfs(x, mul(hd, h[i][c]), i + 1);
    }
  }
}

int solve(long long n) {
  global_n = n;
  ans = 0;
  dfs(1, 1, 1);
  return ans;
}
}  // namespace PNS

int main() {
  PNS::init();
  long long n;
  cin >> n;
  cout << PNS::solve(n) << '\n';
  return 0;
}

「LOJ #6053」简单的函数

Cho \(f(n)\):

\[ f(n) = \begin{cases} 1 & n = 1 \\ p \oplus c & n=p^c \\ f(a)f(b) & n=ab \text{ và } a \perp b \end{cases} \]

Dễ thấy:

\[ f(p) = \begin{cases} p + 1 & p = 2 \\ p - 1 & \text{ngược lại} \\ \end{cases} \]

Xây dựng \(g\) là

\[ g(n) = \begin{cases} 3 \varphi(n) & 2 \mid n \\ \varphi(n) & \text{ngược lại} \\ \end{cases} \]

Dễ chứng minh \(g(p) = f(p)\) và \(g\) tích tính.

Xét tính \(G(n)\).

\[ \begin{aligned} G(n) &= \sum_{i=1}^{n}[i \bmod 2 = 1] \varphi(i) + 3 \sum_{i=1}^{n}[i \bmod 2 = 0] \varphi(i)\\ &= \sum_{i=1}^{n} \varphi(i) + 2\sum_{i=1}^{n} [i \bmod 2 = 0]\varphi(i) \\ &= \sum_{i=1}^{n} \varphi(i) + 2\sum_{i=1}^{\lfloor \frac{n}{2} \rfloor} \varphi(2i) \end{aligned} \]

Đặt \(S_1(n) = \sum_{i=1}^{n} \varphi(i)\), \(S_2(n) = \sum_{i=1}^{n} \varphi(2i)\), ta có \(G(n) = S_1(n) + 2S_2\left(\left\lfloor \dfrac{n}{2} \right\rfloor\right)\).

Khi \(2 \mid n\):

\[ \begin{aligned} S_2(n) &= \sum_{i=1}^{n} \varphi(2i) \\ &= \sum_{i=1}^{\frac{n}{2}} (\varphi(2(2i-1)) + \varphi(2(2i))) \\ &= \sum_{i=1}^{\frac{n}{2}} (\varphi(2i-1) + 2\varphi(2i)) \\ &= \sum_{i=1}^{\frac{n}{2}} (\varphi(2i-1) + \varphi(2i)) + \sum_{i=1}^{\frac{n}{2}} \varphi(2i) \\ &= \sum_{i=1}^{n} \varphi(i) + S_2\left(\frac{n}{2}\right)\\ &= S_1(n) + S_2\left(\left\lfloor \frac{n}{2} \right\rfloor\right)\\ \end{aligned} \]

Khi \(2 \nmid n\):

\[ \begin{aligned} S_2(n) &= S_2(n-1) + \varphi(2n) \\ &= S_2(n-1) + \varphi(n) \\ &= \sum_{i=1}^{n-1} \varphi(i) + S_2\left(\frac{n-1}{2}\right) + \varphi(n)\\ &= S_1(n) + S_2\left(\left\lfloor \frac{n}{2} \right\rfloor\right)\\ \end{aligned} \]

Suy ra \(S_2(n) = S_1(n) + S_2\left(\left\lfloor \dfrac{n}{2} \right\rfloor\right)\).

\(S_1\) có thể tính bằng Du Jiao sieve, \(S_2\) tính trực tiếp theo công thức, từ đó có \(G\).

Mã tham khảo 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
#include <iostream>
#include <map>
using namespace std;

constexpr int MOD = 1e9 + 7;
constexpr int inv2 = (MOD + 1) / 2;

template <typename T>
int mint(T x) {
  x %= MOD;
  if (x < 0) x += MOD;
  return x;
}

int add(int x, int y) {
  return x + y >= MOD ? x + y - MOD : x + y;  // 防止大于模数
}

int mul(int x, int y) { return (long long)1 * x * y % MOD; }

int sub(int x, int y) {
  return x < y ? x - y + MOD : x - y;  // 防负数
}

namespace PNS {
constexpr int N = 2e6 + 5;
constexpr int M = 35;

long long global_n;

int s1[N], s2[N];

int h[N][M];
bool vis_h[N][M];

int ans;

int pcnt, prime[N], phi[N];
bool isp[N];

void sieve(int n) {
  pcnt = 0;
  for (int i = 2; i <= n; ++i) isp[i] = true;  // 判断质数数组
  phi[1] = 1;
  for (int i = 2; i <= n; ++i) {
    if (isp[i]) {
      ++pcnt;
      prime[pcnt] = i;
      phi[i] = i - 1;
    }
    for (int j = 1; j <= pcnt; ++j) {  // 筛去非质数
      long long nxt = (long long)1 * i * prime[j];
      if (nxt > n) break;
      isp[nxt] = false;
      if (i % prime[j] == 0) {  // i是非质数的情况
        phi[nxt] = phi[i] * prime[j];
        break;
      }
      phi[nxt] = phi[i] * phi[prime[j]];
    }
  }

  s1[0] = 0;
  for (int i = 1; i <= n; ++i) s1[i] = add(s1[i - 1], phi[i]);

  s2[0] = 0;
  for (int i = 1; i <= n / 2; ++i) {
    s2[i] = add(s2[i - 1], phi[2 * i]);
  }
}

void init() {
  sieve(N - 1);
  for (int i = 1; i <= pcnt; ++i) h[i][0] = 1;
  for (int i = 1; i <= pcnt; ++i) vis_h[i][0] = true;
}

map<long long, int> mp_s1;

int S1(long long n) {
  if (n < N) return s1[n];
  if (mp_s1.count(n)) return mp_s1[n];

  int ret = mul(mul(mint(n), mint(n + 1)), inv2);
  for (long long i = 2, j; i <= n; i = j + 1) {
    j = n / (n / i);
    ret = sub(ret, mul(mint(j - i + 1), S1(n / i)));
  }
  mp_s1[n] = ret;
  return ret;
}

map<long long, int> mp_s2;

int S2(long long n) {
  if (n < N / 2) return s2[n];
  if (mp_s2.count(n)) return mp_s2[n];
  int ret = add(S1(n), S2(n / 2));
  mp_s2[n] = ret;
  return ret;
}

int G(long long n) { return add(S1(n), mul(2, S2(n / 2))); }

void dfs(long long d, int hd, int pid) {
  ans = add(ans, mul(hd, G(global_n / d)));

  for (int i = pid; i <= pcnt; ++i) {
    if (i > 1 && d > global_n / prime[i] / prime[i]) break;  // 剪枝

    int c = 2;
    for (long long x = d * prime[i] * prime[i]; x <= global_n;
         x *= prime[i], ++c) {
      if (!vis_h[i][c]) {
        int f = prime[i] ^ c, g = prime[i] - 1;

        // p = 2时特判一下
        if (i == 1) g = mul(g, 3);

        for (int j = 1; j <= c; ++j) {
          f = sub(f, mul(g, h[i][c - j]));
          g = mul(g, prime[i]);
        }
        h[i][c] = f;
        vis_h[i][c] = true;
      }

      if (h[i][c]) dfs(x, mul(hd, h[i][c]), i + 1);
    }
  }
}

int solve(long long n) {
  global_n = n;
  ans = 0;
  dfs(1, 1, 1);
  return ans;
}
}  // namespace PNS

int main() {
  PNS::init();  // 预处理函数
  long long n;
  cin >> n;
  cout << PNS::solve(n) << '\n';
  return 0;
}

Bài tập

Tài liệu tham khảo

Last updated on this page:, Update history
Found an error? Want to help improve? Edit this page on GitHub!
Contributors to this page:OI-wiki
All content on this page is provided under the terms of the CC BY-SA 4.0 and SATA license, additional terms may apply