import libnum
import gmpy2
A1=21856963452461630437348278434191434000066076750419027493852463513469865262064340836613831066602300959772632397773487317560339056658299954464169264467234407
B1=21856963452461630437348278434191434000066076750419027493852463513469865262064340836613831066602300959772632397773487317560339056658299954464169264467140596
#根据题意可定义一下C1
#C1=(B1!)%A1
#这里阶乘不可直接使用,所以将其拆分 B1比A1略小一点所以,可以
#C1=(B1!)%A1*(B1+1)%A /(B1+1)
#C1=(B1!)%A1*(B1+1)%A*(B1+2) /(B1+1)*(B1+2)
#C1=(B1!)%A1*(B1+1)%A*(B1+2)....*(A1-1) /(B1+1)*(B1+2)....(A1-1)
#分子直接追加至A1-1的阶乘,以方便使用威尔逊定理
#分母会多出若干项
#C1=(A1-1)!%A1 /(B1+1)*(B1+2)....(A1-1)
#分母用逆元来替换,
#C1=(A1-1)!%A1 * gmpy2.invert((B1+1),A1)*gmpy2.invert(B1+2),A1) ....*gmpy2((A1-1),A1)
#根据威尔逊定理,前一项替换掉,后面若干项保留
#C1=(-1)%A1 * gmpy2.invert((B1+1),A1)*gmpy2.invert(B1+2),A1) ....*gmpy2((A1-1),A1)
#使用for循环,计算出C1
C1=(-1)%A1
for i in range(B1,A1-1):
C1=C1*(gmpy2.invert((i+1),A1)%A1)
C1=C1%A1
print("C1=",C1)
A2=16466113115839228119767887899308820025749260933863446888224167169857612178664139545726340867406790754560227516013796269941438076818194617030304851858418927
B2=16466113115839228119767887899308820025749260933863446888224167169857612178664139545726340867406790754560227516013796269941438076818194617030304851858351026
C2=(-1)%A2
for i in range(B2,A2-1):
C2=C2*(gmpy2.invert((i+1),A2)%A2)
C2=C2%A2
print("C2=",C2)
#同理C2
n=85492663786275292159831603391083876175149354309327673008716627650718160585639723100793347534649628330416631255660901307533909900431413447524262332232659153047067908693481947121069070451562822417357656432171870951184673132554213690123308042697361969986360375060954702920656364144154145812838558365334172935931441424096270206140691814662318562696925767991937369782627908408239087358033165410020690152067715711112732252038588432896758405898709010342467882264362733
#由C1和C2 得到p 和q
p=gmpy2.next_prime(C1)
q=gmpy2.next_prime(C2)
#由pq 计算出r
r=n//(p*q)
print("r=",r)
e=0x1001
c=75700883021669577739329316795450706204502635802310731477156998834710820770245219468703245302009998932067080383977560299708060476222089630209972629755965140317526034680452483360917378812244365884527186056341888615564335560765053550155758362271622330017433403027261127561225585912484777829588501213961110690451987625502701331485141639684356427316905122995759825241133872734362716041819819948645662803292418802204430874521342108413623635150475963121220095236776428
#三个项的欧拉函数使用下面的方法。
phi=(p-1)*(q-1)*(r-1)
d=gmpy2.invert(e,phi)
m=gmpy2.powmod(c,d,n)
print(libnum.n2s(int(m)))