【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;//计算最终答案
}
待到再迷茫时回头望,所有脚印会发出光芒