Powerful Number 篩

定義

Powerful Number（以下簡稱 PN）篩類似於杜教篩，或者説是杜教篩的一個擴展，可以拿來求一些積性函數的前綴和。

要求：

存在一個函數 \(g\) 滿足：
- \(g\) 是積性函數。
- \(g\) 易求前綴和。
- 對於質數 \(p\)，\(g(p) = f(p)\)。

假設現在要求積性函數 \(f\) 的前綴和 \(F(n) = \sum_{i=1}^{n} f(i)\)。

Powerful Number

定義：對於正整數 \(n\)，記 \(n\) 的質因數分解為 \(n = \prod_{i=1}^{m} p_{i}^{e_{i}}\)。\(n\) 是 PN 當且僅當 \(\forall 1 \le i \le m, e_{i} > 1\)。

性質 1：所有 PN 都可以表示成 \(a^{2}b^{3}\) 的形式。

證明：若 \(e_i\) 是偶數，則將 \(p_{i}^{e_{i}}\) 合併進 \(a^{2}\) 裏；若 \(e_i\) 為奇數，則先將 \(p_{i}^{3}\) 合併進 \(b^{3}\) 裏，再將 \(p_{i}^{e_{i}-3}\) 合併進 \(a^{2}\) 裏。

性質 2：\(n\) 以內的 PN 至多有 \(O(\sqrt{n})\) 個。

證明：考慮枚舉 \(a\)，再考慮滿足條件的 \(b\) 的個數，有 PN 的個數約等於

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

那麼如何求出 \(n\) 以內所有的 PN 呢？線性篩找出 \(\sqrt{n}\) 內的所有素數，再 DFS 搜索各素數的指數即可。由於 \(n\) 以內的 PN 至多有 \(O(\sqrt{n})\) 個，所以至多搜索 \(O(\sqrt{n})\) 次。

PN 篩

首先，構造出一個易求前綴和的積性函數 \(g\)，且滿足對於素數 \(p\)，\(g(p) = f(p)\)。記 \(G(n) = \sum_{i=1}^{n} g(i)\)。

然後，構造函數 \(h = f / g\)，這裏的 \(/\) 表示狄利克雷卷積除法。根據狄利克雷卷積的性質可以得知 \(h\) 也為積性函數，因此 \(h(1) = 1\)。\(f = g * h\)，這裏 \(*\) 表示狄利克雷卷積。

對於素數 \(p\)，\(f(p) = g(1)h(p) + g(p)h(1) = h(p) + g(p) \implies h(p) = 0\)。根據 \(h(p)=0\) 和 \(h\) 是積性函數可以推出對於非 PN 的數 \(n\) 有 \(h(n) = 0\)，即 \(h\) 僅在 PN 處取有效值。

現在，根據 \(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{ is PN}}}^{n}h(d) G\left(\left\lfloor \frac{n}{d}\right\rfloor\right) \end{aligned} \]

\(O(\sqrt{n})\) 找出所有 PN，計算出所有 \(h\) 的有效值。對於 \(h\) 有效值的計算，只需要計算出所有 \(h(p^c)\) 處的值，就可以根據 \(h\) 為積性函數推出 \(h\) 的所有有效值。現在對於每一個有效值 \(d\)，計算 \(h(d)G\left(\left\lfloor \dfrac{n}{d} \right\rfloor\right)\) 並累加即可得到 \(F(n)\)。

下面考慮計算 \(h(p^c)\)，一共有兩種方法：一種是直接推出 \(h(p^c)\) 僅與 \(p, c\) 有關的計算公式，再根據公式計算 \(h(p^c)\)；另一種是根據 \(f = g * h\) 有 \(f(p^c) = \sum_{i=0}^c g(p^i)h(p^{c-i})\)，移項可得 \(h(p^c) = f(p^c) - \sum_{i=1}^{c}g(p^i)h(p^{c-i})\)，現在就可以枚舉素數 \(p\) 再枚舉指數 \(c\) 求解出所有 \(h(p^c)\)。

過程

構造 \(g\)
構造快速計算 \(G\) 的方法
計算 \(h(p^c)\)
搜索 PN，過程中累加答案
得到結果

對於第 3 步，可以直接根據公式計算，可以使用枚舉法預處理打表，也可以搜索到了再臨時推。

性質

以使用第二種方法計算 \(h(p^c)\) 為例進行分析。可以分為計算 \(h(p^c)\) 和搜索兩部分進行分析。

對於第一部分，根據 \(O(\sqrt{n})\) 內的素數個數為 \(O\left(\dfrac{\sqrt{n}}{\log n}\right)\)，每個素數 \(p\) 的指數 \(c\) 至多為 \(\log n\)，計算 \(h(p^c)\) 需要循環 \((c - 1)\) 次，由此有第一部分的時間複雜度為 \(O\left(\dfrac{\sqrt{n}}{\log n} \cdot \log n \cdot \log n\right) = O(\sqrt{n}\log{n})\)，且這是一個寬鬆的上界。根據題目的不同還可以添加不同的優化，從而降低第一部分的時間複雜度。

對於搜索部分，由於 \(n\) 以內的 PN 至多有 \(O(\sqrt{n})\) 個，所以至多搜索 \(O(\sqrt{n})\) 次。對於每一個 PN，根據計算 \(G\) 的方法不同，時間複雜度也不同。例如，假設計算 \(G\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)\) 的時間複雜度為 \(O(1)\)，則第二部分的複雜度為 \(O(\sqrt{n})\)。

特別地，若藉助杜教篩計算 \(G\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)\)，則第二部分的時間複雜度為杜教篩的時間複雜度，即 \(O(n^{\frac{2}{3}})\)。因為若事先計算一次 \(G(n)\)，並且預先使用線性篩優化和用支持快速隨機訪問的數據結構（如 C++ 中的 std::map 和 std::unordered_map）記錄較大的值，則杜教篩過程中用到的 \(G\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)\) 都是線性篩中記錄的或者 std::map 中記錄的，這一點可以直接用程序驗證。

對於空間複雜度，其瓶頸在於存儲 \(h(p^c)\)。若使用二維數組 \(a\) 記錄，\(a_{i,j}\) 表示 \(h(p_i^j)\) 的值，則空間複雜度為 \(O\left(\dfrac{\sqrt{n}}{\log n} \cdot \log n\right) = O(\sqrt{n})\)。

例題

Luogu P5325【模板】Min_25 篩

題意：給定積性函數 \(f(p^k) = p^k(p^k-1)\)，求 \(\sum_{i=1}^{n} f(i)\)。

易得 \(f(p) = p(p-1) = \operatorname{id}(p)\varphi(p)\)，構造 \(g(n) = \operatorname{id}(n)\varphi(n)\)。

考慮使用杜教篩求 \(G(n)\)，根據 \((\operatorname{id}\cdot \varphi) * \operatorname{id} = \operatorname{id}_2\) 可得 \(G(n)= \sum_{i=1}^{n} i^2 - \sum_{i=2}^{n} d \cdot G\left(\left\lfloor \dfrac{n}{d} \right\rfloor\right)\)。

之後 \(h(p^k)\) 的取值可以枚舉計算，這種方法不再贅述。

此外，此題還可以直接求出 \(h(p^k)\) 僅與 \(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} \]

再根據 \(h(p) = 0\)，通過累加法即可推出 \(h(p^k) = (k-1)(p-1)p^k\)。

參考代碼

#include <bits/stdc++.h>
using namespace std;

const 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 {
const int N = 2e6 + 5;
const 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, p; 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;
  scanf("%lld", &n);
  printf("%d\n", PNS::solve(n));
  return 0;
}

「LOJ #6053」簡單的函數

給定 \(f(n)\)：

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

易得：

\[ f(p) = \begin{cases} p + 1 & p = 2 \\ p - 1 & \text{otherwise} \\ \end{cases} \]

構造 \(g\) 為

\[ g(n) = \begin{cases} 3 \varphi(n) & 2 \mid n \\ \varphi(n) & \text{otherwise} \\ \end{cases} \]

易證 \(g(p) = f(p)\) 且 \(g\) 為積性函數。

下面考慮求 \(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} \]

記 \(S_1(n) = \sum_{i=1}^{n} \varphi(i)\)，\(S_2(n) = \sum_{i=1}^{n} \varphi(2i)\)，則 \(G(n) = S_1(n) + 2S_2\left(\left\lfloor \dfrac{n}{2} \right\rfloor\right)\)。

當 \(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} \]

當 \(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} \]

綜上，有 \(S_2(n) = S_1(n) + S_2\left(\left\lfloor \dfrac{n}{2} \right\rfloor\right)\)。

\(S_1\) 可以用杜教篩求，\(S_2\) 直接按照公式推，這樣 \(G\) 也可以求出來了。

參考代碼

#include <bits/stdc++.h>
using namespace std;

const int MOD = 1e9 + 7;
const 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 {
const int N = 2e6 + 5;
const 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, p; 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;
  scanf("%lld", &n);
  printf("%d\n", PNS::solve(n));
  return 0;
}

習題

參考資料

本页面最近更新：，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用