自动搬运

查看原文

来自洛谷，原作者为

	moongazer We do not know and will not know

搬运于2025-08-24 22:10:43，当前版本为作者最后更新于2020-01-24 00:35:24，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

一种$O($值域$)$时间预处理$O(1)$时间求最大公约数($\gcd$)的算法

$\mathfrak{View\space it\space on\space my\space Blog}$

Update in 2020.08.07：修正了一处错误，感谢@Kinesis的指正

一些约定

1. $N$为询问的值域
2. $Prime$为全体素数集合
3. 集合$\{a_1,a_2,\cdots,a_m\}$是$n$的分解，当且仅当$a_1\times a_2\times\cdots\times a_m=n$

原理

定理一

内容

可以将值域中的每个$x$分解成$\{a,b,c\}$，满足$a,b,c\leq\sqrt{x}$或$\in Prime$（定义这种分解为合法分解）

证明

不妨设$a\leq b\leq c$若$c\notin Prime$且$c>\sqrt{x}$，则$c$可分解为$\{d,e\}$且$d\leq e$有$d\leq\sqrt{x}$,而$a\times b=\frac{x}{c}<\sqrt{x}$则有$n$的分解$\{d,a\times b,e\}$，若$e>\sqrt{x}$则可按该规律一直分解直到$e\in Prime$或$\leq\sqrt{x}$

定理二

内容

对于询问$\gcd(x,y)$，分别考虑$a,b,c$对答案的贡献，$a$对答案的贡献为$\gcd(a,y)$，再将$y$除以$\gcd(a,y)$（这个因子已经被算过，不能重复计算），再对$b,c$干同样的事，三个贡献相乘即为$\gcd(x,y)$

证明

易得引理若$r\mid p, r\mid q$则$\gcd(p,q)=r\times\gcd(\frac{p}{r},\frac{q}{r})$

分别代入$\left\{ \begin{aligned} &p_1=a\times b\times c,q_1=y,r_1=\gcd(a,q_1) \\ &p_2=b\times c,q_2=\frac{q_1}{r_1},r_2=\gcd(b,q_2) \\ &p_3=c,q_3=\frac{q_2}{r_2},r_3=\gcd(c,q_3) \end{aligned} \right. $即可

实现

我们发现实现的难点在于如何在$O(N)$时间内进行分解，询问部分较为容易

分解

对于$x\geq2$，找到$x$的最小质因子$p$以及$\frac{x}{p}$的合法分解$\{a_0,b_0,c_0\}$且$a_0\leq b_0\leq c_0$,$x$的一种合法分解即为$\{a_0\times p,b_0,c_0\}$的升序排序

考虑证明：

1. $x\in Prime$时显然成立,分解为$\{1,1,x\}$
2. 当$p\le\sqrt[4]{x}$时将$a_0\leq\sqrt[3]{\frac{x}{p}}$带入有$a_0\times p\le\sqrt{x}$
3. 考虑$p>\sqrt[4]{x}$的情况，$\left(1.\right)$ $a_0=1$，显然有$a_0\times p=p\le\sqrt{x}$；$\left(2.\right)$ $a\neq1$，由于$x$不是素数，$\frac{x}{p}$的最小质因子$q$即为$x$的第二小质因子，一定$\geq p$，而$a_0,b_0,c_0$都为$\frac{x}{p}$的非$1$因子，有$p\leq q\leq a_0\leq b_0\leq c_0$,$p\times a_0\times b_0\times c_0>(\sqrt[4]{x})^4=x$与$p\times a_0\times b_0\times c_0=x$相矛盾，故不存在此情况

所以只用跑一次线性筛，用最小质因子更新即可，然后预处理出$\sqrt{n}\times\sqrt{n}$的$\gcd$数组

代码如下

// fac为合法分解，isp表示是否非质数，pri为质数数组，tot为pri的大小，pre为预处理的gcd数组，M为值域，T为sqrt(M)
void work() {
  fac[1][0] = fac[1][1] = fac[1][2] = 1;
  for (int i = 2; i <= M; ++i) {
    if (!isp[i]) {
      fac[i][0] = fac[i][1] = 1;
      fac[i][2] = i;
      pri[++tot] = i;
    }
    for (int j = 1; pri[j] * i <= M; ++j) {
      int tmp = pri[j] * i;
      isp[tmp] = true;
      fac[tmp][0] = fac[i][0] * pri[j];
      fac[tmp][1] = fac[i][1];
      fac[tmp][2] = fac[i][2];
      if (fac[tmp][0] > fac[tmp][1]) {
        fac[tmp][0] ^= fac[tmp][1] ^= fac[tmp][0] ^= fac[tmp][1];
      }
      if (fac[tmp][1] > fac[tmp][2]) {
        fac[tmp][1] ^= fac[tmp][2] ^= fac[tmp][1] ^= fac[tmp][2];
      }
// 对于整数运算a ^= b ^= a ^= b等价于swap(a, b)这里就是手动进行length = 3的排序
      if (i % pri[j] == 0) {
        break;
      }
    }
  }
  for (int i = 0; i <= T; ++i) {
    pre[0][i] = pre[i][0] = i;
  }
  for (int i = 1; i <= T; ++i) {
    for (int j = 1; j <= i; ++j) {
      pre[i][j] = pre[j][i] = pre[j][i % j];
    }
  }
}

查询

若当前枚举的为$a$，只需注意$a>\sqrt{N}$时分$a\mid y$和$a\nmid y$两种情况即可

以下为代码

int gcd(int a, int b) {
// 不想写if-else所以用三目运算符代替并缩进了一下
  int ans = 1;
  for (int i = 0; i < 3; ++i) {
    int tmp = (fac[a][i] > T) ?
                (b % fac[a][i] ?
                   1
                 : fac[a][i]
                )
              : pre[fac[a][i]][b % fac[a][i]];
    b /= tmp;
    ans *= tmp;
  }
  return ans;
}

本题做法

基本上已经没了，就是求$\gcd$然后按题目给的算，放个代码吧

const int N = 5000;
const int M = 1000000;
const int T = 1000;
const int Mod = 998244353;

void work();
int gcd(int, int);

int pre[T + 2][T + 2];
int a[N + 2], b[N + 2];
int fac[M + 2][3];
bool isp[M + 2];
int pri[M / 10], tot;
int n;

int main () {
  work();
  read(n);
  for (int i = 1; i <= n; ++i) {
    read(a[i]);
  }
  for (int i = 1; i <= n; ++i) {
    read(b[i]);
  }
  for (int i = 1; i <= n; ++i) {
    int ans = 0;
    for (int j = 1, now = i; j <= n; ++j, now = now * 1ll * i % Mod) {
      ans = (ans + int(now * 1ll * gcd(a[i], b[j]) % Mod)) % Mod;
    }
    write(ans), EL;
  }
  return 0;
}

void work() {
  fac[1][0] = fac[1][1] = fac[1][2] = 1;
  for (int i = 2; i <= M; ++i) {
    if (!isp[i]) {
      fac[i][0] = fac[i][1] = 1;
      fac[i][2] = i;
      pri[++tot] = i;
    }
    for (int j = 1; pri[j] * i <= M; ++j) {
      int tmp = pri[j] * i;
      isp[tmp] = true;
      fac[tmp][0] = fac[i][0] * pri[j];
      fac[tmp][1] = fac[i][1];
      fac[tmp][2] = fac[i][2];
      if (fac[tmp][0] > fac[tmp][1]) {
        fac[tmp][0] ^= fac[tmp][1] ^= fac[tmp][0] ^= fac[tmp][1];
      }
      if (fac[tmp][1] > fac[tmp][2]) {
        fac[tmp][1] ^= fac[tmp][2] ^= fac[tmp][1] ^= fac[tmp][2];
      }
      if (i % pri[j] == 0) {
        break;
      }
    }
  }
  for (int i = 0; i <= T; ++i) {
    pre[0][i] = pre[i][0] = i;
  }
  for (int i = 1; i <= T; ++i) {
    for (int j = 1; j <= i; ++j) {
      pre[i][j] = pre[j][i] = pre[j][i % j];
    }
  }
}
int gcd(int a, int b) {
  int ans = 1;
  for (int i = 0; i < 3; ++i) {
    int tmp = (fac[a][i] > T) ?
                (b % fac[a][i] ?
                   1
                 : fac[a][i]
                )
              : pre[fac[a][i]][b % fac[a][i]];
    b /= tmp;
    ans *= tmp;
  }
  return ans;
}