Powerful Number Sieve

Definition

Powerful Number (abbreviated as PN) sieve is similar to the Dirichlet convolution sieve, or rather an extension of it, used to compute the prefix sum of certain multiplicative functions.

Requirements:

• There exists a function \(g\) satisfying:
  • \(g\) is a multiplicative function.
  • The prefix sum of \(g\) is easy to compute.
  • For prime \(p\), \(g(p) = f(p)\).

Suppose we want to compute the prefix sum \(F(n) = \sum_{i=1}^{n} f(i)\) of a multiplicative function \(f\).

Powerful Number

Definition: For a positive integer \(n\), let its prime factorization be \(n = \prod_{i=1}^{m} p_{i}^{e_{i}}\). \(n\) is a PN if and only if \(\forall 1 \le i \le m, e_{i} > 1\).

Property 1: All PNs can be expressed in the form \(a^{2}b^{3}\).

Proof: If \(e_i\) is even, merge \(p_{i}^{e_{i}}\) into \(a^{2}\); if \(e_i\) is odd, first merge \(p_{i}^{3}\) into \(b^{3}\), then merge \(p_{i}^{e_{i}-3}\) into \(a^{2}\).

Property 2: There are at most \(O(\sqrt{n})\) PNs not exceeding \(n\).

Proof: Consider enumerating \(a\), then count the number of \(b\) satisfying the condition. The approximate number of PNs is

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

How to find all PNs not exceeding \(n\)? Use a linear sieve to find all primes up to \(\sqrt{n}\), then perform a DFS to search for exponents of these primes. Since there are at most \(O(\sqrt{n})\) PNs not exceeding \(n\), the search will be performed at most \(O(\sqrt{n})\) times.

PN Sieve

First, construct a multiplicative function \(g\) whose prefix sum is easy to compute, and such that for prime \(p\), \(g(p) = f(p)\). Denote \(G(n) = \sum_{i=1}^{n} g(i)\).

Then, define the function \(h = f / g\), where \(/\) denotes Dirichlet convolution division. From the properties of Dirichlet convolution, \(h\) is also a multiplicative function, so \(h(1) = 1\). We have \(f = g * h\), where \(*\) denotes Dirichlet convolution.

For prime \(p\), \(f(p) = g(1)h(p) + g(p)h(1) = h(p) + g(p) \implies h(p) = 0\). From \(h(p)=0\) and the multiplicativity of \(h\), it follows that \(h(n) = 0\) for non-PN numbers \(n\), meaning \(h\) takes non-zero values only at PNs.

Now, from \(f = g * h\), we have

\[ \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} \]

Find all PNs within \(O(\sqrt{n})\) and compute all non-zero values of \(h\). To compute the non-zero values of \(h\), it suffices to compute \(h(p^c)\) for all prime powers, then use the multiplicativity of \(h\) to derive all non-zero values. For each non-zero value \(d\), compute \(h(d)G\left(\left\lfloor \dfrac{n}{d} \right\rfloor\right)\) and accumulate to obtain \(F(n)\).

Now consider computing \(h(p^c)\). There are two methods: one is to derive a formula for \(h(p^c)\) that depends only on \(p\) and \(c\), then compute using the formula; the other is to use \(f = g * h\), which gives \(f(p^c) = \sum_{i=0}^c g(p^i)h(p^{c-i})\). Rearranging yields \(h(p^c) = f(p^c) - \sum_{i=1}^{c}g(p^i)h(p^{c-i})\). Now we can enumerate primes \(p\) and exponents \(c\) to compute all \(h(p^c)\).

Procedure

1. Construct \(g\)
2. Construct a method to compute \(G\) quickly
3. Compute \(h(p^c)\)
4. Search for PNs, accumulating the answer during the search
5. Obtain the result

For step 3, you can compute directly using a formula, precompute a table via enumeration, or derive on-the-fly during the search.

Complexity Analysis

Using the second method to compute \(h(p^c)\) as an example, the analysis can be divided into computing \(h(p^c)\) and the search part.

For the first part, the number of primes up to \(O(\sqrt{n})\) is \(O\left(\dfrac{\sqrt{n}}{\log n}\right)\). For each prime \(p\), the exponent \(c\) is at most \(\log n\), and computing \(h(p^c)\) requires \((c - 1)\) iterations. Thus, the time complexity of the first part is \(O\left(\dfrac{\sqrt{n}}{\log n} \cdot \log n \cdot \log n\right) = O(\sqrt{n}\log{n})\), which is a loose upper bound. Depending on the problem, additional optimizations can be applied to reduce this complexity.

For the search part, since there are at most \(O(\sqrt{n})\) PNs not exceeding \(n\), the search is performed at most \(O(\sqrt{n})\) times. For each PN, the time complexity depends on the method used to compute \(G\). For example, if computing \(G\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)\) takes \(O(1)\) time, the complexity of the second part is \(O(\sqrt{n})\).

In particular, if the Dirichlet convolution sieve is used to compute \(G\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)\), the time complexity of the second part becomes that of the Dirichlet convolution sieve, i.e., \(O(n^{\frac{2}{3}})\). This is because if \(G(n)\) is computed once beforehand, optimized with a linear sieve, and larger values are stored in a data structure supporting fast random access (like std::map or std::unordered_map in C++), then all \(G\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)\) used during the Dirichlet convolution sieve are either recorded in the linear sieve or in the std::map. This can be verified directly with code.

For space complexity, the bottleneck is storing \(h(p^c)\). If a two-dimensional array \(a\) is used, where \(a_{i,j}\) represents \(h(p_i^j)\), the space complexity is \(O\left(\dfrac{\sqrt{n}}{\log n} \cdot \log n\right) = O(\sqrt{n})\).

Example Problems

Luogu P5325【Template】Min_25 Sieve

Problem: Given a multiplicative function \(f(p^k) = p^k(p^k-1)\), compute \(\sum_{i=1}^{n} f(i)\).

We have \(f(p) = p(p-1) = \operatorname{id}(p)\varphi(p)\). Construct \(g(n) = \operatorname{id}(n)\varphi(n)\).

Consider using the Dirichlet convolution sieve to compute \(G(n)\). From \((\operatorname{id}\cdot \varphi) * \operatorname{id} = \operatorname{id}_2\), we get \(G(n)= \sum_{i=1}^{n} i^2 - \sum_{d=2}^{n} d \cdot G\left(\left\lfloor \dfrac{n}{d} \right\rfloor\right)\).

Then \(h(p^k)\) can be computed by enumeration, which we will not elaborate here.

Alternatively, we can derive a formula for \(h(p^k)\) that depends only on \(p\) and \(k\), as follows:

\[ \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} \]

Using \(h(p) = 0\) and summing, we obtain \(h(p^k) = (k-1)(p-1)p^k\).

Reference Code

#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;  // Prevent negative numbers
}

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;  // Prime number flag array
  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) {  // Sieve out non-prime numbers
      long long nxt = (long long)1 * i * prime[j];
      if (nxt > n) break;
      isp[nxt] = false;
      if (i % prime[j] == 0) {  // Case where i is non-prime
        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]);  // Prefix sum of g function
}

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;  // Pruning

    int c = 2;
    for (long long x = d * prime[i] * prime[i]; x <= global_n;
         x *= prime[i], ++c) {  // Calculate f.g function
      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」Simple Function

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

We have:

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

Construct \(g\) as

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

It can be verified that \(g(p) = f(p)\) and \(g\) is multiplicative.

Now consider computing \(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} \]

Let \(S_1(n) = \sum_{i=1}^{n} \varphi(i)\) and \(S_2(n) = \sum_{i=1}^{n} \varphi(2i)\). Then \(G(n) = S_1(n) + 2S_2\left(\left\lfloor \dfrac{n}{2} \right\rfloor\right)\).

When \(2 \mid n\), we have

\[ \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} \]

When \(2 \nmid n\), we have

\[ \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} \]

In summary, \(S_2(n) = S_1(n) + S_2\left(\left\lfloor \dfrac{n}{2} \right\rfloor\right)\).

\(S_1\) can be computed using the Dirichlet convolution sieve, and \(S_2\) can be computed directly using the recurrence. Thus \(G\) can also be computed.

Reference Code

#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;  // Prevent exceeding modulus
}

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;  // Prevent negative numbers
}

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;  // Prime flag array
  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) {  // Sieve out non-primes
      long long nxt = (long long)1 * i * prime[j];
      if (nxt > n) break;
      isp[nxt] = false;
      if (i % prime[j] == 0) {  // Case when i is not prime
        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;  // Pruning

    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;

        // Special case for 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();  // Preprocessing function
  long long n;
  cin >> n;
  cout << PNS::solve(n) << '\n';
  return 0;
}

Exercises

• PE708 Twos are all you need
• PE639 Summing a multiplicative function
• PE484 Arithmetic Derivative

References

• 破壁人五号 - A brief explanation of Powerful number sieve
• command_block - Dirichlet convolution sieve (+ Bell series + powerful number)

---

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