题解 P5435 【基于值域预处理的快速 GCD】

moongazer · 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的指正

一些约定

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.

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;
}