本篇对一些较为复杂的攻击进行研究

$e$、$phi$不互素攻击

题1：

之前已经了解过了模不互素，这里再引入一种不互素的情况：e和phi不互素。

e和phi不互素是指它俩的公约数不为1，下面是一道有关e和phi不互素的题：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import gmpy2
import libnum
import random
import uuid
flag="flag{"+str(uuid.uuid4())+"}"
m=libnum.s2n(flag)

while 1:
    e = random.randint(100,1000)
    p=libnum.generate_prime(1024)
    q=libnum.generate_prime(1024)
    phi_n=(p-1)*(q-1)
    t=gmpy2.gcd(e,phi_n)
    if  gmpy2.invert(e // t, phi_n) and t !=1:
        break
n=p*q
c=pow(m,e,n)
print("p=",p)
print("q=",q)
print("e=",e)
print("c=",c)

'''p= 127577058764408216374028752283743628765651360507566484643526093715329608267323381565274095814069864692746147152580906850350743742856555229701448239882612922698102985146366639955081466129923966803267071097174222576416224094182123529282235807472362341680183683025490897702891081336913842652559163341223338641607
q= 156492273708587234539506501480609692085997989594717058472605523051244522493701609615085173280972894139427194976925940854142835807192417391269823151398439665817176522629246535810290194301862945052149450578938260979300632842291287807430486629994530805358742405299538986591596966945727494262182814875780600646003
e= 750
c= 7029383721249299532521086933490698266831518266762255492452526410777276825803657150303837084263410309063739203644435184397762022380085273363900423091223180151147964276354189658062571415744140073426572149093499560918765793389358300893454490774387180728097370701432534877005948330689495694820361726719418371072834639369078180094444137972424909816959445043108154884587947573054460257114169961823509538355580857411319157089278918107229480661280354242839678709689654304688727345294473487201644985815128413154870914132135222144633969959773621933444285994038028721862094040876152694240238708737727034258171506516394913692187'''

从代码的if语句中我们知道了这里e和phi不互素，但是出题脚本中给了一点提示，我们可以引入一个中间量t作为e和phi的最大公约数，进而求出e的因子里和phi互质的数a，这样我们可以求出a和phi的逆元d。

这个逆元不同于之前的d，因为我们知道d=gmpy2.invert(e,phi)，而这里是把e作为m的指数而言的。但是现在我们把m^b看做整体，d1=gmpy2.invert(a,phi)，这时m由以前的m=pow(c,d,n)变成了m^b=pow(c,d1,n)

我们求出的pow(c,d,n)是m^b，所以我们想获得最后的m再开b次方就可以了，下面是根据上文分析写出的脚本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import gmpy2
from Crypto.Util.number import *

p= 127577058764408216374028752283743628765651360507566484643526093715329608267323381565274095814069864692746147152580906850350743742856555229701448239882612922698102985146366639955081466129923966803267071097174222576416224094182123529282235807472362341680183683025490897702891081336913842652559163341223338641607
q= 156492273708587234539506501480609692085997989594717058472605523051244522493701609615085173280972894139427194976925940854142835807192417391269823151398439665817176522629246535810290194301862945052149450578938260979300632842291287807430486629994530805358742405299538986591596966945727494262182814875780600646003
e= 750
c= 7029383721249299532521086933490698266831518266762255492452526410777276825803657150303837084263410309063739203644435184397762022380085273363900423091223180151147964276354189658062571415744140073426572149093499560918765793389358300893454490774387180728097370701432534877005948330689495694820361726719418371072834639369078180094444137972424909816959445043108154884587947573054460257114169961823509538355580857411319157089278918107229480661280354242839678709689654304688727345294473487201644985815128413154870914132135222144633969959773621933444285994038028721862094040876152694240238708737727034258171506516394913692187

n=p*q
phi=(p-1)*(q-1)
b=gmpy2.gcd(e,phi)
print(b)#b=6
a=e//b
d=gmpy2.invert(a,phi)
m_6=pow(c,d,n)
print(long_to_bytes(gmpy2.iroot(m_6,6)[0]))#b'flag{0e2f9add-31fd-4733-8f25-39297516f4e2}'

题2：

接下来加大难度，引入多组数据的e和phi不互素：

1
2
3
4
5
6
7
8
9
e1 = 14606334023791426
p1 = 121009772735460235364940622989433807619211926015494087453674747614331295040063679722422298286549493698150690694965106103822315378461970129912436074962111424616439032849788953648286506433464358834178903821069564798378666159882090757625817745990230736982709059859613843100974349380542982235135982530318438330859
q1 = 130968576816900149996914427770826228884925960001279609559095138835900329492765336419489982304805369724685145941218640504262821549441728192761733409684831633194346504685627189375724517070780334885673563409259345291959439026700006694655545512308390416859315892447092639503318475587220630455745460309886030186593
c1 = 11402389955595766056824801105373550411371729054679429421548608725777586555536302409478824585455648944737304660137306241012321255955693234304201530700362069004620531537922710568821152217381257446478619320278993539785699090234418603086426252498046106436360959622415398647198014716351359752734123844386459925553497427680448633869522591650121047156082228109421246662020164222925272078687550896012363926358633323439494967417041681357707006545728719651494384317497942177993032739778398001952201667284323691607312819796036779374423837576479275454953999865750584684592993292347483309178232523897058253412878901324740104919248

e2 = 13813369129257838
p2 = 121009772735460235364940622989433807619211926015494087453674747614331295040063679722422298286549493698150690694965106103822315378461970129912436074962111424616439032849788953648286506433464358834178903821069564798378666159882090757625817745990230736982709059859613843100974349380542982235135982530318438330859
q2 = 94582257784130735233174402362819395926641026753071039760251190444144495369829487705195913337502962816079184062352678128843179586054535283861793827497892600954650126991213176547276006780610945133603745974181504975165082485845571788686928859549252522952174376071500707863379238688200493621993937563296490615649
c2 = 7984888899827615209197324489527982755561403577403539988687419233579203660429542197972867526015619223510964699107198708420785278262082902359114040327940253582108364104049849773108799812000586446829979564395322118616382603675257162995702363051699403525169767736410365076696890117813211614468971386159587698853722658492385717150691206731593509168262529568464496911821756352254486299361607604338523750318977620039669792468240086472218586697386948479265417452517073901655900118259488507311321060895347770921790483894095085039802955700146474474606794444308825840221205073230671387989412399673375520605000270180367035526919

给了如上两组数，如果我们使用其中任意一组数求解m，则会在求解d的时候发现e没有phi的模逆，所以用老方法一定是行不通的。

按照第一题的思路，e和phi不互素，所以这两个数就应该有一个大于1的最大公因数，除了公因数，e的其他因数就和phi互素，我们可以让最大公因数为b，用a表示e//b作为和phi互质的数，这样就有了逆元d（注意这个d指代的意义）。

但是指依靠这一组m^14=pow(c,d1,n)无法求出m，尽管两组数里p是相等的，但是q1不等于q2，n是不相同的，m所在的有限域不同，第一组数解出的m不一定满足第二组，所以：

把两组数一起看，我们有：

c1^d1≡m^b mod n1
c2^d2≡m^b mod n2

p1和p2是相等的，根据同余的性质，我们可以把上面式子拆开：
$$
c_1^{d_1}= m^{b} \bmod p①\newline\ \
c_1^{d_1}= m^{b} \bmod q_1②\newline
$$

$$
c_2^{d_2}= m^{b} \bmod p③\newline\ \
c_2^{d_2}= m^{b} \bmod q_2④
$$

利用后两个模q的式子我们可以使用中国剩余定理来求m的特解（由于n1和n2不相等，m需要模n，此时m在两组数据中还无法求解，先作为特解）：

m_b=solve_crt( [c1d1,c2d2], [q1,q2] )

此时，我们可以把解出的m_b、q1和q2几个参数看做一个新的RSA，采用之前的思路，14可以分解成2和7，因为m^14=m^(2*7)，求7和phi（n=q1*q2）的逆元d，我们求出的pow(c,d,n)是m^2，所以我们想获得最后的m再开2次方就可以了

这里由于是模的q1*q2，对题目给的两组数都适合，这里解出的m才是我们真正要的m。

$$
m^2=（m^{14}）^{gmpy2.invert(7,phi)}\ \bmod (q1*q2)
$$

最后直接对m开根号就能得到flag，下面是根据以上推导写出的脚本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
import gmpy2
from libnum import *
from Crypto.Util.number import *
e1 = 14606334023791426
p1 = 121009772735460235364940622989433807619211926015494087453674747614331295040063679722422298286549493698150690694965106103822315378461970129912436074962111424616439032849788953648286506433464358834178903821069564798378666159882090757625817745990230736982709059859613843100974349380542982235135982530318438330859
q1 = 130968576816900149996914427770826228884925960001279609559095138835900329492765336419489982304805369724685145941218640504262821549441728192761733409684831633194346504685627189375724517070780334885673563409259345291959439026700006694655545512308390416859315892447092639503318475587220630455745460309886030186593
c1 = 11402389955595766056824801105373550411371729054679429421548608725777586555536302409478824585455648944737304660137306241012321255955693234304201530700362069004620531537922710568821152217381257446478619320278993539785699090234418603086426252498046106436360959622415398647198014716351359752734123844386459925553497427680448633869522591650121047156082228109421246662020164222925272078687550896012363926358633323439494967417041681357707006545728719651494384317497942177993032739778398001952201667284323691607312819796036779374423837576479275454953999865750584684592993292347483309178232523897058253412878901324740104919248
n1=p1*q1

e2 = 13813369129257838
p2 = 121009772735460235364940622989433807619211926015494087453674747614331295040063679722422298286549493698150690694965106103822315378461970129912436074962111424616439032849788953648286506433464358834178903821069564798378666159882090757625817745990230736982709059859613843100974349380542982235135982530318438330859
q2 = 94582257784130735233174402362819395926641026753071039760251190444144495369829487705195913337502962816079184062352678128843179586054535283861793827497892600954650126991213176547276006780610945133603745974181504975165082485845571788686928859549252522952174376071500707863379238688200493621993937563296490615649
c2 = 7984888899827615209197324489527982755561403577403539988687419233579203660429542197972867526015619223510964699107198708420785278262082902359114040327940253582108364104049849773108799812000586446829979564395322118616382603675257162995702363051699403525169767736410365076696890117813211614468971386159587698853722658492385717150691206731593509168262529568464496911821756352254486299361607604338523750318977620039669792468240086472218586697386948479265417452517073901655900118259488507311321060895347770921790483894095085039802955700146474474606794444308825840221205073230671387989412399673375520605000270180367035526919
n2=p2*q2

phi1=(p1-1)*(q1-1)
phi2=(p2-1)*(q2-1)
b1=gmpy2.gcd(e1,phi1)
b2=gmpy2.gcd(e2,phi2)
a1=e1//b1
a2=e2//b2
d1=gmpy2.invert(a1,phi1)
d2=gmpy2.invert(a2,phi2)
m_b1=pow(c1,d1,n1)
m_b2=pow(c2,d2,n2)
cd1=m_b1%q1
cd2=m_b2%q2
m_b=solve_crt([cd1,cd2],[q1,q2])
phi=(q1-1)*(q2-1)
b=gmpy2.gcd(b1,phi)
a=b1//b
m_2=pow(m_b,gmpy2.invert(a,phi),q1*q2)
print(long_to_bytes(gmpy2.iroot(m_2,2)[0]))#b"flag{gcd_e&\xcf\x86_isn't_1}"

维纳攻击

维纳攻击（Wiener attack）是一种依靠连分数（continued fraction）的求解方式。适用于非常接近某一个值（比如 1）时，求一个比例关系。

适用于维纳攻击的条件：e很大，且有
$$
d<\frac{1}{3}N^{\frac{1}{4}}\ , \ q<p<2q
$$
为了更好的理解维纳攻击，我们需要先了解连分数的有关知识。

密码学数学基础 - 连分数理论

连分数

分数转连分数：

我们可以在python里将分数转换为连分数：

注意到有
$$
\frac{p}{q}=[\frac{p}{q}]+\frac{p \bmod q}{q}=[\frac{p}{q}]+\frac{1}{\frac{q}{p\ \bmod\ q}}
$$
我们需要继续转化分母，得到p,q=q,p%q

1
2
3
4
5
6
7
8
def f2c(p,q):
  a=[]
  while q:#q=0时跳出循环
    a.append(p//q)#求an
    p,q=q,p%q
  return a
  
print(f2c(43,19))#[2,3,1,4]

连分数转分数：

1
2
3
4
5
6
7
8
def c2f(a):#连分数转分数
  x,y=1,0#x为分子y为分母
  for a in a[::-1]:#逆序累加
    x,y=a*x+y,x
  return (x,y)

a=[2,3,1,4]
print(c2f(a))#[43,19]

维纳攻击原理

利用公式推导：
$$
ed \bmod phi = 1
$$

$$
ed=k\cdot phi+1 \ \ ===>\ \ ed\approx kphi
$$

当n较大：phi ≈ n，可以推得：
$$
ed\approx kn\ \ ,\ \ \frac{e}{n}\approx \frac{k}{d}
$$
此时，我们可以从e/n各阶段的连分数中找到能满足条件的一个连分数。

具体证明尚未研究，先记着结论。

我们可以将这个过程得到的所有分数爆破出来：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def c2f(a):#连分数转分数
  x,y=1,0#x为分子y为分母
  for a in a[::-1]:#逆序累加
    x,y=a*x+y,x
  return (x,y)

def sf(a):
  fs=[]
  for i in range(1,len(a)):
    fs.append(c2f(a[:i]))
  return fs

a=[2,3,1,4]
print(sf(a))
#[(1, 2), (3, 7), (4, 9)],我们获得了2,7/3,9/4三个渐进分数对应的分数，注意分母分子顺序，前面是分母，后面是分子。

我们通过找到真正的d和k，进而求phi，破解RSA，这就是维纳攻击。
$$
\frac{e}{n}\approx \frac{k}{d}\ \ ,\ \ ed=k\cdot phi+1
$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def c2f(a):#连分数转分数
  x,y=1,0#x为分子y为分母
  for a in a[::-1]:#逆序累加
    x,y=a*x+y,x
  return (x,y)

def sf(a):
  f=[]
  for i in range(1,len(a)):
    f.append(c2f(a[:i]))
  return f

def wiener_attack(e,n,d):
  q=f2c(e,n)
  for (D,K) in sf(q):
    if D==d and K!=0:
      return (e*d-1)//K  #phi=(ed-1)//k

但是问题是我们并不知道d，然而我们知道n和能爆破出的phi，所以我们可以求p和q：

1
2
3
4
5
6
7
8
9
import gmpy2

def solve_p_q(n,phi):
	x=n-phi+1  #x=p+q=n-p*q+(p+q)-1+1
    y=x*x-4*n  #y=(p-q)^2=(p+q)^2-4n
    z=gmpy2.iroot(y,2)[0]	#z=q-p
    q=(x+z)//2
    p=(x-z)//2
    return p,q

所以现在路线变了，我们可以逐个试d并得到一个phi，进而求pq，如果pq相乘还能得到n，我们就能得到真正的p,q，进而完成维纳攻击。

之前的程序稍作修改并整合一下，一个脚本就写好了：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
import gmpy2
from Crypto.Util.number import *

def f2c(p,q):
 a=[]
 while q:#q=0时跳出循环
  a.append(p//q)#求an
  p,q=q,p%q
 return a

def c2f(a):#连分数转分数
  x,y=1,0#x为分子y为分母
  for a in a[::-1]:#逆序累加
    x,y=a*x+y,x
  return (x,y)

def sf(a):
 f=[]
 for i in range(1,len(a)):
  f.append(c2f(a[:i]))
 return f

def solve_p_q(n,phi):
  x=n-phi+1  #x=p+q=n-p*q+(p+q)-1+1
  y=x*x-4*n  #y=(p-q)^2=(p+q)^2-4n
  z=gmpy2.iroot(y,2)[0]  #z=q-p
  q=(x+z)//2
  p=(x-z)//2
  return p,q

def wiener_attack(e,n):
 a=f2c(e,n)
 for (d,k) in sf(a):
  if d!=0 and k!=0:
    phi=(e*d-1)//k
    (p,q)=solve_p_q(n,phi)
    if p*q==n:
      m=pow(c,d,n)
      print(long_to_bytes(m))
      break#这里break一方面降低运算量，另一方面随着phi增大，x减小，y可能会成负数，z不能开方报错，或者直接return p,q在外面计算也可以

n=10117654819914858266933329374955416887632623769133893090370644563766857602175135282123557069130319485164837923640109707980173187311717714731455048711732890650502393864146567993461691536083330111489342144526034765633893707639465391971659424271400604010051260552831236934617979897198594708643604050329358522040572553492557329327193918343289526476096522685686128709365882540965089876020772451339243051387630968483513100881164719486794479191020450727996212211201807165531274853014517030221518336293845148545671493267736094904720639061287350709209363413742534108909427583009442175135992281621755221466312230819838164838443
e=9821176723459156737162590528498355378679103255669217165700920581299681706733929195884953659518540987340485400582895899813962495604183457377106880733695051072483080763292039235986138262683331839202120494074112671026731661652894069539798773005571447249078493499067926710777981456836165288713067372341722891618455469381820299718375899142104630769808052209736800306823537530231432735329122809506084509365041929994661643608897946882821172042151091436805833261237973879388223150132281413026451714557979869417700194664385291198650864896408143681530963859508908402067374010247738488460155501935400209160082801993945408813513
c=838279327100183842450828959704405383407020841408916285706333834213457887909003957632210005175559669378601653437817015283864372967567045814324446631403131762020243676699045950634510503063630325940620012467503745448306616066932850616255337850567483377961974904557893440882626053258665407295455129124436515237709284339335286446533668177967589716697626618148973094630870394728363810896842456940427809475399274556698585866672673971202275767545143765482579343055060966723452080734906835537296838390574697520016738791840483248208135607782781572593502322902003706653541803004568389346187087997006034664608287414331955367370

wiener_attack(e,n)
#b"SKSEC{Do_y0u_Kn0w_Wi3n3r's_4ttack}"

三素数维纳攻击

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def wiener(e, n):
    m = 12345
    c = pow(m, e, n)
    q0 = 1
    list1 = continued_fraction(Integer(e) / Integer(n))
    conv = list1.convergents()
    for i in conv:
        k = i.numerator()
        q1 = i.denominator()
        for r in range(20):
            for s in range(20):
                d = r * q1 + s * q0
                m1 = pow(c, d, n)
                if m1 == m:
                    return d
        q0 = q1

n = 1051380872316823433883636764584789285333736434524933570584976014001214015195582230652890036553897972507709868668015209942849015953126067065245745528656185713975834472569979390688522786334308265267342203713851528468719085926942222541344057327637925939220941677328780957242592769934672272901444344368132925072596706265951382183962516395989020286284674233929275036820658314907868575906130423049124591076625494555958927499995127864266451094815485831694951169813864579090722449382301556276181043387840715294133205479396272801536638003449150191833236756230648343641542001617352347728075031395009393805016183916046647042759712436405458561551257959973189661621088021485676031370872948254252297718720969276447566410825701049555286286171854286210170310303122812022424180992567945938220732070586725895910160295026510702457152084383868466740490026825254620849716745141621210206706738899418272128830937154184366826368016614460538955646609
# p*r*q
e = 401736266490324552517217032225246450087649254466718075727966287757972026925163477893539047721223580215545684494442007744774900721388024851291175127388784748482162997733029569314385345465782341732251402815344448962923207292061517106665954257137567332560956608813590130048988674456360832136718517222056109287921167927946304453828334895683526480626294649396329610693963520850896266343295069273447109718037817671775114733403154525393468056590686142168244467303584451585069101584798590889529292526446832857035880715351476044875846461493125136227486230514153648388383654116769560511713813660841282782474519612727725715364718394717990336084691567425849567195154254431858216101902309862385524999911892987920197370958210721993073768994565808825036701082088990540391431078200077738350023726344360396487614841476654096629864029557926380538230445278881935443675347998379733770494280427033480933837349115145552011734379818969027809471037
d=wiener(e, n)
print(d)

威尔逊定理相关攻击

下面是威尔逊定理公式：

当且仅当p为素数时。有：
$$
(p-1)!\equiv -1\ (\bmod p)
$$

证明

知道：
$$
p-1\equiv -1(\bmod p)
$$
需要证明$(p-2)!\equiv 1(\bmod p)$，同时，去掉第一个因子$1$，得到剩下的$2,3,…,p-2$，$p$为质数肯定是奇数，所以$p-2$也是奇数，剩下的乘数有偶数个。

对于二次剩余$x^2\equiv 1(\bmod p)$，两个解为$1,p-1$，所以除了这两个数以外的其他数不存在逆元为其本身。除去$1,p-1$，在剩余类中其余数两两配对互为逆元。

题目：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import sympy
import random

def myGetPrime():
    A= getPrime(513)
    print(A)
    B=A-random.randint(1e3,1e5)
    print(B)
    return sympy.nextPrime((B!)%A)
p=myGetPrime()
#A1=21856963452461630437348278434191434000066076750419027493852463513469865262064340836613831066602300959772632397773487317560339056658299954464169264467234407
#B1=21856963452461630437348278434191434000066076750419027493852463513469865262064340836613831066602300959772632397773487317560339056658299954464169264467140596

q=myGetPrime()
#A2=16466113115839228119767887899308820025749260933863446888224167169857612178664139545726340867406790754560227516013796269941438076818194617030304851858418927
#B2=16466113115839228119767887899308820025749260933863446888224167169857612178664139545726340867406790754560227516013796269941438076818194617030304851858351026

r=myGetPrime()

n=p*q*r
#n=85492663786275292159831603391083876175149354309327673008716627650718160585639723100793347534649628330416631255660901307533909900431413447524262332232659153047067908693481947121069070451562822417357656432171870951184673132554213690123308042697361969986360375060954702920656364144154145812838558365334172935931441424096270206140691814662318562696925767991937369782627908408239087358033165410020690152067715711112732252038588432896758405898709010342467882264362733
c=pow(flag,e,n)
#e=0x1001
#c=75700883021669577739329316795450706204502635802310731477156998834710820770245219468703245302009998932067080383977560299708060476222089630209972629755965140317526034680452483360917378812244365884527186056341888615564335560765053550155758362271622330017433403027261127561225585912484777829588501213961110690451987625502701331485141639684356427316905122995759825241133872734362716041819819948645662803292418802204430874521342108413623635150475963121220095236776428
#so,what is the flag?

本题中我们需要注意这个地方：

sympy.nextPrime((B!)%A)

这里是生成了B的阶乘又取余A后的下一个质数，但是B是一个很大的数，用这么大的数计算阶乘基本不可能。所以这里我们可以做一些数学运算：

既然题目里用了p和q两个质数，我们先单看p：

p=sympy.nextPrime((B1!)%A1)

由威尔逊定理可以知道：
$$
(A_1-1)! \equiv -1\ (\bmod A_1)
$$
所以我们可以引入一个k，让：
$$
B_1!\cdot k\equiv -1\ (\bmod A_1)\ ,\ k=\frac{(A_1-1)!}{B_1!}
$$
这里k约分一下就是B1之后到A1-1的数的积，由于B和A就差一个随机在483和485（注意16进制）的值，所以求解k也就是几百个大数的乘积，这里由于我们接下来要对k取模逆，所以k后面取模A1。

两边同时乘k的逆元：
$$
B_1!\equiv-k^{-1}\ (\bmod A_1)
$$
这样就能得到了B1的阶乘，进而就能求解p和q了。

写个简单的代码用来求解：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import gmpy2
from Crypto.Util.number import *

def wei(A,B):
  k=1
  for i in range(B+1,A):
   k=(k*i)%A#模运算定律：(a*b)%p=(a%p*b%p)%p=((a%p)*(b%p))%p进一步简化运算
  K=gmpy2.invert(k,A)
  return gmpy2.next_prime(-K%A)

A1=21856963452461630437348278434191434000066076750419027493852463513469865262064340836613831066602300959772632397773487317560339056658299954464169264467234407
B1=21856963452461630437348278434191434000066076750419027493852463513469865262064340836613831066602300959772632397773487317560339056658299954464169264467140596

A2=16466113115839228119767887899308820025749260933863446888224167169857612178664139545726340867406790754560227516013796269941438076818194617030304851858418927
B2=16466113115839228119767887899308820025749260933863446888224167169857612178664139545726340867406790754560227516013796269941438076818194617030304851858351026

c=75700883021669577739329316795450706204502635802310731477156998834710820770245219468703245302009998932067080383977560299708060476222089630209972629755965140317526034680452483360917378812244365884527186056341888615564335560765053550155758362271622330017433403027261127561225585912484777829588501213961110690451987625502701331485141639684356427316905122995759825241133872734362716041819819948645662803292418802204430874521342108413623635150475963121220095236776428
n=85492663786275292159831603391083876175149354309327673008716627650718160585639723100793347534649628330416631255660901307533909900431413447524262332232659153047067908693481947121069070451562822417357656432171870951184673132554213690123308042697361969986360375060954702920656364144154145812838558365334172935931441424096270206140691814662318562696925767991937369782627908408239087358033165410020690152067715711112732252038588432896758405898709010342467882264362733
e=0x1001

p=wei(A1,B1)
q=wei(A2,B2)
r=n//q//p
phi=(p-1)*(q-1)*(r-1)#注意多素数问题
d=gmpy2.invert(e,phi)
m=pow(c,d,n)
print(long_to_bytes(m))

Pollard’s p-1攻击

光滑数是可以分解为小素数乘积的正整数。

如果$p$是$n$的因子，$p-1$是光滑数，我们对$p-1$进行质因数分解，最大的因子$d^k$就是$B$，而这个质数就是$B-powersmooth$。
例如我们分解$3696=2^4\cdot3\cdot7\cdot11$，那么$3697$就是$16-powersmooth$。

如果$p$是$n$的因子，$p-1$是光滑数，我们可以使用Pollard’s p-1算法分解$n$。

Pollard’s p-1算法主要内容：如果$p-1$正好是一些很小的质数的乘积，那么$p-1$就能整除$n!$。

可以设$p-1=p_1p_2…p_n(\forall 1≤i≤n,p_i≤B)$

如果$p_1,p_2,…,p_n$各不相同，那么得到：
$$
p_1p_2…p_n|B!\Rightarrow (p-1)|B!\Rightarrow B!=k(p-1)
$$
根据上面的公式，利用费马小定理：由于$a^{B!}\equiv a^{k(p-1)}\equiv 1(\bmod p)$，我们计算$gcd(a^{B!}-1,n)$，如果结果不是$1$或者$n$，那就分解出了$p$。

费马小定理中，参数$a$没什么限制，可以先把参数赋值$a=2$。然后开始从小于$B$的范围（含有$p_1,p_2,…,p_n$）里去遍历筛选，大概率能成功分解。

注意，实际RSA中我们并不知道$B$究竟是多少，但是我们很大可能能根据题目的预设知道一个界，我们在这里可以把这个界暂定为$B$。

为了优化算法，我们不需要每一次计算都要计算最大公约数，因为我们是要寻找在$B$这个界最大的那个因子，所以实际情况上不需要每次都算一次$gcd(a^x-1,n)$，最大的素数大概率要在遍历的后半程，所以我们可以在比如$<0.8B$的范围内不进行计算公因数，节省时间。

同时，由于$gcd(a^{B!}-1,n)=gcd(a^{B!}-1 \bmod n,n)$，把数减小，可以进一步优化算法。

实际写代码进行计算时，我们可以采取降幂运算，相当于每次给$a$重新赋值：
$$
a^{b!}\bmod n=(a^{(b-1)!}\bmod n)\bmod n
$$
下面题目为NCTF2019的childRSA：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
from random import choice
from Crypto.Util.number import isPrime, sieve_base as primes
from flag import flag

def getPrime(bits):
while True:
  n = 2
  while n.bit_length() < bits:
      n *= choice(primes)	#primes为前10000个素数的列表
  if isPrime(n + 1):
      return n + 1

e = 0x10001
m = int.from_bytes(flag.encode(), 'big')
p, q = [getPrime(2048) for _ in range(2)]
n = p * q
c = pow(m, e, n)

n = 32849718197337581823002243717057659218502519004386996660885100592872201948834155543125924395614928962750579667346279456710633774501407292473006312537723894221717638059058796679686953564471994009285384798450493756900459225040360430847240975678450171551048783818642467506711424027848778367427338647282428667393241157151675410661015044633282064056800913282016363415202171926089293431012379261585078566301060173689328363696699811123592090204578098276704877408688525618732848817623879899628629300385790344366046641825507767709276622692835393219811283244303899850483748651722336996164724553364097066493953127153066970594638491950199605713033004684970381605908909693802373826516622872100822213645899846325022476318425889580091613323747640467299866189070780620292627043349618839126919699862580579994887507733838561768581933029077488033326056066378869170169389819542928899483936705521710423905128732013121538495096959944889076705471928490092476616709838980562233255542325528398956185421193665359897664110835645928646616337700617883946369110702443135980068553511927115723157704586595844927607636003501038871748639417378062348085980873502535098755568810971926925447913858894180171498580131088992227637341857123607600275137768132347158657063692388249513
c = 26308018356739853895382240109968894175166731283702927002165268998773708335216338997058314157717147131083296551313334042509806229853341488461087009955203854253313827608275460592785607739091992591431080342664081962030557042784864074533380701014585315663218783130162376176094773010478159362434331787279303302718098735574605469803801873109982473258207444342330633191849040553550708886593340770753064322410889048135425025715982196600650740987076486540674090923181664281515197679745907830107684777248532278645343716263686014941081417914622724906314960249945105011301731247324601620886782967217339340393853616450077105125391982689986178342417223392217085276465471102737594719932347242482670320801063191869471318313514407997326350065187904154229557706351355052446027159972546737213451422978211055778164578782156428466626894026103053360431281644645515155471301826844754338802352846095293421718249819728205538534652212984831283642472071669494851823123552827380737798609829706225744376667082534026874483482483127491533474306552210039386256062116345785870668331513725792053302188276682550672663353937781055621860101624242216671635824311412793495965628876036344731733142759495348248970313655381407241457118743532311394697763283681852908564387282605279108

查表可以知道第10000个素数为104729，所有的$p-1$的因子都是从前一万个素数产生的，所以这就可以使用上述算法进行攻击：（实测跑到$i=100000$大概两三分钟）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
from Crypto.Util.number import *
import gmpy2

def pollard(n,B):
    i=2
    a=2
    limit=int(0.8*B)
    for i in range(2,limit):
        a=pow(a,i,n)
    for i in range(limit,B+1):
        a=pow(a,i,n)
        p=gmpy2.gcd(a-1,n)
        if 1<p<n:
            return p,n//p

n=32849718197337581823002243717057659218502519004386996660885100592872201948834155543125924395614928962750579667346279456710633774501407292473006312537723894221717638059058796679686953564471994009285384798450493756900459225040360430847240975678450171551048783818642467506711424027848778367427338647282428667393241157151675410661015044633282064056800913282016363415202171926089293431012379261585078566301060173689328363696699811123592090204578098276704877408688525618732848817623879899628629300385790344366046641825507767709276622692835393219811283244303899850483748651722336996164724553364097066493953127153066970594638491950199605713033004684970381605908909693802373826516622872100822213645899846325022476318425889580091613323747640467299866189070780620292627043349618839126919699862580579994887507733838561768581933029077488033326056066378869170169389819542928899483936705521710423905128732013121538495096959944889076705471928490092476616709838980562233255542325528398956185421193665359897664110835645928646616337700617883946369110702443135980068553511927115723157704586595844927607636003501038871748639417378062348085980873502535098755568810971926925447913858894180171498580131088992227637341857123607600275137768132347158657063692388249513
B=104729
e=0x10001
c = 26308018356739853895382240109968894175166731283702927002165268998773708335216338997058314157717147131083296551313334042509806229853341488461087009955203854253313827608275460592785607739091992591431080342664081962030557042784864074533380701014585315663218783130162376176094773010478159362434331787279303302718098735574605469803801873109982473258207444342330633191849040553550708886593340770753064322410889048135425025715982196600650740987076486540674090923181664281515197679745907830107684777248532278645343716263686014941081417914622724906314960249945105011301731247324601620886782967217339340393853616450077105125391982689986178342417223392217085276465471102737594719932347242482670320801063191869471318313514407997326350065187904154229557706351355052446027159972546737213451422978211055778164578782156428466626894026103053360431281644645515155471301826844754338802352846095293421718249819728205538534652212984831283642472071669494851823123552827380737798609829706225744376667082534026874483482483127491533474306552210039386256062116345785870668331513725792053302188276682550672663353937781055621860101624242216671635824311412793495965628876036344731733142759495348248970313655381407241457118743532311394697763283681852908564387282605279108
p,q=pollard(n,B)
phi=(p-1)*(q-1)
d=gmpy2.invert(e,phi)
m=pow(c,d,n)
print(long_to_bytes(m))#b'NCTF{Th3r3_ar3_1ns3cure_RSA_m0duli_7hat_at_f1rst_gl4nce_appe4r_t0_be_s3cur3}'

不过，据说还有一种方法比这种快很多：

之前我们知道了$a^{B!}\equiv a^{k(p-1)}\equiv 1(\bmod p)$，由于$p-1$是由前一万个素数中的若干个求得，那么我们设前一万个素数乘积为$X$，那么这个$X$就是$p-1$的倍数，让$X=k(p-1)$

根据费马小定理，$a^X-1=k_1p$

由于$gcd(a^X-1,n)=p$，这样就能很快得到$p$。

注意，当$a=2$时，$2^X$也是非常大的，我们需要优化一下算法：

$2^X=1+k_1p$，还有$2^X \bmod n=2^X-k_2n=2^X-k_2pq$

对其进行模$p$，得到$2^X\bmod n \equiv 2^X(\bmod p)$

得到$2^X \bmod n=2^X (\bmod p)$，就有$2^X\bmod n\equiv 1(\bmod p)$

下面是根据以上分析写的代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
from Crypto.Util.number import long_to_bytes,sieve_base#前一万个素数数组
import gmpy2

n=32849718197337581823002243717057659218502519004386996660885100592872201948834155543125924395614928962750579667346279456710633774501407292473006312537723894221717638059058796679686953564471994009285384798450493756900459225040360430847240975678450171551048783818642467506711424027848778367427338647282428667393241157151675410661015044633282064056800913282016363415202171926089293431012379261585078566301060173689328363696699811123592090204578098276704877408688525618732848817623879899628629300385790344366046641825507767709276622692835393219811283244303899850483748651722336996164724553364097066493953127153066970594638491950199605713033004684970381605908909693802373826516622872100822213645899846325022476318425889580091613323747640467299866189070780620292627043349618839126919699862580579994887507733838561768581933029077488033326056066378869170169389819542928899483936705521710423905128732013121538495096959944889076705471928490092476616709838980562233255542325528398956185421193665359897664110835645928646616337700617883946369110702443135980068553511927115723157704586595844927607636003501038871748639417378062348085980873502535098755568810971926925447913858894180171498580131088992227637341857123607600275137768132347158657063692388249513
c = 26308018356739853895382240109968894175166731283702927002165268998773708335216338997058314157717147131083296551313334042509806229853341488461087009955203854253313827608275460592785607739091992591431080342664081962030557042784864074533380701014585315663218783130162376176094773010478159362434331787279303302718098735574605469803801873109982473258207444342330633191849040553550708886593340770753064322410889048135425025715982196600650740987076486540674090923181664281515197679745907830107684777248532278645343716263686014941081417914622724906314960249945105011301731247324601620886782967217339340393853616450077105125391982689986178342417223392217085276465471102737594719932347242482670320801063191869471318313514407997326350065187904154229557706351355052446027159972546737213451422978211055778164578782156428466626894026103053360431281644645515155471301826844754338802352846095293421718249819728205538534652212984831283642472071669494851823123552827380737798609829706225744376667082534026874483482483127491533474306552210039386256062116345785870668331513725792053302188276682550672663353937781055621860101624242216671635824311412793495965628876036344731733142759495348248970313655381407241457118743532311394697763283681852908564387282605279108
e=0x10001

X=1
for i in sieve_base:
    X*=i
p=gmpy2.gcd(pow(2,X,n)-1,n)
q=n//p
phi=(p-1)*(q-1)
d=gmpy2.invert(e,phi)
m=pow(c,d,n)
print(long_to_bytes(m))#b'NCTF{Th3r3_ar3_1ns3cure_RSA_m0duli_7hat_at_f1rst_gl4nce_appe4r_t0_be_s3cur3}'

AMM攻击

原理：
$$
m^e\equiv c(\bmod p)
$$
p在已知e的情况下易分解得到s，最常见的情况就是e是p-1的因子：
$$
p-1=e^ts
$$

具体的计算算法先不记了，有点复杂。这里简单找了几个脚本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
from gmpy2 import *
from Crypto.Util.number import *
import random
import math

def onemod(e, q):
    p = random.randint(1, q-1)
    while(powmod(p, (q-1)//e, q) == 1):  # (r,s)=1
        p = random.randint(1, q)
    return p

def AMM_rth(o, r, q):  # r|(q-1
    assert((q-1) % r == 0)
    p = onemod(r, q)

    t = 0
    s = q-1
    while(s % r == 0):
        s = s//r
        t += 1
    k = 1
    while((s*k+1) % r != 0):
        k += 1
    alp = (s*k+1)//r

    a = powmod(p, r**(t-1)*s, q)
    b = powmod(o, r*a-1, q)
    c = powmod(p, s, q)
    h = 1

    for i in range(1, t-1):
        d = powmod(int(b), r**(t-1-i), q)
        if d == 1:
            j = 0
        else:
            j = (-math.log(d, a)) % r
        b = (b*(c**(r*j))) % q
        h = (h*c**j) % q
        c = (c*r) % q
    result = (powmod(o, alp, q)*h)
    return result

def ALL_Solution(m, q, rt, cq, e):
    mp = []
    for pr in rt:
        r = (pr*m) % q
        # assert(pow(r, e, q) == cq)
        mp.append(r)
    return mp


def calc(mp, mq, e, p, q):
    i = 1
    j = 1
    t1 = invert(q, p)
    t2 = invert(p, q)
    for mp1 in mp:
        for mq1 in mq:
            j += 1
            if j % 1000000 == 0:
                print(j)
            ans = (mp1*t1*q+mq1*t2*p) % (p*q)
            if check(ans):
                return
    return


def check(m):
    try:
        a = long_to_bytes(m).decode('utf-8')
        if 'flag' in a:
            print(a)
            return True
        else:
            return False
    except:
        return False


def ALL_ROOT2(r, q):  # use function set() and .add() ensure that the generated elements are not repeated
    li = set()
    while(len(li) < r):
        p = powmod(random.randint(1, q-1), (q-1)//r, q)
        li.add(p)
    return li


if __name__ == '__main__':
    c=257142110273808470353932790260834385423195836240947804445043168231453970823430490856329111213189231868252755966778005793663199810057119020960651706963993282769973864494156391386075051462934130430940835471240579618826284534276211239463295726838485842840227094026972973568755487969486042883286010886114853697684439422818236539054061509350796956754115566635281319223039237365252887952785342564138484884780518096718205122966618504897723001949304358011007666606602453554821466955282438373689460
    p=19936217237033573753894155829225442996983351260648079359790039921271372502978979089255944115115517326597713075963827804905342890524775515960525156828602110668810581490760693054089408458756578501439658313103640176379107393837271103113860091147479 
    q=15835321764158530652228056326151107407225353281494241384475003369802102321837994236516023188966354559804738234055721039001493117376126537075924755525606506532010338074320834686510734601985723266842927770379824450335028680072267894221203184802411 
    e=14
    cp = c % p
    cq = c % q

    mp = AMM_rth(cp, e, p)
    mq = AMM_rth(cq, e, q)

    rt1 = ALL_ROOT2(e, p)
    rt2 = ALL_ROOT2(e, q)

    amp = ALL_Solution(mp, p, rt1, cp, e)
    amq = ALL_Solution(mq, q, rt2, cq, e)

    calc(amp, amq, e, p, q)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
import random
import math
import libnum
import time
from Crypto.Util.number import bytes_to_long,long_to_bytes
#设置模数
def GF(a):
    global p
    p = a
#乘法取模
def g(a,b):
    global p
    return pow(a,b,p)


def AMM(x,e,p):
    GF(p)
    y = random.randint(1, p-1)
    while g(y, (p-1)//e) == 1:
        y = random.randint(1, p-1)
        print(y)
    #p-1 = e^t*s
    t = 1
    s = 0
    while p % e == 0:
        t += 1
        print(t)
    s = p // (e**t)
    # s|ralpha-1
    k = 1
    while((s * k + 1) % e != 0):
        k += 1
    alpha = (s * k + 1) // e
    #计算a = y^s b = x^s h =1
    #h为e次非剩余部分的积
    a = g(y, (e ** (t - 1) ) * s)
    b = g(x, e * alpha - 1)
    c = g(y, s)
    h = 1
    #
    for i in range(1, t-1):
        d = g(b,e**(t-1-i))
        if d == 1:
            j = 0
        else:
            j = -math.log(d,a)
        b = b * (g(g(c, e), j))
        h = h * g(c, j)
        c = g(c, e)
    #return (g(x, alpha * h)) % p
    root = (g(x, alpha * h)) % p
    roots = set()
    for i in range(e):
        mp2 = root * g(a,i) %p
        #assert(g(mp2, e) == x)
        roots.add(mp2)
    return roots
def check(m):
    if 'flag' in m:
        return True
    else:
        return False
e = 997

p = 169192804045017094881483391290948160084538928031716323749363864952453968973507689162051165395748104110078160856791051809212190939432475142974911541618441458487669050818296365973889691415623806933502603345031427784795571665740530721508383685794846991682950112717404480456329219127191697671498037366841158723543
q = 107516396467746261711633898678341416690878446946218041251896502835689317784482747676107795221812916591321630759086326505565275611515776242892889358779953138176525964380991025435521861396436904104071935067377647496422254521013295763929078451759522826104921925202219553793049032407587608850233803508977340633609
c = 7296955328866123806615327249732627185102404227332181196296735121223965109231156544280256472492779759505533523060928048594910557437933201943976173955148680274140829916070075759044441331615135242760488256932238858269529909634447825461421412145996149026770528870738269768868586920051310346790350630656242240675615378779267818783700730455951708072880647986805110335263926177449091704517836266354071222826319675028232152825498040408774211261689801412297908590166114832939080331783731498956480608994534272354837899909567113733994622681549792329747132730450648055557829163328285671440063040320192447007187073122676185153708


mps = AMM(c,e,q)
# print(mps)
for mpp in mps:
        solution = str(long_to_bytes(mpp))
        if check(solution):
            print(solution)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
import random
from Crypto.Util.number import long_to_bytes

def AMM(o, r, q):

    g = GF(q)
    o = g(o)
    p = g(random.randint(1, q))
    while p ^ ((q-1) // r) == 1:
        p = g(random.randint(1, q))

    t = 0
    s = q - 1
    while s % r == 0:
        t += 1
        s = s // r
    k = 1
    while (k * s + 1) % r != 0:
        k += 1
    alp = (k * s + 1) // r
    a = p ^ (r**(t-1) * s)
    b = o ^ (r*alp - 1)
    c = p ^ s
    h = 1
    for i in range(1, t):
        d = b ^ (r^(t-1-i))
        if d == 1:
            j = 0
        else:
            j = - discrete_log(d, a)
        b = b * (c^r)^j
        h = h * c^j
        c = c^r
    result = o^alp * h
    return result

def findAllPRoot(p, e):
    proot = set()
    while len(proot) < e:
        proot.add(pow(random.randint(2, p-1), (p-1)//e, p))
    return proot

def findAllSolutions(mp, proot, cp, p):
    all_mp = set()
    for root in proot:
        mp2 = mp * root % p
        assert(pow(mp2, e, p) == cp)
        all_mp.add(mp2)
    return all_mp


c = 15433214846771804225704093824935372144929516863829752998270111032551363583267576397009018518790803908369965458162930857063271509296349091229352855725285388975497906340053281554202527432848881160125418406408621879995822551367228501163128699032015069585502994319524445505522625561831240862136447585120010288772692097621553249775117843166714346924868089146429002417223863834435968726551668931140147337199939823985783939085842479154589529244209712172799274024573845157268545992888944742377166586536479490962335287124809557709167220756920767331929168230518135523463578566851417486746667008938122693256033127001185017237773
p = 0xa892eb59b175bcf896be2176598f278437fe10ef032279f06e1092143acfb3c16b31811cca5286699595c2720c652ee64f8adc92c8b16a5601dd981d6f839ce9c0513db30de88c2ec6cae1a726acbd235ea946631bde633707d766287a2f075e9aace1606bd8b4f52d4f5b87dfb81f14fbc5338004575e9430257e180a169eff
q = 0xe3d47225b77e56129dc3fed716181845f89fa15b2eb35453ffdc0f05cdf57c0d90410911d209818e886b202bc4893ebe85a07ef670122f0e70092de1b7963c3b24a58c6a9ec9ed677db3473b1882d10d550e45c18fd57b85a70a5401a074d36760e85c7e6258f0ab08fa69cd433709910fad6e145f7b85f589e83d61d3baf6ad
n = p * q
e = 0x3
cp = c % p
cq = c % q
mp = AMM(cp, e, p)
mq = AMM(cq, e, q)
p_proot = findAllPRoot(p, e)
q_proot = findAllPRoot(q, e)
mps = findAllSolutions(mp, p_proot, cp, p)
mqs = findAllSolutions(mq, q_proot, cq, q)
print("mps=",mps)
print("mqs=",mqs)


flag = []
for mpp in mps:
    for mqq in mqs:
        solution = CRT_list([int(mpp), int(mqq)], [p, q])
        flag.append(solution)
print("flag=",flag)

for i in range(len(flag)):
    assert (((flag[i] ** 3)) % n == c)
    f = long_to_bytes(flag[i] - p)
    if b"DAS" in f:
        print(f)
        break
    else:
        print("No flag.")
        continue
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
from Crypto.Util.number import *
import gmpy2
import time
import random
from tqdm import tqdm
e = 65537
q= 12408795636519868275579286477747181009018504169827579387457997229774738126230652970860811085539129972962189443268046963335610845404214331426857155412988073
p= 12190036856294802286447270376342375357864587534233715766210874702670724440751066267168907565322961270655972226761426182258587581206888580394726683112820379
c= 68960610962019321576894097705679955071402844421318149418040507036722717269530195000135979777852568744281930839319120003106023209276898286482202725287026853925179071583797231099755287410760748104635674307266042492611618076506037004587354018148812584502385622631122387857218023049204722123597067641896169655595
n = p*q
def AMM(o, r, q):
    start = time.time()

    g = GF(q)
    o = g(o)
    p = g(random.randint(1, q))
    while p ^ ((q-1) // r) == 1:
        p = g(random.randint(1, q))
    print('[+] Find p:{}'.format(p))
    t = 0
    s = q - 1
    while s % r == 0:
        t += 1
        s = s // r
    print('[+] Find s:{}, t:{}'.format(s, t))
    k = 1
    while (k * s + 1) % r != 0:
        k += 1
    alp = (k * s + 1) // r
    print('[+] Find alp:{}'.format(alp))
    a = p ^ (r**(t-1) * s)
    b = o ^ (r*alp - 1)
    c = p ^ s
    h = 1
    for i in range(1, t):
        d = b ^ (r^(t-1-i))
        if d == 1:
            j = 0
        else:
            print('[+] Calculating DLP...')
            j = - discrete_log(d, a)
            print('[+] Finish DLP...')
        b = b * (c^r)^j
        h = h * c^j
        c = c^r
    result = o^alp * h
    end = time.time()
    print("Finished in {} seconds.".format(end - start))
    print('Find one solution: {}'.format(result))
    return result

def onemod(p,r):
    t=p-2
    while pow(t,(p-1) // r,p)==1:
        t -= 1
    return pow(t,(p-1) // r,p)

def solution(p,root,e):
    g = onemod(p,e)
    may = set()
    for i in range(e):
        may.add(root * pow(g,i,p)%p)
    return may

cp = c % p
mp = AMM(cp,e,p)
print('-------',mp)

mps = solution(p,mp,e)
for mpp in tqdm(mps):
     ai = [int(mpp)]
     mi = [p]
     m = CRT_list(ai,mi)
     flag = long_to_bytes(m)
     if b'moectf' in flag:
         print(flag)
         exit(0)

# moectf{Oh~Now_Y0u_Kn0W_HoW_RsA_W0rkS!}

高精度浮点实数域计算

题目是2022新生平台的LittleRSA。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
from Crypto.Util.number import *
from gmpy2 import *
from random import *
from flag import flag

m = bytes_to_long(flag)

while True:
    try:
        p = getPrime(512)
        q = next_prime(p+2**420)
        n = p*q
        phi = (p-1)*(q-1)
        d = randint(0,n**0.32)
        e = inverse(d,phi)
        c = pow(m,e,n)
        break
    except:
        continue

print("e = %d"%e)
print("n = %d"%n)
print("c = %d"%c)

'''
e = 101684733522589049376051051576215902510166244234370429058800153902445053536138419222096346715560283781778705047246555278271919928248836576236044123786248907522717751222608113597458768397652361813688176017155353220911686089871315647328303370846954697334521948003485878793121446614220897034652783771882675756065
n = 106490064297459077911162044548396107234298314288687868971249318200714506925762583340058042587392504450330878677254698499363515259785914237880057943786202091010532603853142050802310895234445611880617572636397946757345480447391544962796834842717321639098108976593541239044249391398321435940436125823407760564233
c = 92367575354201067679929326801477992215675304496512806779109227230237905402825022908214026985431756172011616861246881703226244396008088878308925377019775353026444957454196182919500667632574210469783704454438904889268692709062013797002819384105191802781841741128273810101308641357704215204494382259638905571144
'''

由题目可知，q = next_prime(p+2**420)，由此我们可以构造一个一元多项式方程
f = p*q-n = p*(p+2**420)-n，然后求解方程的根，此时p在根附近，于是在将根值往前推直至n%p=0为止，此时p为所求，接下来就是普通RSA解密：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from Crypto.Util.number import *
import gmpy2


e = 101684733522589049376051051576215902510166244234370429058800153902445053536138419222096346715560283781778705047246555278271919928248836576236044123786248907522717751222608113597458768397652361813688176017155353220911686089871315647328303370846954697334521948003485878793121446614220897034652783771882675756065
n = 106490064297459077911162044548396107234298314288687868971249318200714506925762583340058042587392504450330878677254698499363515259785914237880057943786202091010532603853142050802310895234445611880617572636397946757345480447391544962796834842717321639098108976593541239044249391398321435940436125823407760564233
c = 92367575354201067679929326801477992215675304496512806779109227230237905402825022908214026985431756172011616861246881703226244396008088878308925377019775353026444957454196182919500667632574210469783704454438904889268692709062013797002819384105191802781841741128273810101308641357704215204494382259638905571144

RF = RealField(600)
PR.<x> = PolynomialRing(RF)
f=x*(x+2**420)-n
p=int(f.roots()[1][0])

while n%p!=0:
    p-=1

p=10319402322686090423885467950006488404103970273239432095684700570087343967507257994593635913327166893587725950261323349433889479075061548042098460895952847
q=n//p
phi=(p-1)*(q-1)
d=gmpy2.invert(e,phi)
m=pow(c,d,n)
print(long_to_bytes(int(m)))

这里使用的定义域的方式：

1
2
RF = RealField(600)
PR.<x> = PolynomialRing(RF)

使用RealField(n)可以定义具有n位精度的实数域对象，根据题目可以知道f=x*(x+2**420)-n这个方程应该是求不出整数解的，可以求出一个大概的值，然后进行爆破就可以得到真正的p。

有限域开方

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import libnum
import gmpy2
import uuid

flag = "flag{" + str(uuid.uuid4()) + "}"
m = libnum.s2n(flag)
p = libnum.generate_prime(512)
q = libnum.generate_prime(512)
e = 7
while 1:
    p = libnum.generate_prime(512)
    q = libnum.generate_prime(512)
    if (p-1)%e==0 and (q-1)%e==0:
        break
n=p*q
c=pow(m,e,n)
print(flag)
print("n=", n)
print("p=", p)
print("q=", q)
print("c=", c)
print("e=", e)

通常这种情况下提供了素因子，但是$e$是$phi$的因子，发现$e$比较小，尝试在有限域里直接对密文进行开方：

R.<x> = Zmod(p)[]就是R.<x>=PolynomialRing(Zmod(p))构建多项式环的简便写法。

主要思想是从两个因子分别进行求解开方，至于最后的这个中国剩余定理，我的理解是如果对flag进行了pad处理使其二进制位长度大于因子的二进制位长度，进行中国剩余定理计算就可以求得在模$n$有限域下的根，否则求得的只是在缩小的域下的根，不一定会符合最终的要求：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import libnum

n= 111268707014190606153668884936880403903033189636811084190748342013582247181524984027156156751600082022630153178484737583696969737924460430658720077745043722334326797462286738281399709623421449519749033032872576930584696178280691506919932333135019009074258024671340323462243935601901807414139845477984309325673
p= 8829319990290192514343990745487198939057408737809583034593367387309247947172255218117653818474808992817642114513133973664718799951201403850947712246155847
q= 12602183082791809311087186602336708804360795363632019813154454176098834255153353255709663906079935612163450280376805433156799512016580171548991716849918159
c= 55885220678519981429180754853278867559851098342374272398008858602662839679132655399651996257139630339044262644051254178657597024166376700436329323381069534057238654849750900059528411876399512449220295427108651138742456755791576064932113111520299761812000274344470586300391658703812002906514018421358826196448
e= 7

R.<x> = Zmod(p)[]
f = x ^ e - c
f = f.monic()
res1 = f.roots()

R.<x> = Zmod(q)[]
f = x ^ e - c
f = f.monic()
res2 = f.roots()

for i in res1:
    for j in res2:
        m = crt(int(i[0]),int(j[0]),p,q)
        flag = libnum.n2s(int(m))
        if flag.startswith(b'flag'):
            print(flag)

位数变高之后，有限域开根就变得困难了，512位的小e还是容易开出来的。

Franklin-Reiter相关消息攻击

如果两个信息之间存在已知的固定差异（$f(x)=ax+b$），并且在相同的模数$n$下进行RSA加密，那么就有可能恢复出这两个消息。

有：
$$
C_1\equiv M_1^e(\bmod N)
$$

$$
C_2\equiv M_2^e(\bmod N)
$$

$$
M_1\equiv f(M_2)(\bmod N)
$$

可以知道$M_2$是$f(x)^e-C_1\equiv 0(\bmod N)$的一个解，也是$x^e-C_2\equiv 0(\bmod N)$的解，对两个多项式分别进行辗转相除法求得多项式的最大公因子，如果最大公因子是线性的（次数为1），即最大公因子可以表示为一个一次多项式（$x-M_2$），那么在模数下仍满足$x-M_2=0$，就求得了$M_2$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
from Crypto.Util.number import *
from tqdm import trange

n = 18993579800590288733556762316465854395650778003397512624355925069287661487515652428099677335464809283955351330659278915073219733930542167360381688856732762552737791137784222098296804826261681852699742456526979985201331982720936091963830799430264680941164508709453794113576607749669278887105809727027129736803614327631979056934906547015919204770702496676692691248702461766117271815398943842909579917102217310779431999448597899109808086655029624478062317317442297276087073653945439820988375066353157221370129064423613949039895822016206336117081475698987326594199181180346821431242733826487765566154350269651592993856883
c1 = 3089900890429368903963127778258893993015616003863275300568951378177309984878857933740319974151823410060583527905656182419531008417050246901514691111335764182779077027419410717272164998075313101695833565450587029584857433998627248705518025411896438130004108810308599666206694770859843696952378804678690327442746359836105117371144846629293505396610982407985241783168161504309420302314102538231774470927864959064261347913286659384383565379900391857812482728653358741387072374314243068833590379370244368317200796927931678203916569721211768082289529948017340699194622234734381555103898784827642197721866114583358940604520
c2 = 6062491672599671503583327431533992487890060173533816222838721749216161789662841049274959778509684968479022417053571624473283543736981267659104310293237792925201009775193492423025040929132360886500863823523629213703533794348606076463773478200331006341206053010168741302440409050344170767489936681627020501853981450212305108039373119567034948781143698613084550376070802084805644270376620484786155554275798939105737707005991882264123315436368611647275530607811665999620394422672764116158492214128572456571553281799359243174598812137554860109807481900330449364878168308833006964726761878461761560543284533578701661413931


def GCD(a,b):
    if b == 0:
        return a.monic()
    else:
        return GCD(b,a % b)
    
P.<x> = PolynomialRing(Zmod(n))

for e in trange(2**10):
    if isPrime(e):
        i=e
        f1 = (114*x+2333) ^ i - c1
        f2 = (514*x+4555) ^ i - c2
        if GCD(f1,f2)[0] != 1:
            print(long_to_bytes(int(n - GCD(f1,f2)[0])))

Boneh Durfee攻击

原理暂时先不了解了，直接记一下使用条件吧：

1
2
3
4
5
6
n=0xbadd260d14ea665b62e7d2e634f20a6382ac369cd44017305b69cf3a2694667ee651acded7085e0757d169b090f29f3f86fec255746674ffa8a6a3e1c9e1861003eb39f82cf74d84cc18e345f60865f998b33fc182a1a4ffa71f5ae48a1b5cb4c5f154b0997dc9b001e441815ce59c6c825f064fdca678858758dc2cebbc4d27L
d=random.getrandbits(1024*0.270)
e=invmod(d,phin)
hex(e)=0x11722b54dd6f3ad9ce81da6f6ecb0acaf2cbc3885841d08b32abc0672d1a7293f9856db8f9407dc05f6f373a2d9246752a7cc7b1b6923f1827adfaeefc811e6e5989cce9f00897cfc1fc57987cce4862b5343bc8e91ddf2bd9e23aea9316a69f28f407cfe324d546a7dde13eb0bd052f694aefe8ec0f5298800277dbab4a33bbL
m=random.getrandbits(512)
c=pow(m,e,n)=0xe3505f41ec936cf6bd8ae344bfec85746dc7d87a5943b3a7136482dd7b980f68f52c887585d1c7ca099310c4da2f70d4d5345d3641428797030177da6cc0d41e7b28d0abce694157c611697df8d0add3d900c00f778ac3428f341f47ecc4d868c6c5de0724b0c3403296d84f26736aa66f7905d498fa1862ca59e97f8f866cL

知道$e$,还知道$d≤n^{0.27}$，可以使用Boneh Durfee攻击：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
import time

"""
Setting debug to true will display more informations
about the lattice, the bounds, the vectors...
"""
debug = True

"""
Setting strict to true will stop the algorithm (and
return (-1, -1)) if we don't have a correct
upperbound on the determinant. Note that this
doesn't necesseraly mean that no solutions
will be found since the theoretical upperbound is
usualy far away from actual results. That is why
you should probably use `strict = False`
"""
strict = False

"""
This is experimental, but has provided remarkable results
so far. It tries to reduce the lattice as much as it can
while keeping its efficiency. I see no reason not to use
this option, but if things don't work, you should try
disabling it
"""
helpful_only = True
dimension_min = 7 # stop removing if lattice reaches that dimension

############################################
# Functions
##########################################

# display stats on helpful vectors
def helpful_vectors(BB, modulus):
    nothelpful = 0
    for ii in range(BB.dimensions()[0]):
        if BB[ii,ii] >= modulus:
            nothelpful += 1

    print (nothelpful, "/", BB.dimensions()[0], " vectors are not helpful")

# display matrix picture with 0 and X
def matrix_overview(BB, bound):
    for ii in range(BB.dimensions()[0]):
        a = ('%02d ' % ii)
        for jj in range(BB.dimensions()[1]):
            a += '0' if BB[ii,jj] == 0 else 'X'
            if BB.dimensions()[0] < 60:
                a += ' '
        if BB[ii, ii] >= bound:
            a += '~'
        print (a)

# tries to remove unhelpful vectors
# we start at current = n-1 (last vector)
def remove_unhelpful(BB, monomials, bound, current):
    # end of our recursive function
    if current == -1 or BB.dimensions()[0] <= dimension_min:
        return BB

    # we start by checking from the end
    for ii in range(current, -1, -1):
        # if it is unhelpful:
        if BB[ii, ii] >= bound:
            affected_vectors = 0
            affected_vector_index = 0
            # let's check if it affects other vectors
            for jj in range(ii + 1, BB.dimensions()[0]):
                # if another vector is affected:
                # we increase the count
                if BB[jj, ii] != 0:
                    affected_vectors += 1
                    affected_vector_index = jj

            # level:0
            # if no other vectors end up affected
            # we remove it
            if affected_vectors == 0:
                print ("* removing unhelpful vector", ii)
                BB = BB.delete_columns([ii])
                BB = BB.delete_rows([ii])
                monomials.pop(ii)
                BB = remove_unhelpful(BB, monomials, bound, ii-1)
                return BB

            # level:1
            # if just one was affected we check
            # if it is affecting someone else
            elif affected_vectors == 1:
                affected_deeper = True
                for kk in range(affected_vector_index + 1, BB.dimensions()[0]):
                    # if it is affecting even one vector
                    # we give up on this one
                    if BB[kk, affected_vector_index] != 0:
                        affected_deeper = False
                # remove both it if no other vector was affected and
                # this helpful vector is not helpful enough
                # compared to our unhelpful one
                if affected_deeper and abs(bound - BB[affected_vector_index, affected_vector_index]) < abs(bound - BB[ii, ii]):
                    print ("* removing unhelpful vectors", ii, "and", affected_vector_index)
                    BB = BB.delete_columns([affected_vector_index, ii])
                    BB = BB.delete_rows([affected_vector_index, ii])
                    monomials.pop(affected_vector_index)
                    monomials.pop(ii)
                    BB = remove_unhelpful(BB, monomials, bound, ii-1)
                    return BB
    # nothing happened
    return BB

"""
Returns:
* 0,0   if it fails
* -1,-1 if `strict=true`, and determinant doesn't bound
* x0,y0 the solutions of `pol`
"""
def boneh_durfee(pol, modulus, mm, tt, XX, YY):
    """
    Boneh and Durfee revisited by Herrmann and May

    finds a solution if:
    * d < N^delta
    * |x| < e^delta
    * |y| < e^0.5
    whenever delta < 1 - sqrt(2)/2 ~ 0.292
    """

    # substitution (Herrman and May)
    PR.<u, x, y> = PolynomialRing(ZZ)
    Q = PR.quotient(x*y + 1 - u) # u = xy + 1
    polZ = Q(pol).lift()

    UU = XX*YY + 1

    # x-shifts
    gg = []
    for kk in range(mm + 1):
        for ii in range(mm - kk + 1):
            xshift = x^ii * modulus^(mm - kk) * polZ(u, x, y)^kk
            gg.append(xshift)
    gg.sort()

    # x-shifts list of monomials
    monomials = []
    for polynomial in gg:
        for monomial in polynomial.monomials():
            if monomial not in monomials:
                monomials.append(monomial)
    monomials.sort()

    # y-shifts (selected by Herrman and May)
    for jj in range(1, tt + 1):
        for kk in range(floor(mm/tt) * jj, mm + 1):
            yshift = y^jj * polZ(u, x, y)^kk * modulus^(mm - kk)
            yshift = Q(yshift).lift()
            gg.append(yshift) # substitution

    # y-shifts list of monomials
    for jj in range(1, tt + 1):
        for kk in range(floor(mm/tt) * jj, mm + 1):
            monomials.append(u^kk * y^jj)

    # construct lattice B
    nn = len(monomials)
    BB = Matrix(ZZ, nn)
    for ii in range(nn):
        BB[ii, 0] = gg[ii](0, 0, 0)
        for jj in range(1, ii + 1):
            if monomials[jj] in gg[ii].monomials():
                BB[ii, jj] = gg[ii].monomial_coefficient(monomials[jj]) * monomials[jj](UU,XX,YY)

    # Prototype to reduce the lattice
    if helpful_only:
        # automatically remove
        BB = remove_unhelpful(BB, monomials, modulus^mm, nn-1)
        # reset dimension
        nn = BB.dimensions()[0]
        if nn == 0:
            print ("failure")
            return 0,0

    # check if vectors are helpful
    if debug:
        helpful_vectors(BB, modulus^mm)

    # check if determinant is correctly bounded
    det = BB.det()
    bound = modulus^(mm*nn)
    if det >= bound:
        print ("We do not have det < bound. Solutions might not be found.")
        print ("Try with highers m and t.")
        if debug:
            diff = (log(det) - log(bound)) / log(2)
            print ("size det(L) - size e^(m*n) = ", floor(diff))
        if strict:
            return -1, -1
    else:
        print ("det(L) < e^(m*n) (good! If a solution exists < N^delta, it will be found)")

    # display the lattice basis
    if debug:
        matrix_overview(BB, modulus^mm)

    # LLL
    if debug:
        print ("optimizing basis of the lattice via LLL, this can take a long time")

    BB = BB.LLL()

    if debug:
        print ("LLL is done!")

    # transform vector i & j -> polynomials 1 & 2
    if debug:
        print ("looking for independent vectors in the lattice")
    found_polynomials = False

    for pol1_idx in range(nn - 1):
        for pol2_idx in range(pol1_idx + 1, nn):
            # for i and j, create the two polynomials
            PR.<w,z> = PolynomialRing(ZZ)
            pol1 = pol2 = 0
            for jj in range(nn):
                pol1 += monomials[jj](w*z+1,w,z) * BB[pol1_idx, jj] / monomials[jj](UU,XX,YY)
                pol2 += monomials[jj](w*z+1,w,z) * BB[pol2_idx, jj] / monomials[jj](UU,XX,YY)

            # resultant
            PR.<q> = PolynomialRing(ZZ)
            rr = pol1.resultant(pol2)

            # are these good polynomials?
            if rr.is_zero() or rr.monomials() == [1]:
                continue
            else:
                print ("found them, using vectors", pol1_idx, "and", pol2_idx)
                found_polynomials = True
                break
        if found_polynomials:
            break

    if not found_polynomials:
        print ("no independant vectors could be found. This should very rarely happen...")
        return 0, 0

    rr = rr(q, q)

    # solutions
    soly = rr.roots()

    if len(soly) == 0:
        print ("Your prediction (delta) is too small")
        return 0, 0

    soly = soly[0][0]
    ss = pol1(q, soly)
    solx = ss.roots()[0][0]

    #
    return solx, soly

def example():
    ############################################
    # How To Use This Script
    ##########################################

    #
    # The problem to solve (edit the following values)
    #

    # the modulus
    N = 0xbadd260d14ea665b62e7d2e634f20a6382ac369cd44017305b69cf3a2694667ee651acded7085e0757d169b090f29f3f86fec255746674ffa8a6a3e1c9e1861003eb39f82cf74d84cc18e345f60865f998b33fc182a1a4ffa71f5ae48a1b5cb4c5f154b0997dc9b001e441815ce59c6c825f064fdca678858758dc2cebbc4d27
    # the public exponent
    e = 0x11722b54dd6f3ad9ce81da6f6ecb0acaf2cbc3885841d08b32abc0672d1a7293f9856db8f9407dc05f6f373a2d9246752a7cc7b1b6923f1827adfaeefc811e6e5989cce9f00897cfc1fc57987cce4862b5343bc8e91ddf2bd9e23aea9316a69f28f407cfe324d546a7dde13eb0bd052f694aefe8ec0f5298800277dbab4a33bb

    # the hypothesis on the private exponent (the theoretical maximum is 0.292)
    delta = 0.280 # this means that d < N^delta

    #
    # Lattice (tweak those values)
    #

    # you should tweak this (after a first run), (e.g. increment it until a solution is found)
    m = 4 # size of the lattice (bigger the better/slower)

    # you need to be a lattice master to tweak these
    t = int((1-2*delta) * m)  # optimization from Herrmann and May
    X = 2*floor(N^delta)  # this _might_ be too much
    Y = floor(N^(1/2))    # correct if p, q are ~ same size

    #
    # Don't touch anything below
    #

    # Problem put in equation
    P.<x,y> = PolynomialRing(ZZ)
    A = int((N+1)/2)
    pol = 1 + x * (A + y)

    #
    # Find the solutions!
    #

    # Checking bounds
    if debug:
        print ("=== checking values ===")
        print ("* delta:", delta)
        print ("* delta < 0.292", delta < 0.292)
        print ("* size of e:", int(log(e)/log(2)))
        print ("* size of N:", int(log(N)/log(2)))
        print ("* m:", m, ", t:", t)

    # boneh_durfee
    if debug:
        print ("=== running algorithm ===")
        start_time = time.time()

    solx, soly = boneh_durfee(pol, e, m, t, X, Y)

    # found a solution?
    if solx > 0:
        print ("=== solution found ===")
        if False:
            print ("x:", solx)
            print ("y:", soly)

        d = int(pol(solx, soly) / e)
        print ("private key found:", d)
    else:
        print ("=== no solution was found ===")

    if debug:
        print("=== %s seconds ===" % (time.time() - start_time))

if __name__ == "__main__":
    example()

Fault Injection攻击

参见近期比赛复现Vol.2题目Fault!Fault!

大步小步算法（中间相遇攻击）

大步小步算法是一种解决离散对数问题的方法，主要思路是中间相遇攻击。

对于求解离散对数：$y=g^x$，可以定义$x=im+j$，$m$的值为$[\sqrt n]$（取最接近的整数），$i,j$两个未知量的值在$m$的范围内，$m$的值决定算法计算出结果所需要的时间。

上式可变为：
$$
y=g^x=g^{im+j}\Rightarrow yg^{-mi}=g^j
$$
可以枚举$j$并计算$g^j$，之后枚举$i$计算$yg^{-mi}$的值，寻找$g^j$和$yg^{-mi}$相等的时候，得到我们想要的$i,j$值。代入$x=im+j$，求得$x$。

例如：
$$
g=17,p=509,y=g^x\bmod p=438
$$
求解$x$可以先枚举$j$：

$j$ $g^j\bmod p$
1 17
13 238
25 278
44 356
59 438

然后再对$i$枚举：

$i$ $yg^{-mi}\bmod p$
1 199
2 238
25 185
44 368
59 445

发现当$i=2,j=13$时求出的值相同，$x=2m+13$，$m=[\sqrt p]=23$，代入得到$x=59$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import gmpy2

g = 17
p = 509
y = 438
m = gmpy2.iroot(p, 2)[0] + 1

def BSGS(g, y, p):
    m = gmpy2.iroot(p, 2)[0] + 1
    s = {pow(g, j, p): j for j in range(m)}
    gs = pow(g, -m, p)
    for i in range(m):
        if y in s:
            return m*i + s[y]
        y = y * gs % p
    return None

print(BSGS(g, y, p))

一道例题：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from libnum import *
from gmpy2 import *

flag = "flag{big_and_small}"
a =  s2n(flag)
print(a)
m = a >> 16 << 16
p = generate_prime(1024)
g = 5
y = pow(g,a,p)
print(f"m = {m}")
print("g =", g)
print("p =", p)
print("y =", y)

#m = 2284117282071477104642039231536944446836572160
#g = 5
#p = 117061102100923822494374971271356428838245797897131713833497324490347029897509453234867144323074426358251468387009784964409490629474760505325012358657612358218630936761844853452752207513376553400570437173183245325363865380385444332069521612219917516076129636378470207520184345636826950827640031147884327602403
#y = 36330338308028857799449375552374448525921689801126735837219113818543840068942804412602004817464522788143367479566994608891201675854663902669408381421532026724950780811697860458785100785315287105361742435587380380417609918284799177467034247188834721615753928207192485082356827804549590792117891230279436213258

这题对$a$进行了移位，我们只求低位即可：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from libnum import *
from gmpy2 import *

m = 2284117282071477104642039231536944446836572160
g = 5
p = 104902226273988324576427043280139630983677434981789024088951356016939994526304119354435543071985576434554541580765980475339863817923069063092822047759376193281443519615446894334112711636144801809678018855363243259870518701678872861580764859373668820885488992214504069164524801155987240703746883303678741931657
y = 52132387535704137390458836331516728932227502073670482441552732370459557430160087735998314933737812855984248770511516860301143061293647971484306910330283327130899876785232553721318312581135200914802966676159312701878015611579070659605675344114334474279389271370807188711609485615280774712598619602033548583710

def BSGS(g, y, p, bound):
    x = iroot(bound, 2)[0] + 1
    s = {pow(g, j, p): j for j in range(x)}
    sg = pow(g, -m, p)
    for i in range(m):
        if y in s:
            return i*x + s[y]
        y = y * sg % p
    return None

a = BSGS(g, y*invert(pow(g,m,p), p)%p, p, 2**32) + m
print(n2s(int(a)))

#b'flag{big_and_small}'

这里顺带sage整理几个求解离散对数的函数：

这里求解的是以base为底，a的对数，ord为base的阶。

1
2
3
4
discrete_log(a,base,ord,operation)#通用求离散对数办法，纯暴力
discrete_log_rho(a,base,ord,operation)#Pollard-Rho算法
discrete_log_lambda(a,base,bounds,operation)#Pollard-kangaroo算法
bsgs(base,a,bounds,operation)#大步小步法