#緑以下B #中国剰余定理 ''' ■参考文献 (1)中国剰余定理 (CRT) の解説と、それを用いる問題のまとめ https://qiita.com/drken/items/ae02240cd1f8edfc86fd (2)けんちょんさんの中国剰余定理の補足 https://sen-comp.hatenablog.com/entry/2021/03/05/194205 (3)【数論】オイラーの定理とその2通りの証明 https://mathlandscape.com/euler-theorem/ (4)任意modでの畳み込み演算をO(n log(n))で https://math314.hateblo.jp/entry/2015/05/07/014908 (5)もう一度中国剰余定理 https://joho-ka.mints.ne.jp/chinese-remainder-theorem2?doing_wp_cron=1677505078.4430840015411376953125 (6)拡張互除法と中国剰余定理(Garnerのアルゴリズム) https://www.creativ.xyz/ect-gcd-crt-garner-927/ (7)Algorithms with Python : 拡張ユークリッドの互除法 http://www.nct9.ne.jp/m_hiroi/light/pyalgo70.html ''' #①素因数の振り分け def CRT_SnukeDistribute(A,B): #N≡P mod A≡Q mod B, AとBが互いに素になるよう振り直し G=Euclid(A,B)[0] A1,B1=A//G,B//G #A,Bそれぞれに固有の素因数 A2=Euclid(A1,G)[0] #G: 最大公約数の割り振りを決める B2=G//A2 G1=Euclid(A2,B2)[0] while G1>1: A2,B2=A2*G1,B2//G1 G1=Euclid(A2,B2)[0] return A1*A2, B1*B2 #②Garnerのアルゴリズム def Garner(Xlist,Ylist): #N≡X1 mod Y1≡X2 mod Y2≡ ... #X=[X1,X2, ...], Y=[Y1,Y2, ...] のとき、最小のNを求めよ。なければ-1を出力せよ。 if len(Xlist)!=len(Ylist): return -1 #解なしの判定、法の素因数振り分け for i in range(len(Ylist)): for j in range(i+1,len(Ylist)): G=Euclid(Ylist[i],Ylist[j])[0] if G>1 and (Xlist[i]-Xlist[j])%G!=0: return -1 #解なし elif G>1 and (Xlist[i]-Xlist[j])%G==0: Ylist[i],Ylist[j]=CRT_SnukeDistribute(Ylist[i],Ylist[j]) Xlist[i],Xlist[j]=Xlist[i]%Ylist[i],Xlist[j]%Ylist[j] #前から順にGarner Xg=Xlist[0] #初期条件: 漸化式のN≡Xi mod Yi まで満たすN値を格納する Yg=Ylist[0] for i in range(1,len(Xlist)): Xi,Yi=Xlist[i],Ylist[i] Xg=(Xg+Garnerv1(Xg,Yg,Xi,Yi)*Yg)%(Yg*Yi) Yg*=Yi #N=Xg+v1*Yg mod Yg*Yi return Xg def MODGarner(Xlist,Ylist,M): #N≡X1 mod Y1≡X2 mod Y2≡ ... とする。 #X=[X1,X2, ...], Y=[Y1,Y2, ...] のとき、最小のN mod Mを求めよ。 #適切なNがなかったり、Mが0ならば-1を出力せよ。 if len(Xlist)!=len(Ylist): return -1 MYihantei=False #そもそもYi=M を満たすならその値を記録。解なしなら出力。 MYi=0 for i in range(len(Ylist)): if Ylist[i]==M: MYihantei=True MYi=Xlist[i] NoAns=False for i in range(len(Ylist)): for j in range(i+1,len(Ylist)): G=Euclid(Ylist[i],Ylist[j])[0] if G>1 and (Xlist[i]-Xlist[j])%G!=0: NoAns=True elif G>1 and (Xlist[i]-Xlist[j])%G==0: Ylist[i],Ylist[j]=CRT_SnukeDistribute(Ylist[i],Ylist[j]) Xlist[i],Xlist[j]=Xlist[i]%Ylist[i],Xlist[j]%Ylist[j] if NoAns: #ここまでの判定を先に行う return -1 if MYihantei: return MYi #前から順にGarner。N=Xg[i]+v1*Yg[i] の式を念頭に。 #Xg[i] : i項のX mod Yi(順に計算), Yg[i] : Y0~Yi-1の累積積 mod Yi #ここで、Xg・Ygともに末尾は mod Mの値を格納することとする。 Ylist.append(M) #Xlist+1=Ylist となる Xg=[Xlist[0]]*(len(Ylist)) for i in range(len(Ylist)): Xg[i]%=Ylist[i] Yg=[1]*(len(Ylist)) for i in range(1,len(Xlist)): #N≡Xi mod Yi≡Xg[i]+v1*Yg mod Yi for j in range(i,len(Ylist)): #Yg[i]の更新 Yg[j]*=Ylist[i-1] Yg[j]%=Ylist[j] Xi,Yi=Xlist[i],Ylist[i] v1=(Xi-Xg[i])*EuclidMODInv(Yg[i],Yi)%Yi #v1*Yg[i]≡(Xi-Xg[i]) mod Yi for j in range(i,len(Ylist)): #N=X[i]+v1*Y[i] mod X Xg[j]+=v1*Yg[j] Xg[j]%=Ylist[j] return Xg[len(Ylist)-1] def Garnerv1(P,A,Q,B): #N≡P mod A≡Q mod B, N=P+v1*A (mod AB), 最適なv1を求めよ return (Q-P)*EuclidMODInv(A,B)%B #A*v1=(Q-P) mod B #③逆元計算 def Euclid(A,B): #Ax+By=gcd(A,B) を満たす gcd,x,yの組 #http://www.nct9.ne.jp/m_hiroi/light/pyalgo70.html から剽窃 Xs=(A,1,0) #A=A*1+B*0 Ys=(B,0,1) #B=A*0+B*1 while Ys[0]!=0: Q,Z=Xs[0]//Ys[0],Xs[0]%Ys[0] #A÷Bの商,余りを格納 Xs,Ys=Ys,(Z,Xs[1]-Q*Ys[1],Xs[2]-Q*Ys[2]) return Xs def EuclidMODInv(A,M): #A^(-1) mod M G,x,y=Euclid(A,M) if G!=1: return 0 else: return x%M def EulerPhi(N): #φ(N) if N<=0: return N CheckNumber=int(N) #素因数分解のライブラリから SoinsuList=[] #素因数分解の結果。(素数,次数)の形でtuple型に格納する for Soinsu in range(2,CheckNumber): if Soinsu*Soinsu>CheckNumber: break if CheckNumber%Soinsu!=0: continue SoinsuCount=0 while CheckNumber%Soinsu==0: SoinsuCount+=1 CheckNumber//=Soinsu SoinsuList.append((Soinsu,SoinsuCount)) if CheckNumber!=1: SoinsuList.append((CheckNumber,1)) EulerNo=int(N) for Prime,Order in SoinsuList: EulerNo=round(EulerNo*(1-(1/Prime))) return EulerNo def EulerMODInv(A,M): #A^(-1)≡A^(φ(M)-1) mod M if Euclid(A,M)[0]==1: return pow(A,EulerPhi(M)-1,M) else: return -1 #入力受取 A,B,a,b = map(int,input().split()) X = [a,b] Y = [A,B] N = Garner(X,Y) print(N)