Powerful Number Sieve

Definition

The Powerful Number (hereinafter PN) sieve is similar to the Du sieve, or more precisely an extension of it, and can be used to compute prefix sums of some multiplicative functions.

Requirement:

• There exists a function $g$ satisfying:
  • $g$ is a multiplicative function.
  • $g$ has easily computed prefix sums.
  • For every prime $p$, $g(p) = f(p)$.

Assume that we want the prefix sum of the multiplicative function $f$, $F(n) = \sum_{i=1}^{n} f(i)$.

Powerful Number

Definition: For a positive integer $n$, write the prime factorization of $n$ as $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: Every PN can be represented 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}$, and then merge $p_{i}^{e_{i}-3}$ into $a^{2}$.

Property 2: There are at most $O(\sqrt{n})$ PNs no greater than $n$.

Proof: Enumerate $a$, then count the number of $b$ satisfying the condition. The number of PNs is approximately

So how do we find all PNs no greater than $n$? Use a linear sieve to find all primes no greater than $\sqrt{n}$, then search the exponents of the primes with DFS. Since there are at most $O(\sqrt{n})$ PNs no greater than $n$, we need search at most $O(\sqrt{n})$ times.

PN Sieve

First, construct a multiplicative function $g$ whose prefix sums are easy to compute and which satisfies $g(p) = f(p)$ for every prime $p$. Let $G(n) = \sum_{i=1}^{n} g(i)$.

Next, construct $h = f / g$, where $/$ denotes Dirichlet-convolution division. From the properties of Dirichlet convolution, $h$ is also multiplicative, so $h(1) = 1$. We have $f = g * h$, where $*$ denotes Dirichlet convolution.

For a 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 fact that $h$ is multiplicative, we can conclude that $h(n) = 0$ for every non-PN number $n$; that is, $h$ is nonzero only at PNs.

Now, from $f = g * h$, we have

Find all PNs in $O(\sqrt{n})$ time and compute all effective values of $h$. To compute all effective values of $h$, it is enough to compute the values at all $h(p^c)$; the multiplicativity of $h$ then gives all the effective values of $h$. For every effective value $d$, compute $h(d)G\left(\left\lfloor \dfrac{n}{d} \right\rfloor\right)$ and add them to obtain $F(n)$.

Now consider computing $h(p^c)$. There are two methods: one is to derive directly a formula for $h(p^c)$ that depends only on $p, c$, and then compute $h(p^c)$ from the formula; the other uses $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})$, so we can enumerate primes $p$ and then exponents $c$ to solve for all $h(p^c)$.

Procedure

1. Construct $g$.
2. Construct a method for quickly computing $G$.
3. Compute $h(p^c)$.
4. Search for PNs, accumulating the answer along the way.
5. Obtain the result.

For step 3, we can compute directly from a formula, use enumeration to precompute a table, or derive the value temporarily when it is reached during the search.

Properties

Take the second method for computing $h(p^c)$ as an example. The analysis can be divided into computing $h(p^c)$ and searching.

For the first part, the number of primes no greater than $O(\sqrt{n})$ is $O\left(\dfrac{\sqrt{n}}{\log n}\right)$, and the exponent $c$ of each prime $p$ is at most $\log n$. Computing $h(p^c)$ requires $(c - 1)$ iterations, so 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})$. This is a loose upper bound; depending on the problem, different optimizations can further reduce the time complexity of the first part.

For the search part, there are at most $O(\sqrt{n})$ PNs no greater than $n$, so we search at most $O(\sqrt{n})$ times. For each PN, the time complexity also depends on how $G$ is computed. For example, if computing $G\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)$ takes $O(1)$ time, then the second part has complexity $O(\sqrt{n})$.

In particular, if the Du sieve is used to compute $G\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)$, then the time complexity of the second part is that of the Du sieve, namely $O(n^{\frac{2}{3}})$. If we compute $G(n)$ once in advance, optimize with a linear sieve, and record larger values in a data structure supporting fast random access (such as std::map and std::unordered_map in C++), then every $G\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)$ used by the Du sieve is either recorded by the linear sieve or stored in std::map; this can be verified directly by a program.

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

Examples

Luogu P5325 【Template】 Min_25 Sieve

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

It is easy to obtain $f(p) = p(p-1) = \operatorname{id}(p)\varphi(p)$, so construct $g(n) = \operatorname{id}(n)\varphi(n)$.

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

The values of $h(p^k)$ can then be computed by enumeration; this method is omitted.

In addition, this problem allows us to derive a formula for $h(p^k)$ that depends only on $p, k$. The derivation is as follows:

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

#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;  // prevents negative values
}

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;  // mark the prime 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 composite numbers
      long long nxt = (long long)1 * i * prime[j];
      if (nxt > n) break;
      isp[nxt] = false;
      if (i % prime[j] == 0) {  // the case where i is composite
        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
}

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) {  // compute the 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)$:

It is easy to obtain:

Construct $g$ as

It is easy to prove that $g(p) = f(p)$ and $g$ is multiplicative.

Now consider computing $G(n)$.

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

When $2 \nmid n$, we have

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

$S_1$ can be computed with the Du sieve, while $S_2$ follows directly from the formula, so $G$ can also be computed.

#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;  // prevents values from exceeding the 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;  // prevents negative values
}

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;  // mark the prime 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 composite numbers
      long long nxt = (long long)1 * i * prime[j];
      if (nxt > n) break;
      isp[nxt] = false;
      if (i % prime[j] == 0) {  // the case where i is composite
        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;

        // handle p = 2 as a special case
        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

• Wallbreaker No. 5 - A concise solution to the Powerful Number sieve
• command_block - Du sieve (+ Bell series + Powerful Number)