結果
| 問題 |
No.186 中華風 (Easy)
|
| コンテスト | |
| ユーザー |
 navel_tos
|
| 提出日時 | 2023-02-28 13:30:43 |
| 言語 | Python3 (3.13.1 + numpy 2.2.1 + scipy 1.14.1) |
| 結果 |
AC
|
| 実行時間 | 35 ms / 2,000 ms |
| コード長 | 20,243 bytes |
| コンパイル時間 | 121 ms |
| コンパイル使用メモリ | 13,184 KB |
| 実行使用メモリ | 11,264 KB |
| 最終ジャッジ日時 | 2024-09-15 15:44:22 |
| 合計ジャッジ時間 | 1,998 ms |
|
ジャッジサーバーID （参考情報） |
judge3 / judge4 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 23 |
ソースコード
#yukicoder 中華風(Easy)
'''
yukicoder No.186
中華風(Easy)
https://yukicoder.me/problems/no/186
■問題設定
N≡X1 mod Y1≡X2 mod Y2≡X3 mod Y3 ･･･①
を満たす最小の正整数を求めよ。存在しなければ-1を出力せよ。
制約：0<=Xk,Yk<=10**6
■まえがき
中国剰余定理の問題。
逆元を使ったごり押し解法もあるが、今回はそれの練習と、
「法が互いに素でない場合への拡張」などについても考える。
■ユークリッドの互除法
N≡P mod A≡Q mod B
を満たすNが存在する必要十分条件は、
P≡Q mod gcd(A,B) ･･･②
が成立することである。ここで、gcdは最大公約数を表す。
AとBが互いに素である場合の実装例は
N≡P mod A≡Q mod B ･･･③
が与えられたとき、
N=Bx+Ay+ABk (kは整数) ･･･④
とすると、③④から
Bx≡P mod A
Ay≡Q mod B
ここで、AとBが互いに素であるなら、xはBの逆元のP倍である。yも同様なので
x≡P*(B**-1) mod A≡P*pow(B,A-2,A)
y≡Q*(A**-1) mod B≡Q*pow(A,B-2,B)
となる(逆元の計算はフェルマーの小定理とモジュラ逆数を用いる)。
これを④に代入して
N=B*P*pow(B,A-2,A)+A*Q*pow(A,B-2,B)+A*B*k (kは定数) ･･･⑤
とNの一般解が得られた。
この方法をGaussのアルゴリズムと呼び、③を反復することで①の解を求められる。
[問題点]
Gaussのアルゴリズムの実装上・仕様上、いくつかの問題が生じる。
(1)Aが素数でない場合
Bx≡1 mod A
と与えられた場合にBの逆元xをモジュラ逆数で計算できる条件は、
Aが素数であること が必要十分条件である。
(2)AとBが互いに素でない場合
AとBが互いに素でない場合、逆数判定以外にもトラブルがある。
具体的には、法の素因数の割り振りが必要。
(3)桁数のoverflow
⑤で得られた値は N≡X mod AB とおけるが、これを更に反復することで
N≡X' mod ABC(DEF...) と法の累積積が際限なく大きくなってしまう。
Pythonは多倍長整数(2^64以上の整数)が扱えるので気にしなくてよいが、
C++などの言語はoverflowしてしまう。
そうでなくてもモジュラ逆数の計算が大変。
■法が合成数の場合の逆数計算
問題点1を解消する。
(1)オイラーの定理
フェルマーの小定理を一般化したもので、オイラーのφ(ファイ)関数を用いる。
φ(n) : nと互いに素である1以上n以下の自然数の個数
pが素数なら φ(p)=p-1
p^kの中でpの倍数はp^(k-1)個あるので φ(p^k)=p^(k-1)*(p-1)=p^k*(1-(1/p))
mとnが互いに素なら φ(mn)=φ(m)*φ(n)
ここで、AとMが互いに素であれば、オイラー関数を用いて
A^(φ(M))≡1 mod M ･･･⑥
が成り立つ。上記のφ(M)の定義から、M=X1^Y1*X2^Y2* ... としたとき
φ(M)=φ(X1^Y1*X2^Y2* ...)
=φ(X1^Y1)*φ(X2^Y2)* ...
=M*(1-(1/X1))*(1-(1/X2))* ... ･･･⑦
となる。拡張ユークリッドの互除法よりも計算が遅いため、やや不利。
(2)拡張ユークリッドの互除法
Ax+By=1 あるいは Ax+By=gcd(A,B) となるx,yの組を見つけるアルゴリズム。
やっていることはAとBの剰余をAとBを用いて表す、をgcdに至るまでくり返すこと。
難しいので具体例を与える。
A,B=120,31 のとき、 Ax+By=1 となる x,yの組を求める。
ユークリッドの互除法でやるようにAとBを互いに割ってみる。
A,B=120,31 から順番にユークリッドの互除法を行うと
120=3*31+27 から 27=A-3*B
31=27*1+4 から 4=B-27= -A+ 4*B
27= 4*6+3 から 3=27-4=7*A-27*B
4= 3*1+1 から 1=4-3=-8*A+31*B (答)
となる。当例ではgcd(A,B)=1 だったが、そうでなければgcdに至るまで反復する。
さて、AとBが互いに素の場合、拡張ユークリッドの互除法では
Ax+By=1 となるx,yの組を求めることができた。するとAを法として
By≡1 mod A となる。これからyはBの逆元と求められる。
実装上は Ax+By=1 を変形することで
Axa+Bxb=gcd(A,B), Aya+Byb=gcd(A,B) (ただし xa,xb=1,0 ya,yb=0,1)
として、A/B=Q ･･･ Z となるZを求め順次更新を行えばよい。
[補足]
オイラーの定理の証明
φ(mn)=φ(m)φ(n) は、mn以下の正整数を mod mで分類し、各々でmnと互いに素なものを数える。
⑥は、M以下の自然数集合を集合X, Mと互いに素な自然数をAとする。
集合X={ x1, x2, ... , xφ(M)}
集合Xの各要素をA倍した集合Yを作成すると、法Mの元で集合Yと集合Xは合同である。
証明: 集合Xの任意の要素をxi, xjとする。
(1)集合Yのすべての要素は、集合Xの要素のいずれかと合致する。
(A・xiはMと互いに素のため、Axi≡xj mod M となるxjが存在する)
(2)集合Xだけでなく、集合Yのすべての要素は互いに異なる。
(Axi≡Axj mod M とすると、Aで除算して xi≡xj mod M となるが、
集合Xの定義上 xi≠xj かつ 0<xi,xj<M のため矛盾する)
(3)集合Xと集合Yの要素数は同一である。
集合Y={Ax1,Ax2, ... ,Axφ(M)}
≡{ x1, x2, ... , xφ(M)}≡集合X mod M
とおける。ここで集合Xの積和をZとして、集合Xと集合Yの積和をとると
Z≡A^(φ(M))*Z mod M → A^(φ(M))≡1 mod M (終)
■素因数の割り振り
問題点2を解消する。
(1)けんちょん式
互いに素でない自然数A,Bに対して、
N≡P mod A≡Q mod B ･･･⑧
を満たすNを求める。
AとBの最大公約数をG=math.gcd(A,B)とする。
まず、⑧が解を持つ必要十分条件を判定する。②を変形して
P-Q≡0 mod G ･･･⑨
を満たすか確認する。続いてAとBを素因数分解して、
A=X1^Y1a * X2^Y2a * * X4^Y4a * ...
B=X1^Y1b * * X3^Y3b * X4^Y4b * ...
と変形する。各素因数に対して、次数が多い方を残してA,Bを再構成すればよい。
具体例:
N≡3 mod 84≡15 mod 120 ･･･⑩ (答: 255)
G=math.gcd(84,120)=12
P-Q≡3-15≡-12≡0 mod G
A= 84=2^2*3^1* *7^1
B=120=2^3*3^1*5^1
次数が多いものを残してA,Bを再構成する。
3^1のみ共通項となっているが、好きな方に戻せばよい。次数上げは不要。
(補足: N≡X mod Y→N=Yk+X→N≡X+Tk mod TY となり次数上げは面倒)
N≡ 3 mod 21 : Aのうち、次数がB以上だった3,7を戻した
≡15 mod 40 : Bのうち、次数がA未満だった8,5を戻した
となる。オイラーの定理から逆元を得ると、φ(21)=12, φ(40)=16 なので
N≡P mod A≡Q mod B ･･･③'
x≡P*(B**-1) mod A≡P*pow(B,φ(A)-1,A)
y≡Q*(A**-1) mod B≡Q*pow(A,φ(B)-1,B)
から、A,B,P,Q=21,40,3,15 を代入すると
N=B*P*pow(B,φ(A)-1,A)+A*Q*pow(A,φ(B)-1,B)+A*B*k (kは定数) ･･･⑤'
=8655≡255 mod AB
として⑩の解を求められた。
当アルゴリズムでも十分高速だが、snuke式のほうが高速に求まる。
(2)snuke式
素因数分解を行わず、最大公約数のみを利用して割り振りを行う。
N≡P mod A≡Q mod B ･･･⑧
において、G=gcd(A,B)とする。GにはAとBに共通の素因数をすべて含む。
a1=A/G, b1=B/G として、AとBに共通する因数を削除する。
→a1とb1はAとBの固有の因数といえる。
A1=gcd(a1,G) とすると、G側に残されたa1の素因数のみを抽出できる。
→a1*A1 で、a1に含まれる素因数の次数を復元できる。
b1側には、Gからa1の素因数を「完全に」除去したものを振り分けたい。
b1側に振り分ける候補B1をまず求める。
B1=G/A1 として、gcd(A1,B1) を計算する。
G=A1*B1 である点に注意すると、振り分けが完全なら gcd(A1,B1)=1 となり、
逆に振り分けが不完全なら gcd(A1,B1)=G1 >1 となる。
以降、G1=gcd(A1,B1)=1 となるまで以下の「G1をA1に振る」操作を繰り返す。
A2=A1*G1, B2=B1/G1 と修正する。
G2=gcd(A2,B2) を計算し、G2=1 ならA1,B1,G1=A2,B2,G2 とする。
最後に A=a1*A1, B=b1*B1 で素因数の振り分けを完了する。
■Garnerのアルゴリズム
問題点3を解消するアルゴリズムである。
(1)「解を mod 10**9+7 で求めよ」のように任意modを求めたり、
(2)多倍長整数の保存に必要なメモリを削減したりできる。
X≡P mod A≡Q mod B≡R mod C(≡S mod D≡ ... )
のとき、X mod ABC(D...) の最小値を求める。
Garnerのアルゴリズムでは、求めるXを
X=v0 + v1*A + v2*AB (+ v3*ABC + ...) mod ABC(D...) ･･･(A)
として、Xの満たす条件式を順番に増やすようにXを加算してゆく。
具体例:
X≡4 mod 5≡1 mod 7≡2 mod 11 を満たす最小のXを求める。
(1)X≡4 mod 5 から、v0=4, X=4
(2)X≡1 mod 7 を満たすように、X≡4 mod 5を保つようXに5を加算してゆく。
X=v0+v1*5=4+v1*5≡1 mod 7 → 5*v1≡-3≡4 mod 7
5^(-1)≡3 mod 7 なので、v1≡12≡5 mod 7, X=29
(3)X≡2 mod 11 の条件を追加すると、
X=v0+v1*5+v2*5*7≡2 mod 11 → 29+35*v2≡(7+22)+(2+33)*v2≡7+2*v2≡2 mod 11
2^(-1)≡6 mod 11 なので、v2≡36≡3 mod 11
(4)以上から X=134 mod (5*7*11) (答)
最終的にX mod M (Mは素数)を出力する場合について。
法の累積積 ABC(D...) について、Mで割った余りが必要になる。
現時点までで X≡P mod A≡Q mod B≡R mod C を満たすXが既に求まっており、
X=v0 + v1*A + v2*AB=Xc (mod ABC) ･･･(B)
と計算できているものとする。これに 条件式 X≡S mod D を追加する状況を考える。
X=v0 + v1*A + v2*AB + v3*ABC
=Xc + v3*ABC ･･･(C)
(C)とX≡S mod D を連立させると
Xc + v3*ABC≡S mod D → v3*ABC≡(S-Xc) mod D ･･･(D)
(D)において、ABC^(-1) mod D を計算し、(C)に代入することでXが求まる。
mod M とする場合の計算だが、(C)について
Xc (mod ABC)≡Yc mod M ･･･(E)
とおく。条件式 X≡S mod D を追加すると
X=Xc + v3*ABC≡Yc + v3*ABC mod M ･･･(F)
となる。ABC(D...) mod M は予見して前計算しておけるため、
(D)で求めたv3を(F)に代入することで、X mod Mの答えが求まる。
■まとめ
X≡P mod A≡Q mod B≡R mod C(≡S mod D≡ ... )
のとき、X mod ABC(D...) の最小値を求めるGarnerのアルゴリズムを示す。
(1)法 A,B,C,(D, ...) が互いに素になるように素因数を振り分ける。
ここで、P≡Q mod gcd(A,B) が成立しなければ解なしである。
(2)Garnerのアルゴリズムを用いて、左側の条件式から順番にXを構成する。
X=v0 + v1*A + v2*AB (+ v3*ABC + ...) mod ABC(D...)
とXをおき、法の累積積の逆数を毎回計算してv0,v1,v2, ... と求める。
(3)整数の逆数は、オイラーの定理か拡張ユークリッドの互除法で求められる。
オイラーの定理は、φ(n) : nと互いに素である1以上n以下の自然数の個数 と定義して、
φ(M)=φ(p1^k1*p2^k2* ...)=round(M*(1-(1/p1))*(1-(1/p2))* ...)
A^(φ(M))≡1 mod M
となる性質を利用して、Mを素因数分解した後、Aの逆数を以下のように求める。
A^(-1)≡A^(φ(M)-1) mod M
拡張ユークリッドの互除法は、Ax+By=gcd(A,B) となるx,yの組を求める。
これから、AとMが互いに素の場合、
Ax+My=1 → Ax≡1 mod M
となるため、xがAの逆元として与えられる。
p=gcd(A,M)>1 の場合は逆数が存在しない。
(A=pa,M=pm とすると、 Ax=Mk+1 (kは定数) → pax=pmk+1 → 0≡0+1 mod p
となり矛盾する。Axは法M=pmの元、pの倍数しか取れないことから反証する)
(4)X mod Mを求める場合、あらかじめ法ABC(D...)の各累積積についてmod Mを前計算する。
左側の条件式から順にXa,Xb,Xc, ... と求め、毎回 mod Mを取って計算を進めればよい。
■実装
①素因数の振り分け
def CRT_SnukeDistribute(A,B):
N≡P mod A≡Q mod B とする。AとBが互いに素になるように適切に振り直せ。
②Garnerのアルゴリズム
def Garner(Xlist,Ylist):
N≡X1 mod Y1≡X2 mod Y2≡ ... とする。
X=[X1,X2, ...], Y=[Y1,Y2, ...] のとき、最小のNを求めよ。なければ-1を出力せよ。
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を出力せよ。
def Garnerv1(P,A,Q,B):
N≡P mod A≡Q mod B, N=P+v1*A (mod AB)とする。最適なv1を求めよ。
③逆元計算
def Euclid(A,B):
Ax+By=gcd(A,B) を満たすx,y,gcd(A,B)の組を1組出力せよ。
def EuclidMODInv(A,M):
AとMが互いに素なら、Ax+My=1 から Ax≡1 mod M を求め、Aの逆元x mod Mを出力せよ。
AとMが互いに素でなければ、0を出力せよ。
def EulerPhi(N):
φ(N)を出力せよ。
def EulerMODInv(A,M):
AとMが互いに素なら、A^(φ(M))=1 mod M を利用し、Aの逆元 A^(φ(M)-1) mod Mを出力せよ。
AとMが互いに素でなければ、0を出力せよ。
■参考文献
(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
#そもそもYi=M を満たすならその値を記録。解なしなら出力。
MYihantei=False
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
'''
ここから回答
'''
X=[]
Y=[]
for i in range(3):
x,y=list(map(int,input().split()))
X.append(x)
Y.append(y)
ans=Garner(X,Y)
if ans==0: #最小公倍数の計算
ans=Y[0]*Y[1]//Euclid(Y[0],Y[1])[0]
ans=ans*Y[2]//Euclid(ans,Y[2])[0]
print(ans)