【笔记】数论：筛法 2023.8.3

8.1 数论 + 8.3 类欧

以下记 \(F(n)=\sum_{i=1}^n f(i)\) 为数论函数 \(f\) 的前缀和。\(n/m\) 默认表示 \(\left\lfloor\frac{n}{m}\right\rfloor\)。

这边进行一波保存：https://www.cnblogs.com/black-swallow/p/16504685.html 和 https://zhuanlan.zhihu.com/p/521413862?utm_id=0。

不会的东西：min_25 筛，块筛复杂度和对应做法，类欧几里得算法。

狄利克雷卷积

积性函数

  • 数论函数 f :: Int -> a 可以认为是定义域为正整数的函数。
  • 积性函数：满足 \(\forall(a,b)=1\Rightarrow f(ab)=f(a)f(b)\) 的数论函数。每个积性函数都应该有 \(f(1)=1\)。
  • 完全积性函数：满足 \(\forall a,b,f(ab)=f(a)f(b)\) 的数论函数。

线性筛积性函数

积性函数可以在线性时间内进行线性筛得到每个点的点值，从而计算前缀和。设线性筛积性函数 \(f\)。

  1. 默认 \(f(p^c)\) 可以关于 \(c\) 线性求。
  2. 记录 \(low_x\) 表示 \(x\) 的最小质因子 \(p\) 的出现的幂次，\(f(x)=f(x/p^{low_u})f(p^{low_u})\)。

例子：线性筛 \(d,g=d*d\)。

点击查看代码 
int pri[1 << 26], cnt, low[1 << 26];
LL D[1 << 26], G[1 << 26];
void sieve(int N) {
  memset(pri, 1, sizeof pri), pri[0] = pri[1] = 0;
  D[1] = G[1] = low[1] = 1;
  for (int i = 1; i <= N; i++) {
    if (pri[i]) {
      pri[++cnt] = low[i] = i;
      for (LL x = i, j = 1; x <= N; x *= i, j++) {
        D[x] = j + 1;
        for (int k = 0; k <= j; k++) G[x] += (k + 1) * (j - k + 1);
      }
    }
    for (LL j = 1, t; j <= cnt && (t = i * pri[j]) <= N; j++) {
      pri[t] = 0;
      if (i % pri[j] == 0) {
        low[t] = low[i] * pri[j];
        D[t] = D[i / low[i]] * D[low[t]], G[t] = G[i / low[i]] * G[low[t]];
        break;
      }
      low[t] = pri[j];
      D[t] = D[i] * D[pri[j]], G[t] = G[i] * G[pri[j]];
    }
  }
  for (int i = 1; i <= N; i++) D[i] += D[i - 1], G[i] += G[i - 1];
}

狄利克雷卷积

对于两个数论函数 \(f,g\)，我们定义它们的狄利克雷卷积 (*) :: (Int -> Int) -> (Int -> Int) -> (Int -> Int) 为 \((f*g)(n)=\sum_{d|n}f(d)g(n/d)\)。（注意这里没有需要积性）

  • 交换律 \(f*g=g*f\)。
  • 结合律 \(f*(g*h)=(f*g)*h\)。
  • 对加法分配律 \((f+g)*h=f*h+g*h\)。
  • 对积性函数封闭（积性函数怎么卷都是积性函数）。
  • 对点乘分配律 \((f*g)\cdot (h*g)=(f*h)\cdot g\)，注意 \(g\) 必须是完全积性函数。
  • 一个由 狄利克雷卷积 构造出的数论函数 \(h=f*g\)，如果 \(f,g\) 可以求得前 \(n\) 项的值，则 \(h\) 的前 \(n\) 项可以在 \(O(n\log n)\) 的时间内计算，这是显然的。

常见积性函数

以下是一些非常重要的完全积性函数：

  • 元函数 \(\epsilon(n)=[n=1]\) 是单位元，有 \(f*\epsilon=f\)。
  • 恒等函数 \(I(n)=1\)。
  • 单位函数 \(id(n)=n\)。
  • 幂函数 \(id_k(n)=n^k\)。

积性函数：

  • 莫比乌斯函数 \(\mu\) 的定义是：有平方因子的数的 \(\mu\) 为 \(0\)，否则为 \((-1)^c\) 其中 \(c\) 是质因子个数。\(\mu*I=\epsilon\)，意思是 \(\mu\) 是 \(I\) 的逆元。
  • 欧拉函数 \(\varphi=\sum_{1\leq i\leq n}[(i,n)=1]\)。\(\varphi*I=id\)，推论是 \(id*\mu=\varphi\)，\(\varphi*\mu=id*\mu^2\)（迫真）。
  • 因数个数函数 \(d=\sigma_0=I*I=\sum_{d|n}1\)。
  • 除数函数 \(\sigma_k=id_k*I=\sum_{d|n}d^k\)。

筛法

整除分块

  1. \(n/i\) 取值只有 \(O(\sqrt{n})\) 种，形如：\(1,2,3,\cdots,\sqrt n,n/\sqrt n,n/(\sqrt n-1),\cdots,n\)。
  2. 封闭性：令 \(S=\{n/x|x\in\mathbb Z\}\)，那么 \(S\) 中的任意一个数 \(x\) 除 \(n\) 的结果 \(n/x\) 仍在 \(S\) 中，且：
    • \(x\leq \sqrt n\Rightarrow n/x>\sqrt n\)。
    • \(x>\sqrt n\Rightarrow n/x\leq \sqrt n\)。
  3. \((a/b)/c=a/(bc)\)。
  4. 满足 \(n/i=n/j\) 的 \(j\) 的最大值为 \(n/(n/i)\)。

于是在计算形如 \(\sum\limits_{i=1}^n f(i)g(n/i)\) 的式子，可以对每一个相同的 \(n/i\) 合并计算，将问题转化为对 \(F,g\) 块筛。

块筛

对 \(f\) 块筛定义为对固定的 \(n\) 的求出数论函数 \(f\) 的前缀和 \(F\) 的所有 \(n/x(x\in\mathbb Z)\) 处的取值，以下的所有筛法都是在想方设法的在亚线性时间复杂度内解决这个问题，因为线性做法是非常明显的。

  • 数论函数 \(f,g\) 可块筛，\(f*g\) 可线性筛，则 \(f*g\) 可在 \(O(\sqrt n)\) 内单点求前缀和，\(O(n^\frac 2 3)\) 块筛。
  • 数论函数 \(f,g\) 可块筛，\(f/g\) 可线性筛，则 \(f/g\) 可在 \(O(n^\frac 2 3)\) 内块筛。
  • 数论函数 \(f,g\)，\(g\) 可块筛，且 \(g(p)=f(p)\)，则 \(f\) 的前缀和可在 \(O(\sqrt n)\) 内单点求值。

杜教筛前置

若 \(f*g=h\)，则 \(H(n)=\sum_{i=1}^nf(i)G(n/i)\)，即若 \(f,g\) 已完成块筛，则 \(H\) 可以在 \(O(\sqrt n)\) 的时间内完成单点求值。以下证明：

\[H(n)=\sum_{i=1}^n\sum_{d|i}g(d)f(i/d)=\sum_{d=1}^n\sum_{d|i}^{n}g(d)f(i/d)\\ =\sum_{d=1}^ng(d)\sum_{t=1}^{n/d}f(t)=\sum_{d=1}^ng(d)F(n/d). \]

杜教筛

杜教筛可以在低于线性的时间（\(O(n^{\frac 2 3})\)）内块筛数论函数 \(f\)，这需要我们找到另外两个数论函数 \(g,h\) 满足 \(f*g=h\)，并且 \(g,h\) 可以块筛。以下是做法：

\[g(1)F(n)=\sum_{i=1}^n g(i)F(n/i)-\sum_{i=2}^n g(i)F(n/i)\\ =H(n)-\sum_{i=2}^ng(i)F(n/i). \]

时间复杂度

通过一些微积分技巧得到复杂度为（这里 \(T(n/i)\) 被放缩为 \(\sqrt{n/i}\)）。

\[T(n)=O(\sqrt{n})+O\left(\sum_{i=1}^\sqrt{n}T(n/i)\right)\\ =O(\sqrt{n})+O\left(\sum_{i=1}^\sqrt{n}\sqrt{\frac{n}{i}}\right)\\ =O\left(\int_0^\sqrt{n}\sqrt{\frac{n}{x}}\mathrm dx\right)=O(n^{3/4}). \]

如果预处理 \(O(n^{2/3})\) 的函数 \(F\) 点值（注意一定要预处理）那么复杂度为 \(O(n^{2/3})\)。

Powerful Number 筛

Powerful Number

定义 Powerful Number 为所有满足所有质因子指数都 \(> 1\) 的数，可以证明这样的数一定能被写成 \(a^2b^3\) 的形式，那么这样的数一共不超过 \(O(\int_0^\sqrt{n}\sqrt[3]{n/x^2}\textrm dx)=O(\sqrt{n})\) 个。

想要求出所有 Powerful number，只需要线性筛 \(O(\sqrt n)\) 以内的质数然后 DFS 枚举质数指数。

Powerful Number 筛

Powerful number 筛可以求出积性函数 \(f\) 的前缀和（单点），需要构造积性函数 \(g\) 使得 \(g(p)=f(p)\)（以下 \(p\) 为质数）且 \(g\) 可以块筛。令 \(h*g=f\)，那么有：

  • \(h\) 是积性函数。\(h(1)=1\)。
  • \(h(p)=0\)，这是因为 \(f(p)=g(1)h(p)+g(p)h(1)=h(p)+g(p)\)，又有 \(g(p)=f(p)\)，所以 \(h(p)=0\)。
  • \(h\) 只可能在 Powerful number 处有值。因为如果 \(n\) 不是 Powerful number，那么 \(h(n=pc)=h(p)h(c)=0\)。

从而有：

\[F(n)=\sum_{i=1}^nh(i)G(n/i)\quad(\text{杜教筛的引理})=\sum_{i=1,i\in PN}^nh(i)G(n/i). \]

所以预处理每个 \(h(p^c)\)，爆搜所有 Powerful number 求值即可。

\(h\) 的求法

因为 \(h\) 是积性函数，在 DFS 的过程中我们可以将若干个 \(h(p^c)\) 乘起来就能计算 Powerful number 的点值，所以现在要算 \(h\) 的质数幂点值。两种做法：

  1. 大力解析 \(h(p^c)\) 使其只和 \(p,c\) 有关。
  2. \(f=g*h\to f(p^c)=\sum_{k=0}^cg(p^k)h(p^{c-k})\to g(1)h(p^c)=f(p^c)-\sum_{k=1}^cg(p^k)h(p^{c-k})\)。

时间复杂度

对于 \(h\) 第二种求法，复杂度为 \(O(\sqrt{n}\log n)\)。

对于爆搜部分，复杂度为 \(O(\sqrt{n})\)。

例题：SPOJ divcnt3

这题可以使用素数拟合（PN 筛），观察到 \(f(p)=4\)，那么构造 \(g=d*d\)，那么 \(g(p)=4\)。考虑块筛 \(g\)？

\(G(n)=\sum_{i=1}^nd(i)D(n/i)\)，所以要块筛 \(d\)；

\(D(n)=\sum_{i=1}^n d(i)=\sum_{i=1}^n\sum_{d|i}1=\sum_{d=1}^n(n/d)\)。

两个都可以整除分块，都可以线性筛，那么做完。

代码实现：

点击查看代码 
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <vector>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define debug(...) void(0)
#endif
typedef long long LL;
int pri[1 << 26], cnt, low[1 << 26];
LL D[1 << 26], G[1 << 26];
void sieve(int N) {
  memset(pri, 1, sizeof pri), pri[0] = pri[1] = 0;
  D[1] = G[1] = low[1] = 1;
  for (int i = 1; i <= N; i++) {
    if (pri[i]) {
      pri[++cnt] = low[i] = i;
      for (LL x = i, j = 1; x <= N; x *= i, j++) {
        D[x] = j + 1;
        for (int k = 0; k <= j; k++) G[x] += (k + 1) * (j - k + 1);
      }
    }
    for (LL j = 1, t; j <= cnt && (t = i * pri[j]) <= N; j++) {
      pri[t] = 0;
      if (i % pri[j] == 0) {
        low[t] = low[i] * pri[j];
        D[t] = D[i / low[i]] * D[low[t]], G[t] = G[i / low[i]] * G[low[t]];
        break;
      }
      low[t] = pri[j];
      D[t] = D[i] * D[pri[j]], G[t] = G[i] * G[pri[j]];
    }
  }
  for (int i = 1; i <= N; i++) D[i] += D[i - 1], G[i] += G[i - 1];
}
LL n;
LL rD[1 << 13], rG[1 << 13];
/*
auto remember=[](LL*F,LL*rF){return [&](LL x)->LL&{
        debug("call %lld\n",x);
        if(x<=6e7) return debug("return F[%lld]=%lld\n",x,F[x]),F[x];
        else return debug("return rF[n/%lld]=%lld\n",x,rF[n/x]),rF[n/x];
};};
auto gD=remember(D,rD),gG=remember(G,rG);
*/
LL& gD(LL x) {
  if (x <= 6e7)
    return D[x];
  else
    return rD[n / x];
}
LL& gG(LL x) {
  if (x <= 6e7)
    return G[x];
  else
    return rG[n / x];
}
LL h[50];
LL dfs(LL x, int now, LL hv) {
  LL res = hv * gG(n / x);
  for (int i = now; i <= cnt; i++) {
    LL p = pri[i];
    if (x > n || x * p > n || x * p * p > n) break;
    for (LL j = 2, t = x * p * p; t <= n; j++, t *= p)
      res += dfs(t, i + 1, hv * h[j]);
  }
  return res;
}
int main() {
  vector<LL> g(60);
  for (int i = 0; i <= 50; i++) {
    h[i] = 3 * i + 1;
    for (int j = 0; j <= i; j++) g[i] += (j + 1) * (i - j + 1);
    for (int j = 1; j <= i; j++) h[i] -= g[j] * h[i - j];
  }
  int T;
  scanf("%d", &T);
  if (T <= 300)
    sieve(6e7);
  else
    sieve(1e4);
  while (T--) {
    scanf("%lld", &n);
    for (LL i = 1; i * i <= (n); i++) {
      LL x = n / i;
      if (x <= 6e7) break;
      gD(x) = 0;
      for (LL l = 1, r; l <= x; l = r + 1)
        r = x / (x / l), gD(x) += (r - l + 1) * (x / l);
    }
    for (LL i = 1; i * i <= (n); i++) {
      LL x = n / i;
      if (x <= 6e7) break;
      gG(x) = 0;
      for (LL l = 1, r; l <= x; l = r + 1)
        r = x / (x / l), gG(x) += (gD(r) - gD(l - 1)) * gD(x / l);
    }
    printf("%lld\n", dfs(1, 1, 1));
  }
  return 0;
}
/*
LL D(LL n) {
  LL res = 0;
  for (LL l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    res += (r - l + 1) * (n / l);
  }
  return res;
}
LL G(LL n) {
  LL res = 0;
  for (LL l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    res += (D(r) - D(l - 1)) * D(n / l);
  }
  return res;
}
*/

后面的内容，有 min_25 筛、类欧几里得算法，都学不动了。标记了。

类欧几里得算法

以下纯口胡。

直线下整点计数

即计算 \(\sum\limits_{i=0}^n\left\lfloor\frac{ai+b}{c}\right\rfloor\)。

see -> https://www.cnblogs.com/caijianhong/p/18290748/template-euclid

无理数有理逼近

\(\sqrt n\) 逼近为 \(\frac p q\)？

\(a_0=\left\lfloor\sqrt n\right\rfloor\)。

\(a_i=\left\lfloor\frac{1}{\sqrt n-a_{i-1}}\right\rfloor\)。

就是写成了连分数的形式 \(\sqrt n=a_0+\dfrac{1}{a_1+\frac{1}{a_2+\cdots}}\)。

作用是可以计算 \(\left\lfloor i\sqrt r\right\rfloor\) 之类的东西了，有理逼近 \(\sqrt r\) 后用类欧几里得那一套。

模意义有理逼近？

抄了 ix35 课件。

给定素数 \(p\leq 10^{18}\) 和 \(a, b, l, r\)，求 \(x \in [l, r]\) 使得 \((ax + b) \bmod p\) 最小。以下全部在模 \(p\) 意义下讨论。

  • 如果 \(a(r − l) < p\)，那么 \((ax+b)\bmod p\) 分为至多两个单调段，暴力即可。
  • 否则，一定存在一个 \(x\) 使得 \(ax + b < a\)，于是答案一定不超过 \(a\)，接着让我们关注那些满足 \(ax + b \in[0, a)\)，将它们按照 \(x\) 从小到大排序，则下一项在 \(\pmod a\) 意义下等于上一项减去 \(p \bmod a\)。于是我们把原来的问题 \((a, p)\) 变为了问题 \((−p\bmod a, a)\)，这里 \(b, l, r\) 需要重新计算一下，剩下的无非就是欧几里得过程递归。
  • 就是一直跳到 \(al+b\) 前面，然后做分讨？不会。

万能欧几里得

see -> https://www.cnblogs.com/caijianhong/p/18290748/template-euclid

posted @ 2023-08-01 20:35  caijianhong  阅读(31)  评论(0)  编辑  收藏  举报