【洛谷P5435】基于值域预处理的快速 GCD

题目

题目链接：https://www.luogu.com.cn/problem/P5435
给定 \(n\) 个正整数 \(a_1,a_2,\dots,a_n\)，再给定 \(n\) 个正整数\(b_1,b_2,\dots,b_n\)，你需要对每对 \((i,j)\) 求出 \(a_i\) 与 \(b_j\) 的最大公因数。
不难发现你的输出应有 \(n^2\) 个正整数。为了减少输出对程序的运行效率的影响，你只需要输出 \(n\) 行，每行一个整数 \(A_i\)。
其中对于 \(i\in[1,n]\)，\(A_i=\sum_{j=1}^{n}i^j\gcd(a_i,b_j)\)。由于答案可能过大，你只需要输出模 \(998,244,353\) 后的结果即可。
\(n\leq 5000\)，\(1\leq a_i\leq 10^6\)。

思路

一种挺有趣的科技。
对于任意一个数 \(x\)，肯定有一种方法可以把 \(x\) 分解成 \(x_1\times x_2\times x_3\)，且 \(x_1,x_2,x_3\) 都 \(\leq \sqrt{x}\) 或是质数。
具体的构造方法，我们可以线性筛出值域内所有数的最小质因子。对于一个数 \(p\)，如果它是质数，显然可以分解为 \(1\times 1\times p\)；如果它不是质数，设它的最小质因子是 \(d\)，且 \(q=\frac{p}{d}\) 分解为了 \(q_1\times q_2\times q_3\)（\(q_1\leq q_2\leq q_3\)），那么可以把 \(p\) 分解为 \(q_1\times d,q_2,q_3\)。
证明：

• 如果 \(d\leq \sqrt[4]{p}\)，因为 \(q_1\leq \sqrt[3]{\frac{p}{d}}\)，所以 \(d\times q_1\leq d\times \sqrt[3]{\frac{p}{d}}=p^{\frac{1}{3}}\times d^{\frac{2}{3}}\)。而 \(\sqrt{p}\div p^{\frac{1}{3}}=p^{\frac{1}{6}}\)，\(d^{\frac{2}{3}}\leq \left(\sqrt[4]{p}\right)^{\frac{2}{3}}\leq p^{\frac{1}{6}}\)。故 \(d\times q_1\leq \sqrt{p}\)。
• 如果 \(d>\sqrt[4]{p}\)
  • 如果 \(q_1=1\)，显然 \(1\times d\leq \sqrt{p}\)。
  • 如果 \(q_1>1\)，则 \(d\leq q_1\leq q_2\leq q_3\)，则 \(d\times q_1\times q_2\times q_3>\left(\sqrt[4]{p}\right)^4=p\)，矛盾。

这样把值域内所有的数分解好后，预处理出一个 \(\sqrt{v}\times \sqrt{v}\) 的 \(\gcd\) 数组（\(v\) 是值域）。对于每一次询问 \(a,b\)，把 \(a\) 分解成 \(a_1\times a_2\times a_3\)。

• 如果 \(a_1\) 是质数：如果 \(b\) 是 \(a_1\) 的倍数返回 \(a_1\)，否则返回 \(1\)。
• 如果 \(a_1\) 不是质数，\(\gcd(a_1,b)=\gcd(a_1,b\bmod a_1)\)，此时两个数都不超过 \(\sqrt{v}\)，直接返回预处理的值。

对 \(a_1,a_2,a_3\) 分别做一次乘起来，就是 \(\gcd(a,b)\) 了。注意每次 \(b\) 需要除掉返回的 \(\gcd\)。
时间复杂度 \(O(v+n^2)\)。

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int N=5010,M=1000010,MM=1010,MOD=998244353;
int n,m,a[N],b[N],prm[M],v[M],f[M][3],g[MM][MM];

void findprm(int n)
{
	v[1]=1;
	for (int i=2;i<=n;i++)
	{
		if (!v[i]) prm[++m]=i,v[i]=i;
		for (int j=1;j<=m;j++)
		{
			if (i>n/prm[j] || prm[j]>v[i]) break;
			v[i*prm[j]]=prm[j];
		}
	}
}

void prework()
{
	for (int i=1;i<MM;i++) g[0][i]=g[i][0]=i;
	for (int i=1;i<MM;i++)
		for (int j=1;j<=i;j++)
			g[i][j]=g[j][i]=g[j][i%j];
	for (int i=1;i<M;i++)
		if (v[i]==i)
			f[i][0]=1,f[i][1]=1,f[i][2]=i;
		else
		{
			f[i][0]=f[i/v[i]][0]*v[i]; f[i][1]=f[i/v[i]][1]; f[i][2]=f[i/v[i]][2];
			if (f[i][0]>f[i][1]) swap(f[i][0],f[i][1]);
			if (f[i][1]>f[i][2]) swap(f[i][1],f[i][2]);
		}
}

int gcd(int x,int &y)
{
	int d=(x>MM)?((y%x)?1:x):g[x][y%x];
	y/=d;
	return d;
}

int main()
{
	findprm(M-1);
	prework();
	scanf("%d",&n);
	for (int i=1;i<=n;i++) scanf("%d",&a[i]);
	for (int i=1;i<=n;i++) scanf("%d",&b[i]);
	for (int i=1;i<=n;i++)
	{
		ll ans=0,res=1;
		for (int j=1,k;j<=n;j++)
		{
			res=res*i%MOD; k=b[j];
			ans=(ans+res*gcd(f[a[i]][0],k)*gcd(f[a[i]][1],k)*gcd(f[a[i]][2],k))%MOD;
		}
		cout<<ans<<"\n";
	}
	return 0;
}