Luogu P5435 基于值域预处理的快速 GCD

突然想起来曾经偶然在陈指导的博客看到过这个 $O(1)$ 做 $\gcd$ 的方法，其实理解了之后还是比较简单的，以下设数的值域为 $S$

首先我们定义对于一个数 $x$ 的一种分解方式为三元组 $(a,b,c)$ ，满足 $a\le b\le c\and a\times b\times c=x$ ，并且若 $c$ 为合数则 $c\le \sqrt x$

证明的话也很简单，考虑若 $c>\sqrt x$ ，则 $a\times b<\sqrt x$ ，考虑找到 $c$ 的一种合法分解 $c=d\times e$ 满足 $d\le e\le \sqrt x$ 那么必然可以把三元组调整为 $(d,a\times b,e)$

有了这个性质我们就很好办了，考虑询问 $\gcd(x,y)$ 的话我们可以分别考虑 $x$ 的分解 $(a,b,c)$ 和 $y$ 的贡献

$a$ 和 $y$ 的贡献就是将答案乘上 $\gcd(a,y)$ ，然后将 $y$ 除以 $\gcd(a,y)$ 来消去这个因子的贡献，再对 $b,c$ 做同样的事即可

由于一般情况下 $a,b,c\le \sqrt x$ ，因此我们可以记录下 $\sqrt S$ 范围内的 $\gcd$ ，这样就可以 $O(1)$ 处理了，当然质数的情况可以简单讨论掉

现在就要考虑怎么 $O(S)$ 求所有数的分解了，其实很简单，只要在欧拉筛的时候通过 $x$ 的最小质因子 $p$

设 $\frac{x}{p}$ 的分解为 $(a_0,b_0,c_0)$ ，那么 $x$ 的合法分解为 $\{a_0\times p,b_0,c_0\}$ 的不降排序，考虑证明：

当 $x$ 为质数是分解显然成立
当 $p\le \sqrt[4] x$ 时，将 $a_0\le \sqrt[3]{\frac{x}{p}}$ 代入有 $a_0\times p\le \sqrt x$
当 $p>\sqrt[4]x$ 时，分情况考虑：
- 若 $a_0=1$ ，显然 $a_0\times p\le \sqrt x$
- 若 $a_0\ne 1$ ，由于 $x$ 不是质数，那么 $\frac{x}{p}$ 的最小质因子 $q$ 就是 $x$ 的第二小质因子，一定 $\ge p$ 。而 $a_0,b_0,c_0$ 都为 $\frac{x}{p}$ 的非 $1$ 因子，故有 $p\le q\le a_0\le b_0\le c_0,p\times a_0\times b_0\times c_0>(\sqrt[4]x)^4=x$ ，矛盾。故不存在这种情况

因此代码就很简单了

#include<cstdio>
#include<algorithm>
#define RI register int
#define CI const int&
using namespace std;
const int N=1000005,S=1000,mod=998244353;
int n,a[N],b[N],prm[N],cnt,G[S+5][S+5],frac[N][3]; bool vis[N];
inline void init(CI n)
{
	RI i,j; for (i=0;i<3;++i) frac[1][i]=1;
	for (i=2;i<=n;++i)
	{
		if (!vis[i]) for (prm[++cnt]=frac[i][2]=i,j=0;j<2;++j) frac[i][j]=1;
		for (j=1;j<=cnt&&i*prm[j]<=n;++j)
		{
			vis[i*prm[j]]=1; for (RI k=0;k<3;++k) frac[i*prm[j]][k]=frac[i][k];
			frac[i*prm[j]][0]*=prm[j]; sort(frac[i*prm[j]],frac[i*prm[j]]+3);
			if (i%prm[j]==0) break;
		}
	}
	for (i=1;i<=S;++i) G[i][0]=G[0][i]=i;
	for (i=1;i<=S;++i) for (j=1;j<=i;++j) G[i][j]=G[j][i]=G[i%j][j];
}
inline int gcd(int x,int y,int ret=1)
{
	for (RI i=0;i<3;++i)
	{
		int cur; if (frac[x][i]>S)
		{
			if (y%frac[x][i]) cur=1; else cur=frac[x][i];
		} else cur=G[frac[x][i]][y%frac[x][i]];
		y/=cur; ret*=cur;
	}
	return ret;
}
int main()
{
	//freopen("CODE.in","r",stdin); freopen("CODE.out","w",stdout);
	RI i,j; for (init(1000000),scanf("%d",&n),i=1;i<=n;++i) scanf("%d",&a[i]);
	for (i=1;i<=n;++i) scanf("%d",&b[i]); for (i=1;i<=n;++i)
	{
		int ans=0,cur=i; for (j=1;j<=n;++j,cur=1LL*cur*i%mod)
		(ans+=1LL*cur*gcd(a[i],b[j])%mod)%=mod; printf("%d\n",ans);
	}
	return 0;
}