跳转至

Powerful Number 筛

定义

Powerful Number（以下简称 PN）筛类似于杜教筛，或者说是杜教筛的一个扩展，可以拿来求一些积性函数的前缀和。

要求：

  • 存在一个函数 满足：
    • 是积性函数。
    • 易求前缀和。
    • 对于质数 ，。

假设现在要求积性函数 的前缀和 。

Powerful Number

定义：对于正整数 ，记 的质因数分解为 。 是 PN 当且仅当 。

性质 1：所有 PN 都可以表示成 的形式。

证明：若 是偶数，则将 合并进 里；若 为奇数，则先将 合并进 里，再将 合并进 里。

性质 2： 以内的 PN 至多有 个。

证明：考虑枚举 ，再考虑满足条件的 的个数，有 PN 的个数约等于

那么如何求出 以内所有的 PN 呢？线性筛找出 内的所有素数，再 DFS 搜索各素数的指数即可。由于 以内的 PN 至多有 个，所以至多搜索 次。

PN 筛

首先，构造出一个易求前缀和的积性函数 ，且满足对于素数 ，。记 。

然后，构造函数 ，这里的 表示狄利克雷卷积除法。根据狄利克雷卷积的性质可以得知 也为积性函数，因此 。，这里 表示狄利克雷卷积。

对于素数 ，。根据 和 是积性函数可以推出对于非 PN 的数 有 ，即 仅在 PN 处取有效值。

现在，根据 有

找出所有 PN，计算出所有 的有效值。对于 有效值的计算，只需要计算出所有 处的值，就可以根据 为积性函数推出 的所有有效值。现在对于每一个有效值 ，计算 并累加即可得到 。

下面考虑计算 ，一共有两种方法：一种是直接推出 仅与 有关的计算公式，再根据公式计算 ；另一种是根据 有 ，移项可得 ，现在就可以枚举素数 再枚举指数 求解出所有 。

过程

  1. 构造
  2. 构造快速计算 的方法
  3. 计算
  4. 搜索 PN，过程中累加答案
  5. 得到结果

对于第 3 步，可以直接根据公式计算，可以使用枚举法预处理打表，也可以搜索到了再临时推。

性质

以使用第二种方法计算 为例进行分析。可以分为计算 和搜索两部分进行分析。

对于第一部分，根据 内的素数个数为 ，每个素数 的指数 至多为 ，计算 需要循环 次，由此有第一部分的时间复杂度为 ，且这是一个宽松的上界。根据题目的不同还可以添加不同的优化，从而降低第一部分的时间复杂度。

对于搜索部分，由于 以内的 PN 至多有 个，所以至多搜索 次。对于每一个 PN，根据计算 的方法不同，时间复杂度也不同。例如，假设计算 的时间复杂度为 ，则第二部分的复杂度为 。

特别地，若借助杜教筛计算 ，则第二部分的时间复杂度为杜教筛的时间复杂度，即 。因为若事先计算一次 ，并且预先使用线性筛优化和用支持快速随机访问的数据结构（如 C++ 中的 std::map 和 std::unordered_map）记录较大的值，则杜教筛过程中用到的 都是线性筛中记录的或者 std::map 中记录的，这一点可以直接用程序验证。

对于空间复杂度，其瓶颈在于存储 。若使用二维数组 记录， 表示 的值，则空间复杂度为 。

例题

Luogu P5325【模板】Min_25 筛

题意：给定积性函数 ，求 。

易得 ，构造 。

考虑使用杜教筛求 ，根据 可得 。

之后 的取值可以枚举计算，这种方法不再赘述。

此外，此题还可以直接求出 仅与 有关的公式，过程如下：

再根据 ，通过累加法即可推出 。

参考代码 
#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」简单的函数

给定 ：

易得：

构造 为

易证 且 为积性函数。

下面考虑求 。

记 ，，则 。

当 时，有

当 时，有

综上，有 。

可以用杜教筛求， 直接按照公式推，这样 也可以求出来了。

参考代码 
#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;
}

习题

参考资料

最后更新: July 31, 2023
创建日期: May 29, 2021
回到页面顶部