基于值域预处理的快速 GCD
P5435 基于值域预处理的快速 GCD
这道题是快速gcd的模板题，首先题中数据如果直接对于每一个\(gcd(a_i,b_j)\)都用欧几里得算法很明显是会\(T\)的，所以我们就需要快速的\(gcd\)了。
首先感谢洛谷上的@クトリ大佬，我是看他的题解学会的快速\(gcd\)的，本文有部分证明不会讲到，主要我不会推，可以在クトリ大佬的题解中看到。
（本文中的大多数情况是以学会这种算法为主，所以并不会讲其所以然，也就是为什么这样是对的，即证明过程省略，所以定理就当作是对的即可。）
正式题解
前置知识：欧几里得算法，快速筛
引入
我们计算很多最大公约数时，每次使用欧几里得的时间复杂度是很高的，那么我们可不可以减少调用欧几里得的次数，很多算法都是通过提前计算一些值来优化算法的计算速度，快速gcd也是这样的。
一步一步看，很多过程中不理解的东西，会慢慢解释
定理1
我们需要求\(gcd(x,y)\),可以把\(x\)分解，这样我们可以通过提前计算小部分的最大公约数存储在\(gcd[x][y]\)中，这样就可以减少调用欧几里得算法的次数。
假设要询问\(gcd(x,y)\),将\(x\)分解为\(x_1,x_2...x_n\),先算\(gcd(x_1,y)\),然后加到\(ans\)中，再将\(y\)除以\(gcd(x_1,y)\),因为这个因子已经计算，不能重复计数，然后是\(x_2\)直到全部的贡献都算完然后加到\(ans\)中。过程大致如下：
int ans=1;
for(int i=1;i<=n;++i){
int tmp=gcd[fac[x][i]][b%fac[x][i]];//gcd为预处理的最大公约数，fac[x][i]为x的第i个因子
y/=tmp;
ans*=tmp;
}
cout<<ans<<endl;
定理2
很显然，如果是这样计算我们就必须要考虑两个问题，gcd数组需要算多少数对的最大公约数，x应该分解成几项（分解的过多会劣化时间复杂度，也会使程序更加复杂）。
任何一个\(x∈{1,n}\)都可以分解成\({a,b,c},且{a,b,c}∈{1,\sqrt x}或者∈Prime\)，Prime是素数。
分解的实现过程
基本就是线性筛的过程中添加了一些步骤，在线性筛标记合数时，存储它的因子,需要注意的是在这个过程中保持\(a<b<c\),这个过程描述不清楚，直接看代码吧。
void init(){
num[0]=num[1]=1;
fen[1][1]=fen[1][2]=fen[1][3]=1;
for(int i=2;i<=M;++i){
if(num[i]==0){
fen[i][1]=fen[i][2]=1;
fen[i][3]=i;
prime[++cnt]=i;
}
for(int j=1;j<=cnt&&i*prime[j]<=M;++j){
int x=i*prime[j];
num[x]=1;
fen[x][1]=fen[i][1]*prime[j];
fen[x][2]=fen[i][2];
fen[x][3]=fen[i][3];
if(fen[x][1]>fen[x][2])swap(fen[x][1],fen[x][2]);
if(fen[x][2]>fen[x][3])swap(fen[x][2],fen[x][3]);
if(i%prime[j]==0)break;
}
}
}
查询的实现过程
上边提到，提前计算gcd时是有上界的，并且尽可能的小，其实计算到\(gcd[\sqrt x][\sqrt x]\)即可，但是有一个问题，我们分解的因子中可能存在素数，这个地方就很高明，因为他是素数，所以根本不需要再次欧几里得，我们取余判断是否可以整除，就能算出这个因子的贡献，因为如果\(y\)可以被这个素数因子整除，那么说明\(x,y\)都有这个因子，它的贡献就是1，否则说明，\(y\)就没有这个因子，那么最大公约数的因子必定有这个素数因子。
inline long long get_gcd(int x,int y){
long long res=1ll;
for(int i=1,tmp=1;i<=3;++i){
if(fen[x][i]>T){
if(y%fen[x][i])tmp=1;
else tmp=fen[x][i];
}
else tmp=gcd[fen[x][i]][y%fen[x][i]];
y/=tmp,res*=tmp;
}
return res;
}
最后需要注意的点，在求最后的答案时是需要开longlong的，但是不能直接#define int long long 会爆内存，别问我为什么知道
完整代码
#include<bits/stdc++.h>
using namespace std;
const int N=5e3+2;
const int M=1e6+2;
const int T=1e3+2;
const int mod=998244353 ;
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch))x=(x<<3)+(x<<1)+(ch^'0'),ch=getchar();
return x*f;
}
int n,a[N],b[N],gcd[T+1][T+1];
int num[M],prime[M/10],fen[M][4],cnt;
inline int GCD(int x,int y){
if(y==0)return x;
else return GCD(y,x%y);
}
void init(){
num[0]=num[1]=1;
fen[1][1]=fen[1][2]=fen[1][3]=1;
for(int i=2;i<=M;++i){
if(num[i]==0){
fen[i][1]=fen[i][2]=1;
fen[i][3]=i;
prime[++cnt]=i;
}
for(int j=1;j<=cnt&&i*prime[j]<=M;++j){
int x=i*prime[j];
num[x]=1;
fen[x][1]=fen[i][1]*prime[j];
fen[x][2]=fen[i][2];
fen[x][3]=fen[i][3];
if(fen[x][1]>fen[x][2])swap(fen[x][1],fen[x][2]);
if(fen[x][2]>fen[x][3])swap(fen[x][2],fen[x][3]);
if(i%prime[j]==0)break;
}
}
for(int i=0;i<=T;++i){
for(int j=i;j<=T;++j){
gcd[i][j]=gcd[j][i]=GCD(i,j);
}
}
}
inline long long get_gcd(int x,int y){
long long res=1ll;
for(int i=1,tmp=1;i<=3;++i){
if(fen[x][i]>T){
if(y%fen[x][i])tmp=1;
else tmp=fen[x][i];
}
else tmp=gcd[fen[x][i]][y%fen[x][i]];
y/=tmp,res*=tmp;
}
return res;
}
int main(){
n=read();
for(int i=1;i<=n;++i)a[i]=read();
for(int i=1;i<=n;++i)b[i]=read();
init();
for(int i=1;i<=n;++i){
long long ans=0,tmp=i;
for(int j=1;j<=n;++j,tmp=tmp*i%mod){
ans=(ans+tmp*get_gcd(a[i],b[j])%mod)%mod;
}
printf("%d\n",ans%mod);
}
return 0;
}
`
浙公网安备 33010602011771号