【CF571E】Geometric Progressions（分类讨论细节题）

点此看题面

• 给定\(n\)个序列，第\(i\)个序列为\(\{a_i,a_ib_i,a_ib_i^2,a_ib_i^3...\}\)。
• 问最小的出现在所有序列中的数（模\(10^9+7\)）。
• \(n\le100,a_i,b_i\le10^9\)

质因数分解

显然我们第一步是把每个数质因数分解掉，这样乘法就变成了幂次的加法。

然后考虑我们维护前\(i-1\)个数的已知答案\(P\)和最小的数\(Q\)满足其为所有\(b\)的幂，然后尝试用第\(i\)个元素\(A,B\)去更新\(P,Q\)。

那么我们就是要找到两个最小自然数\(x,y\)，使得\(PQ^x=AB^y\)。

于是就迎来了一大波分类讨论。

我们枚举每一个质因子，分为下面几类情况：

• \(B,Q\)中都不存在，则除非\(A=P\)，否则显然无解。
• \(B,Q\)中某一个存在（假设是\(B\)中存在），则我们可以直接通过\(A,P\)中该质因子个数的差值除以\(B\)中该质因子的个数，得到\(x\)，然后另找一个\(Q\)中存在的质因子就可以计算出\(y\)了（如果\(Q\)没有质因子，即\(Q=1\)，显然\(y\)是多少都无所谓）。
• \(B,Q\)中都存在，这类才是最麻烦的情况，下面会重点讨论。

显然有第二类时就可以直接过掉了，否则我们需要考虑上面的第三类情况。

\(B,Q\)中都存在的质因子

对于每个\(B,Q\)中都存在的质因子我们可以列出一个关于\(x,y\)的二元一次方程，先去重消去其中等价的那些。

如果剩下方程超过一个，可以直接任取两个解二元一次方程组，然后代回去检验一下。

如果只剩一个，众所周知可以用\(exgcd\)解出二元一次方程的一组解，然后把它们转化成最小自然数解就可以了。

口胡起来就这么简单，写起来有点复杂。。。

代码：\(O(n\sqrt VlogV)\)

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 100
#define S 50000
#define X 1000000007
#define int long long
#define NA() (puts("-1"),exit(0),0)//无解，输出-1并直接结束程序
using namespace std;
int n,a[N+5],b[N+5],A[S+5],B[S+5],P[S+5],Q[S+5];
I int QP(RI x,RI y) {RI t=1;W(y) y&1&&(t=t*x%X),x=x*x%X,y>>=1;return t;}
I int gcd(CI x,CI y) {return y?gcd(y,x%y):x;} 
I void exgcd(CI x,CI y,int& a,int& b) {y?(exgcd(y,x%y,b,a),b-=x/y*a):(a=1,b=0);}//exgcd解二元一次方程
namespace Prime
{
	int Pt,Pr[S+5];I int IsP(CI x) {for(RI i=2;i*i<=x;++i) if(!(x%i)) return 0;return 1;}//判质数
	int cnt;I void Init() {for(RI i=2;i<=31622;++i) IsP(i)&&(Pr[++Pt]=i);cnt=Pt;}//预处理小于等于sqrt(1e9)的质数
	map<int,int> id;I int ID(CI x) {return !id[x]&&(Pr[id[x]=++cnt]=x),id[x];}//给大于sqrt(1e9)的质数标号
	I void Work(int* V,RI x) {for(RI i=1;i<=Pt;++i) W(!(x%Pr[i])) x/=Pr[i],++V[i];x^1&&++V[ID(x)];}//分解质因数
}
namespace Equation//解方程
{
	int tot;struct Data
	{
		int a,b,c;I Data(CI x=0,CI y=0,CI z=0):a(x),b(y),c(z){}//二元一次方程ax+by=c
		I bool operator < (Con Data& o) Con {return a^o.a?a<o.a:(b^o.b?b<o.b:c<o.c);}
	}s[S+5];set<Data> vis;
	I void Add(CI a,CI b,CI c)//新增一个方程
	{
		RI g=abs(gcd(a,b));c%g&&NA();Data t(a/g,b/g,c/g);//同除以公约数化简
		!vis.count(t)&&(vis.insert(s[++tot]=t),0);//去重
	}
	I void Solve(int& x,int& y)//解二元一次方程组
	{
		RI i,p=s[1].c*s[2].b-s[1].b*s[2].c,q=s[1].a*s[2].b-s[1].b*s[2].a;
		(!q||p%q)&&NA(),x=p/q,(s[1].c-s[1].a*x)%s[1].b&&NA(),y=(s[1].c-s[1].a*x)/s[1].b;//解出x,y
		for(i=3;i<=tot;++i) (s[1].a*x+s[1].b*y)^s[1].c&&NA();//代入检验
	}
}
signed main()
{
	using namespace Prime;using namespace Equation;
	RI i,j,x,y;for(scanf("%lld",&n),i=1;i<=n;++i) scanf("%lld%lld",a+i,b+i);
	RI u,v,t,fg;for(Init(),i=1;i<=n;++i)
	{
		memset(A,0,sizeof(A)),memset(B,0,sizeof(B)),Work(A,a[i]),Work(B,b[i]);//质因数分解
		if(i==1) {for(j=1;j<=Pt+cnt;++j) P[j]=A[j],Q[j]=B[j];continue;}//第一个直接存储并跳过
		for(u=v=-1,fg=0,j=1;j<=Pt+cnt;++j) if(B[j]||Q[j])
		{
			fg=1;if(!B[j]) {((t=A[j]-P[j])<0||t%Q[j])&&NA(),v=t/Q[j];break;}
			if(!Q[j]) {((t=P[j]-A[j])<0||t%B[j])&&NA(),u=t/B[j];break;}
		}else A[j]^P[j]&&NA();
		if(!fg) continue;if(~u||~v)//如果存在第二类
		{
			if(!~u) for(j=1;j<=Pt+cnt;++j) if(B[j])
				{((t=P[j]+v*Q[j]-A[j])<0||t%B[j])&&NA(),u=t/B[j];break;}
			if(!~v) for(j=1;j<=Pt+cnt;++j) if(Q[j])
				{((t=A[j]+u*B[j]-P[j])<0||t%Q[j])&&NA(),v=t/Q[j];break;}
			for(j=1;j<=Pt+cnt;++j) (A[j]+u*B[j])^(P[j]+=v*Q[j])&&NA();goto End;//更新P，同时检验
		}
		for(tot=0,vis.clear(),j=1;j<=Pt+cnt;++j) (B[j]||Q[j])&&(Add(B[j],-Q[j],P[j]-A[j]),0);//列出方程
		if(tot==1)//只有一个方程
		{
			exgcd(s[1].a,s[1].b*=-1,u,v),u*=s[1].c,v*=-s[1].c;//exgcd解出一组解
			t=max(u<0?(-u-1)/s[1].b+1:0,v<0?(-v-1)/s[1].a+1:0),u+=t*s[1].b,v+=t*s[1].a;//使其非负
			(u<0||v<0)&&NA(),t=min(u/s[1].b,v/s[1].a),u-=t*s[1].b,v-=t*s[1].a;//使其最小
			for(j=1;j<=Pt+cnt;++j) P[j]+=v*Q[j];goto End;//更新P
		}
		for(Solve(u,v),j=1;j<=Pt+cnt;++j) P[j]+=v*Q[j];//更新P
		End:for(j=1;j<=Pt+cnt;++j) (B[j]||Q[j])&&(Q[j]*=B[j]/gcd(B[j],Q[j]));//更新Q
	}
	for(t=i=1;i<=Pt+cnt;++i) (t*=QP(Pr[i],P[i]))%=X;return printf("%lld\n",t),0;//计算最终答案
}