RSA学习笔记（六）

RSA学习是不能停止的，不是吗？

优秀资源：风二西/rsa_f2x - 码云 - 开源中国 (gitee.com)

共素数RSA

摘自文章：Common Prime RSA 笔记 | 独奏の小屋 (hasegawaazusa.github.io)

Wiener提出，在RSA中，当$p,q$为素数，当$p-1$和$q-1$有一个大数因子，可以反制维纳攻击，将这种$p-1$和$q-1$含有一个共同大素数因子的RSA称为共素数RSA。

对于公因子大素数$g$，让$p=2ga+1,q=2gb+1$，这里让$a,b$互素，且$h=2gab+a+b$也是素数，这样就可以确保$\gcd(p-1,q-1)=2g$，且$\frac {pq-1}2=gh$和$N=pq$大小接近。

发现：

$$\lambda(pq)=lcm(p-1,q-1)=lcm(2ga,2gb)=2gab$$ $$\phi(pq)=(p-1)(q-1)=2ga2gb=2g\lambda(pq)$$

所以还需要规定$e,d$还要和$\lambda(pq)$互素。

$$N=pq=(2ga+1)(2gb+1)=2g(2gab+a+b)+1=2gh+1$$ $$N-1=2gh$$

定义$\gamma$为公因子$g$相对于$N$的大小，$g=N^{\gamma}$，考虑$g≤N^{\frac 12}$，得到$0≤\gamma≤\frac 12$。

攻击方式

$\gamma$接近$\frac 12$

当$\gamma$接近$\frac 12$时，$g=N^\gamma$接近$N^{\frac 12}$，可以通过修改Pollard’s rho method的$x_i$函数进行攻击。

Pollard’s rho算法：能快速找到一个整数的非1非自身因子。

给定一个多项式$f(x)$和初始值$x_0$，计算数列$x_k,k≥1$，可以迭代计算$x_{k+1}\equiv f(x_k)\bmod n$，通常使用函数$f(x)=x^2+1$，寻找$k$满足$1＜(x_{2k}-x_k,n)＜n$，从而得到一个因数（不一定是质因数）。

这里使用的Pollard’s rho是修改迭代函数为$f(x)=x^{N-1}+3\pmod p$。

import random
from math import gcd

def rho(N):
    f = lambda x: (pow(x, N-1, N) + 3) % N
    while True:
        # 加快入环速度
        t = random.randint(2, N)
        h = f(t)
        step_times = 0
        step_limit = 2
        while True:
            if not step_times < step_limit:
                step_times = 0
                step_limit *= 2
                t = h
                h = f(h)
            p = gcd(abs(int(t) - int(h)), N)
            if p == N:
                break
            elif p > 1:
                return (p, N // p)
            else:
                h = f(h)
                step_times += 1

print(rho(84236796025318186855187782611491334781897277899439717384242559751095347166978304126358295609924321812851255222430530001043539925782811895605398187299748256080526691975084042025794113521587064616352833904856626744098904922117855866813505228134381046907659080078950018430266048447119221001098505107823645953039))

所以为了安全起见，$\gamma$的选取不能过于接近$\frac 12$。

已知$a,b$

已知$N=2g(2gab+a+b)+1$，构造方程：

$$4abg^2+2(a+b)g-N+1=0$$

可以在$log(N)$时间内分解$N$：

a = 1185431345934512
b = 1989628969125971
N = 54692260436051338814890781701826055707958209029414126894070449935683071253184867947357262267840171428710181955973010913204514025135188192484651672240708141692701667242130748316666406528479191422804307020656050201187035928715833163999813216597718706449260040885862566373392398826670863398295350419792842640631
P.<g> = ZZ[]
f = 4 * a * b * g ^ 2 + 2 * (a + b) * g - N + 1
g = f.roots()
if g:
    g = g[0][0]
    p = 2 * g * a + 1
    q = 2 * g * b + 1
    assert p * q == N

上面是sage的代码，在python中也可以使用sympy库进行求解。

已知$g$

当$g≥a+b$时，也可以在多项式$log(N)$时间内分解$N$。

由于$a,b$相对差距较小，所以这个实现条件约等于$g＞N^{\frac 14}$。

简单证明：

①如果$g>a+b$，那么给$N,g$，令$M=\frac {N-1}{2g}$，$c=a+b$，可变换方程：

$$N=2g(2gab+a+b)+1\Rightarrow M=2gab+c$$

因为$c=a+b<g$，所以可以规约为：

$$c=M\pmod g$$

$M$是可以通过已知求得的，所以可以求得$c=a+b$。

将$b=c-a$代入$N=2g(2gab+a+b)+1$，整理得到二次方程：

$$2ga^2-2gac+\frac{N-1}{2g}-c=0$$

可以求解$a,b$。

g = 2056971706333850947354991471886113601423457483931388832864204860321308350537317091564919029078296379733989138742162694786565228112885684303
N = 67324909911911622626246005558967775211455024820506932698435813321567574468019013664789401988015894964099052816176029553245881317276340043887466584645914352982274378611180595397686920214079479901514703963131435008906250160656759300390805929849374653321934393399433471228218819498373221757779799476717494079667
M = (N - 1) // (2 * g)
c = M % g
P.<a> = ZZ[]
f = 2 * g * a ^ 2 - 2 * g * c * a + M - c
a = f.roots()
if a:
    a, b = a[0][0], a[1][0]
    p = 2 * g * a + 1
    q = 2 * g * b + 1
    assert p * q == N

②如果$g=a+b$，将方程$N=2g(2gab+a+b)+1$变换，得到：

$$\frac {N-1}{4g^2}=ab+\frac 12$$

替换$g=a+b$，同样得到二次方程：

$$a^2-ga+\frac{N-1}{4g^2}-\frac 12=0$$

可以求解$a,b$。

g = 2855372645569408464444580237486670388029956719716115953907612135874419892154982850222965560661211729647325085879529571229774148545656169021
N = 159549169988238873893531105042878385551537587717347282632324748268846735710748763722602882823022008548774298858161130258369850715542192739582830583643642436399008902770027668038725347353393047833875066622910131525247842517372845617227325882916166114361718015983671803859502931814932543107911548450229250776542101141849788751722460468073974316977656001286989710480324512919121409123619799426221232443036698458643438020098037548757403
M = (N - 1) // (2 * g)
P.<a> = ZZ[]
f = 2 * a ^ 2 - 2 * g * a + (N - 1) // (2 * g ^ 2) - 1
a = f.roots()
if a:
    a, b = a[0][0], a[1][0]
    p = 2 * g * a + 1
    q = 2 * g * b + 1
    assert p * q == N

因式分解$N-1$

当$\gamma$过小时，也就是$g=N^\gamma$过小，已知$N-1=2gh$，所以容易分解出$g$。

可以使用yafu对$\frac {N-1}2$进行分解，当$\gamma$约为0.1分解迅速。

dp泄露优化运算

在之前对dp泄露的学习中，主要思路是寻找一个$(0,e)$范围内的数求出的$p$能被$n$整除，得到$p,q$，但是实际解题中可能会出现$e$较大无法实现遍历的情况。这里就需要换个思路，从公式开始：

$$dp=d\bmod (p-1)$$ $$dpe=de\bmod (p-1)$$

我们还知道在RSA中满足$ed\equiv 1\pmod {phi}$，所以就有：

$$dpe=k_1(p-1)(q-1)+1+k_2(p-1)$$ $$dpe\bmod (p-1)=1$$ $$dpe=1+k(p-1)$$

在这里想要得到$p$，可以引入一个小值，根据欧拉定理：

$$a^{dpe}\bmod p=a$$ $$a^{dpe}-a=kp$$

可以观察到这里我们得到了$kp$，猜测可以求最大公因数来分解$p$，但是等式左边数太大没法求，需要化简一下，对等式同时模$n$：

$$a^{dpe}\bmod n-a\bmod n=kp\bmod n$$

由于$a$为小值，所以模$n$可以忽略，式子可以继续约减：

$$a^{dpe}\bmod n-a=kp+tpq=p(k+tq)$$

这下就得到$p$了：

import gmpy2
from Crypto.Util.number import *

n=21082670655765549884141819030007477108817658491958447251093801606827378015193183556560209676301855782713945178801713329661806980350484336319501413676462721753617996650487260006277044569746515551628506171091864342450615204710909617213978934393796221059684745808609899726269440298467864210863676979548190338924347634582692877394725374050517129243153142378346805606271337107330079283903711934907556258371771296950387976466583088754428126811641919095093323861779757831735110475049041450467360956711753885749069213977988577474980908072886657312709638534439681970984892821295832244323484076015299706780730054259776289020073
e=227377575879236824204281410135954861057
c=6368844360770044065150863396025662365155502227885525414015899617873465570290962316793808652847966917288200533610888150296569398786221968313804433853973895050500047290594629145605691322774722288754210780839109143486593482056923144392302466689379808824428405899211768990461310994387371799453723737500551143011177693577553528635057224811215176751994816694977332094928055171013645814003852621025492016382232115360052709760826719719603090013890469847895044614844574827180553216574311861332786344893596925513203381508229215196083858585314897856706740524495343023412105682917436282434071209745628995859549842004616614138608
dp=100150669785686503003332919730638996199363050791478547879134303971915339826903078879298353164463026902022293972967277473554856528177131280448828956070411149953672091355923938438702058340437083863887955193747852625107455812377390167673142537349913995448677427628216181698994450948775040561560290230800139124049
p=gmpy2.gcd(n,pow(6,dp*e,n)-6)
q=n//p
phi=(p-1)*(q-1)
d=gmpy2.invert(e,phi)
m=pow(c,d,n)
print(long_to_bytes(m))

这种做法也可以应用于普通的dp泄露攻击中，我觉得这种做法比目前流行的做法方便很多。

一点想法

当然，如果我们能确定明文$m<p$，那么也不需要再进行后续的模板计算$m$了，因为：

$$dp\equiv d\pmod {p-1}$$

如果我们把计算定义到模$p$域上，那么给了$dp$就相当于已知了私钥$d$，所以可以直接拿来解$m$：

import gmpy2
from Crypto.Util.number import *

n=21082670655765549884141819030007477108817658491958447251093801606827378015193183556560209676301855782713945178801713329661806980350484336319501413676462721753617996650487260006277044569746515551628506171091864342450615204710909617213978934393796221059684745808609899726269440298467864210863676979548190338924347634582692877394725374050517129243153142378346805606271337107330079283903711934907556258371771296950387976466583088754428126811641919095093323861779757831735110475049041450467360956711753885749069213977988577474980908072886657312709638534439681970984892821295832244323484076015299706780730054259776289020073
e=227377575879236824204281410135954861057
c=6368844360770044065150863396025662365155502227885525414015899617873465570290962316793808652847966917288200533610888150296569398786221968313804433853973895050500047290594629145605691322774722288754210780839109143486593482056923144392302466689379808824428405899211768990461310994387371799453723737500551143011177693577553528635057224811215176751994816694977332094928055171013645814003852621025492016382232115360052709760826719719603090013890469847895044614844574827180553216574311861332786344893596925513203381508229215196083858585314897856706740524495343023412105682917436282434071209745628995859549842004616614138608
dp=100150669785686503003332919730638996199363050791478547879134303971915339826903078879298353164463026902022293972967277473554856528177131280448828956070411149953672091355923938438702058340437083863887955193747852625107455812377390167673142537349913995448677427628216181698994450948775040561560290230800139124049
p=gmpy2.gcd(n,pow(6,dp*e,n)-6)
m=pow(c,dp,p)
print(long_to_bytes(m))

$eq\equiv 1\pmod{p}$攻击

实现是这样的：

while 1:
    p = libnum.generate_prime(512)
    q = libnum.invmod(e, p)
    if gmpy2.is_prime(q):
        break

简单推理一下：

$$eq= 1+k p$$ $$epq=p+kp^2$$ $$4ken+1=4k^2p^2+4kp+1=(2kp+1)^2$$

$k$不大的时候，可以爆破，求得$p$：

import gmpy2
from Crypto.Util.number import *

n = 151913753330829779363367789673597253978048533620182497747132604072824540323836504426239804081241567027736270622449855295828160124814438423488754623643933533204599629980943193643551097211543531050709295786991507076525522113596966181359611779370941657496310723772521162794885147507996772774824707716432570452403
c = 81337035006499494768796081417947319039576994033747417935302379695580101083380316241731843797604372408639889020709527188720865744673275465620672796534287478878369491626837235522013863942055110385349674739290041899890532139474221925394370465243005480529009377659505930236007519685951328899152586367092905739547
e = 65537

for k in range(e,1,-1):
    tmp=1+4*e*n*k
    r,s=gmpy2.iroot(tmp,2)
    if s:
        print(k)
        p=(r-1)//(2*k)
        break
    
q=n//p
phi=(p-1)*(q-1)
d=gmpy2.invert(e,phi)
m=pow(c,d,n)
print(long_to_bytes(int(m)))

$invert(p,q),invert(q,p)$泄露

pinv = gmpy2.invert(p, q)
qinv = gmpy2.invert(q, p)

已知$phi$

$$\begin{cases}p^{-1}p\equiv 1\pmod q\\q^{-1}q\equiv 1\pmod p\end{cases}$$

可以构建方程组求解：

import gmpy2
import libnum
import sympy

e = 65537
phi = 98229890689284912089157469320766533791995454910192245906925225670195818406213091522014617153508961926275340781138620279250871290490873820930414041953789279215193143006258942396087516751221787778485548659324564021606857623517821094138942269722883193348321303203233297445863729579811837657596605227772947768080
c = 20962649632855841088231373617057756055795073764826640468523117901947795570191538819434503909156672490273164571512993513456958532656058586024929659925126633909116178539998046336126716552474329852733455142305774551889295732464058882055828240152313078876041502758035842268176085433166898402827107253041748761839
p1 = 9278107621912111502897373552245845971739215588186465664496052995069919521787577753719495078780428864673954897367324663969013147023292682927178810784000630
q1 = 1888254482465757731585051308443185269531054824465971018646930721602120686960471254468405491101133069657301849498853283496490197246867160626826664577479894

p = sympy.symbols('p')
q = sympy.symbols('q')
f1 = p1 * p + q1 * q - 1 - p * q
f2 = (p - 1) * (q - 1) - phi
pq = sympy.solve([f1, f2], [p, q])
p = (pq[1][0])
q = (pq[1][1])
n = p * q
d = gmpy2.invert(e, phi)
m = pow(c, int(d), int(n))
print(libnum.n2s(int(m)))

已知$e,d$

遍历求解：

p1= 97656077071270197914785321263555295398683131376741279762023269752748032446917564784973530125497364995807144029368514414686812526434051342203148587548869714362392149049199571504107521002327725946079985444111124295505554135534791115828253787424842366876671199369039838469527700802048214832651973914615014793434
q1= 39346433326508892813150536257919182675225716975877498778468626666265319210911093864684217868389854790505495572529483241732413112619372606685748980784688553398278315703981192828280813023954306657912647634729098972136705540024619937335024692325575092594446158807029214788285050359047781519982931005193436553164
d= 13982206046522221689006948762375515892286671914597370002120896052266280914021885507889789598565474633660266121772385535949954654478591294089438758529565129571215458972099907852294228150256147359389834287744392536374255187239416806275744570662248893522074706262730636918053134499811014564470545819778015584753659205965073426501902382941707278051440033966771995946005135307665514929445057737703258338176975384665564222178081103000139437728386604903444661523794053627011660691133434801451813006697587015285256584456452850523012639216643290392821724776006584995172392016576699118057963479401425915020921720042294733838337
e= 65537
c= 13235028416286929520160810669079534965765212020283157183506000327481359053022089166886460096491962822206497841520395685863417454564782563002357235184688540460362871929729166716636143577909568477095287025119038866855396071566223268818065261766734924398179202129826294778702597618099032226876112654250403885134458713788086798004345808257538019652522209856828455756694071035571007463724635674562294436700546760823969026149755043139120494280665479551005965472492371422616161244117935709924215027494208339801474837499105640735072011664365891939864038690374913413240430803230704628673411107455299647523018041707726365601893
phi=16033591784555690839033601641933863819862675274128076669740291933394723758788009352788583792713912217703542497595829207543956085361202366334328595985304274508481724780489075813896396106581346750635700757907059309496772854950127016253035273087412614339939285140337656629662098889174735118835844089284571097739371559723771996652001268037000802800553357819154820449211550823391567258031927998247803081416301018071582489307196620368843405837158330864310134068534625079672858460155583646814590373389203531298114865193213718937686858488611970228064763729093357377320357215677266414126622909951204685623007887150263682313184

import sympy
import gmpy2
import libnum

p = sympy.symbols('p')
q = sympy.symbols('q')
for k in range(e,1,-1):
    if (e*d-1)%k==0:
        phi=(e*d-1)//k
        print(phi)
        print(k)
        try:
            f1 = (p - 1) * (q - 1) - phi
            f2 = q1 * q + p1 * p - p * q - 1
            c1 = sympy.solve([f1, f2], [p, q])
            p = int(c1[1][0])
            print(p)
            q = int(c1[1][1])
            print(q)
            break
        except:
            continue


n=p*q
m=pow(c,d,n)
print(libnum.n2s(m))

变形

p1 = gmpy2.invert(p, q)-q
q1 = gmpy2.invert(q, p)-p

知道$phi$：

import gmpy2
import libnum
import sympy


e= 65537
phi= 153526208973286457920479376447622546334064100044507994947932972947056601063592943200823498761273906685365632931851181655405163001559274037884326450823919640514634182433839200397012200852498394907718538509040702321554785017853186412670980153149366201235845635320919980848682144430326738538775142630420821822960
c= 45404866893480020010799585996788902747901025575540655893265307771493992595444010006643443639614244962689890716458286889196613281948316354607540970243071548624687668875908716690311985915197944902857988253464493234021796216923039652003206533181461815990153623141243415592031190212825941414063970906141991979873
p1= -11065332660497699825690181811194036599102527446002445107664152730915570503984416572510322521934240155183790860963748546243230674204488210677533190963726687
q1= -1516147835082043403956384639695921948575165356837676821387726774997348251688472575043896909177257102383857065378906848506987546488480618072779943395893217

p = sympy.symbols('p')
q = sympy.symbols('q')
f1 = p1 * p + q1 * q - 1 +p * q
f2 = (p - 1) * (q - 1) - phi
pq = sympy.solve([f1, f2], [p, q])
p = (pq[1][0])
q = (pq[1][1])
n = p * q
d = gmpy2.invert(e, phi)
m = pow(c, int(d), int(n))
print(libnum.n2s(int(m)))

$c\equiv m\pmod p$类型

从小鸡块大佬那里发现的几个题，简单复现一下。

这里同余的前提是$m>>p$，所以无法爆破求解。

easy_mod

from Crypto.Util.number import *
from random import *

table = "01234567"
p = getPrime(328)
flag = b"NSSCTF{" + "".join([choice(table) for i in range(70)]).encode() + b"}"
c = bytes_to_long(flag) % p

print("p =",p)
print("c =",c)

'''
p = 501785758961383005891491265699612686883993041794260611346802080899615437298977076093878384543577171
c = 327005346153237517234971706274055111857447948791422192829214537757745905845319188257204611848165263
'''

flag的构成为0-7的随机组成，ASCII的范围为48-53，可以将flag的头尾去掉，只保留中间的ASCII较小部分，然后去造格求解：

我们根据位数位置去掉前缀premix和后缀suffix：

$$m=256^{71}pre+256m_0+suf$$

根据$c\equiv m\pmod p$，联立化简一下得到：

$$m_0\equiv 256^{-1}(c-256^{71}pre-suf)\pmod p$$

这里我们求得的$m_0$是flag内部的0-7随机组合得到的，但是我们未知某个位置上为哪个数，可以写为多项式形式：

$$m_0 = \sum_{i=0}^{69}256^is_i =s_0 + 256s_1 + 256^2s^2 + ... + 256^{69}s^{69} \pmod p$$

造格：

$$(s_0,s_1,...,s_{69},1,k) \left(\begin{matrix} 1 &&& &&1\\ &1&& &&256\\ &&\ddots& &&\vdots\\ &&&1 &&256^{69}\\ &&& &1&-m_0\\ &&& &&p\\ \end{matrix}\right) = (s_0,s_1,...,s_{69},1,0)$$

这样事实上得不到最终由48-53数值构成的向量，我自己试了一次，确实不行。

可以尝试继续减小向量的大小，设：

$$t_i=s_i-48$$

同时处理一下刚刚的同余式：

$$m_0\equiv \sum^{69}_{i=0}256^i(t_i+48)\pmod p$$ $$m_1=m_0-48\sum^{69}_{i=0}256^i=\sum^{69}_{i=0}256^it_i+kp$$

用这个新的式子造格，还需要进行配平操作：

from Crypto.Util.number import *

p = 501785758961383005891491265699612686883993041794260611346802080899615437298977076093878384543577171
c = 327005346153237517234971706274055111857447948791422192829214537757745905845319188257204611848165263
prefix = b"NSSCTF{"
suffix = b"}"
length = 78 - len(prefix) - len(suffix)

c -= 256^(len(suffix) + length) * bytes_to_long(prefix)
c -= bytes_to_long(suffix)
c = c * inverse(256,p) % p

L = Matrix(ZZ,length+2,length+2)
for i in range(length):
    L[i,i] = 1
    L[i,-1] = 256^i
    c -= 256^i*48

L[-2,-2] = 4
L[-2,-1] = -c
L[-1,-1] = p
L[:,-1:] *= p
res = L.BKZ()
flag = ""
for i in res[:-1]:
    if(abs(i[-2]) == 4 and all(abs(j) < 8 for j in i[:-2])):
        for j in i[:-2][::-1]:
            flag += chr(48 + abs(j))
print(flag)