[CF571E] Geometric Progressions——数论、质因数分解

CF571E Geometric Progressions

题解

我几乎是完全看着这篇题解打的。

乘积很大不方便处理，所以我们把所有 $a_i,b_i$ 都做质因数分解。

然后考虑合并两个数列，结果比较显然：要么无解，要么得到一个数，要么得到一个新的等比数列。我们的任务是对 $n$ 个数列两两合并，最后输出结果（等比数列就输出第一项）。这题中为了方便，我们可以把一个数看成是公比为1的等比数列。

记 $c_p(x)$ 表示 $x$ 的质因数分解中质数 $p$ 的次数，那么我们合并两个数列 $a\ast b^{?}$ 和 $c\ast d^{?}$ （目的是找到自然数系数 $x,y$ 使得 $a\ast b^x=c\ast d^y$）时，可以枚举每个质因子 $p$，然后分情况讨论：

• $c_p(b)$ 和 $c_p(d)$ 都等于0，此时如果 $c_p(a)\ne c_p(c)$ 必定无解，否则不管它；
• $c_p(b)$ 和 $c_p(d)$ 两者其一等于0，此时非零的那边的系数已经确定（也可能无自然数解，那么无解），为零的那边系数待定；
• $c_p(b)$ 和 $c_p(d)$ 都非0，此时得到二元一次方程 $c_p(b)x-c_p(d)y=c_p(c)-c_p(a)$，先记录下来放一边。

处理完了后，会出现几种情况：

• 无解，那么不用管了；
• 两个系数都没确定，此时若不是两个数列都只有一个值，那么一定剩下了若干个二元一次方程。把这些方程中线性相关的去重过后，剩下如果有两个及以上方程，那么系数有唯一解，代入判断即可；若只有一个方程，那么需要用扩展欧几里得算法解出 $x,y$ 都为非负整数，且 $x$ 或 $y$ 最小的解，然后可以合并出一个新的等比数列。设这个新的等比数列为 $a^{\prime }\ast b^{\prime ?}$，那么有$$\begin{gathered} c_p(a^{\prime })=c_p(a)+x\times c_p(b) \\ {c}_p(b^{\prime })=\text{lcm}(c_p(b),c_p(d)) \end{gathered}$$
• 其中一个系数定了，那么可以把某个数列处理成一个数，然后再重复上面过程。这个时候枚举质因子值可能出现1、2两种情况，所以最后两个系数都能确定，然后接着讨论：
• 两个系数都定了，那么带进去检验即可。

为了方便实现去重，我用了大常数 map。最后根据质因数分解输出模意义下 $a$ 的值即可。

代码看起来长，其实有大段都是复读机。

代码

#include<bits/stdc++.h>//JZM yyds!!
#define ll long long
#define lll __int128
#define uns unsigned
#define fi first
#define se second
#define IF (it->fi)
#define IS (it->se)
#define END putchar('\n')
#define lowbit(x) ((x)&-(x))
#define inline jzmyyds
using namespace std;
const int MAXN=114514;
const ll INF=1e18;
ll read(){
	ll x=0;bool f=1;char s=getchar();
	while((s<'0'||s>'9')&&s>0){if(s=='-')f^=1;s=getchar();}
	while(s>='0'&&s<='9')x=(x<<1)+(x<<3)+(s^48),s=getchar();
	return f?x:-x;
}
int readuns(){
	int x=0;char s=getchar();
	while(s<'0'||s>'9')s=getchar();
	while(s>='0'&&s<='9')x=(x<<1)+(x<<3)+(s^48),s=getchar();
	return x;
}
int ptf[50],lpt;
void print(ll x,char c='\n'){
	if(x<0)putchar('-'),x=-x;
	ptf[lpt=1]=x%10;
	while(x>9)x/=10,ptf[++lpt]=x%10;
	while(lpt>0)putchar(ptf[lpt--]^48);
	if(c>0)putchar(c);
}
const ll MOD=1e9+7;
ll ksm(ll a,ll b,ll mo){
	ll res=1;
	for(;b;b>>=1,a=a*a%mo)if(b&1)res=res*a%mo;
	return res;
}
ll exgcd(ll a,ll b,ll&x,ll&y){
	if(!b)return x=1,y=0,a;
	ll d=exgcd(b,a%b,y,x);
	return y-=a/b*x,d;
}
ll gcd(ll a,ll b){return !b?a:gcd(b,a%b);}
map<ll,ll> getp(ll x){
	map<ll,ll>res;
	for(ll i=2;i*i<=x;i++)if(x%i==0){
		ll&d=res[i]=0;
		while(x%i==0)x/=i,d++;
	}if(x>1)res[x]=1;
	return res;
}
#define pll pair<ll,ll>
struct itn{
	ll a,b,c;itn(){}
	itn(ll A,ll B,ll C){
		ll d=gcd(gcd(A,B),C);
		if(A<0)d=-d;
		if(d)a=A/d,b=B/d,c=C/d;
		else a=b=c=0;
	}
	bool operator<(const itn&y)const{
		return a<y.a||(a==y.a&&(b<y.b||(b==y.b&&c<y.c)));
	}
	pll rsl()const{
		ll x,y,d=exgcd(a,abs(b),x,y),p=b/d,q=a/d;
		if(b<0)y=-y;
		if(c%d)return pll(-1,-1);
		x*=c/d,y*=c/d;
		if(q<0)q=-q,p=-p;
		if(y<0)x-=p*((-y+q-1)/q),y+=q*((-y+q-1)/q);
		if(p<0)p=-p,q=-q;
		if(x<0)y-=q*((-x+p-1)/p),x+=p*((-x+p-1)/p);
		if(x>0){
			if(q>0)y+=q*(x/p),x-=p*(x/p);
			else{
				q=-q;ll md=min(x/p,y/q);
				x-=p*md,y-=q*md;
			}
		}
		if(x<0||y<0)return pll(-1,-1);
		return pll(x,y);
	}
};
pll getcro(const itn&p,const itn&q){
	ll D=p.a*q.b-p.b*q.a;
	if(!D)return pll(-1,-1);
	ll D1=p.c*q.b-p.b*q.c,D2=p.a*q.c-p.c*q.a;
	if(D<0)D=-D,D1=-D1,D2=-D2;
	if(D1<0||D2<0||D1%D||D2%D)return pll(-1,-1);
	return pll(D1/D,D2/D);
}
map<ll,ll>a,b;
bool ok=1,hv=0;
int main()
{
	for(int N=read();N--;){
		ll A=read(),B=read();
		if(!ok)continue;
		map<ll,ll>c=getp(A),d=getp(B);
		for(auto it:d)c[it.fi];
		if(!hv){swap(a,c),swap(b,d),hv=1;continue;}
		for(auto it:c)a[it.fi];
		set<itn>st;
		ll u=-1,v=-1;
		for(auto mmp:a){
			const ll p=mmp.fi,ci=b[p],cj=d[p];
			// printf("\n %lld  %lld %lld %lld %lld\n",p,a[p],c[p],ci,cj);
			if(!ci&&!cj){
				if(a[p]==c[p])continue;
				ok=0;break;
			}if(!ci){
				ll dt=a[p]-c[p];
				if(dt<0||dt%cj){ok=0;break;}
				v=dt/cj;continue;
			}if(!cj){
				ll dt=c[p]-a[p];
				if(dt<0||dt%ci){ok=0;break;}
				u=dt/ci;continue;
			}st.insert(itn(ci,-cj,c[p]-a[p]));
		}
		if(!ok)continue;
		if(u<0&&v<0){
			if(st.empty())continue;
			if(st.size()==1){
				pll rs=(*st.begin()).rsl();
				if(rs.fi<0){ok=0;continue;}
				for(auto mmp:a){
					const ll p=mmp.fi,ci=b[p],cj=d[p];
					a[p]+=rs.fi*ci;
					if(!ci||!cj)b[p]=0;
					else b[p]=ci/gcd(ci,cj)*cj;
				}continue;
			}
			itn f=*st.begin(),g=*st.rbegin();
			pll rs=getcro(f,g);
			if(rs.fi<0){ok=0;continue;}
			for(auto x:st)if(rs.fi*x.a+rs.se*x.b!=x.c)ok=0;
			if(!ok)continue;
			u=rs.fi,v=rs.se;
		}
		if(u<0||v<0){
			for(auto mmp:a){
				const ll p=mmp.fi;
				if(u>=0)a[p]+=u*b[p],b[p]=0;
				if(v>=0)c[p]+=v*d[p],d[p]=0;
			}u=v=0;
			for(auto mmp:a){
				const ll p=mmp.fi,ci=b[p],cj=d[p];
				if(!ci&&!cj){
					if(a[p]==c[p])continue;
					ok=0;break;
				}if(!ci){
					ll dt=a[p]-c[p];
					if(dt<0||dt%cj){ok=0;break;}
					v=dt/cj;continue;
				}if(!cj){
					ll dt=c[p]-a[p];
					if(dt<0||dt%ci){ok=0;break;}
					u=dt/ci;continue;
				}
			}if(!ok)continue;
		}
		if(u>=0&&v>=0){
			for(auto mmp:a){
				const ll p=mmp.fi;
				a[p]+=u*b[p],c[p]+=v*d[p],b[p]=0;
				if(a[p]^c[p]){ok=0;break;}
			}
			if(!ok)continue;
		}
	}ll ans=1;
	if(!ok)return printf("-1\n"),0;
	for(auto it:a)(ans*=ksm(it.fi,it.se,MOD))%=MOD;
	print(ans);
	return 0;
}