multi-power-RSA 笔记

\(N = p^rq\) Attack (\(p\) 多次幂攻击)

加解密代码

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

def hensel_lifting(m, c, e, p, k):
    """
    Require pow(m, e, p) == c % p

    Get m % (p ** k)
    """
    assert pow(m, e, p) == int(c % p)
    mk = int(m % p)
    t = int(m * pow(c, -1, p) % p)
    pk = int(p)
    for _ in range(1, k):
        pk *= p
        x = int((c - pow(mk, e, pk)) % pk)
        y = int((x * t * int(pow(e, -1, p))) % pk)
        mk = mk + y
    return mk

def genkey(pbits: int, ebits: int, r: int, coprime: bool = True):
    r"""
    Generate a pair public key (n, e) and private key (p, q, d)

    where $n=p^rq, ed=1\pmod{\varphi}$

    Param:
        pbits - p bits and q bits, where p and q is balanced
        ebits - e bits
        r - power of p
    
    Return:
        (n, e)
    """
    p = getPrime(pbits)
    q = getPrime(pbits)
    phi = p ** r * (p - 1) * q
    n = p ** r * q
    e = random.getrandbits(ebits)
    while coprime and gcd(e, phi) != 1:
        e = random.getrandbits(ebits)
    return (n, e), (p, q)

def encrypt(m, e, n):
    return pow(m, e, n)

def decrypt(c, e, p, q, r):
    dp = pow(e, -1, p - 1)
    dq = pow(e, -1, q - 1)
    mp = pow(c, dp, p)
    mq = pow(c, dq, q)
    mpr = hensel_lifting(mp, c, e, p, r)
    m = crt([mpr, mq], [p ** r, q])
    return int(m)

pbits = 128
ebits = 64
r = 5
(n, e), (p, q) = genkey(128, 64, r)
m = 1234567890123456789012345678901234567890
c = encrypt(m, e, n)
_m = decrypt(c, e, p, q, r)
print(f'{_m = }')

The First Attack

Lu, Zhang and Lin theorem (线性模多项式方程定理)

设 \(N\) 整除 \(p^u\)，使得 \(p\geq N^\beta\ (0<\beta\leq 1)\)，

设 \(f(x,y)\in\Z[x,y]\) 为齐次线性多项式，

如果满足条件： \[ \gamma1+\gamma2<uv\beta^2 \] 我们就可以在多项式时间内找到方程 \(f(x,y)=0\pmod{p^v}\) 所有的解 \((x,y)\)，其中 \(\gcd(x,y)=1,|x|<N^{\gamma1},|y|<N^{\gamma2}\)。

Legendre's theorem (连分数定理)

如果存在 \(\xi \in \Q^*\)，\(a,b\in \Z\) 且 \(\gcd(a,b)=1\) 满足 \[ \left|\xi-\frac{a}{b}\right|<\frac{1}{2b^2} \] 那么 \(\displaystyle\frac{a}{b}\) 就是 \(\xi\) 的一个有理近似。

原理

设 \(N=p^rq\) 是 RSA 中的模数，且公钥 \(e\) 满足方程 \(ex-\varphi(N)y=z\)，如果 \[ |xz|<N^{\frac{r(r-1)}{(r+1)^2}} \] 那么我们可以在多项式时间内分解 \(N\)。

证明忽略，最后得到结论 \[ \gcd(ex-z,N)=\gcd(p^{r-1}(p-1)(q-1)y,p^rq)=g \] 可以确定 \(p\)，因此将 \(N\) 按如下情况分解：

• 当 \(g=p^{r-1}\) 时，\(p=g^{\frac{1}{r-1}}\)
• 当 \(g=p^r\) 时，\(p=g^{\frac{1}{r}}\)
• 当 \(g=p^{r-1}q\) 时，\(p=\frac{N}{g}\)

因此我们分解了 \(N\)。

在算法中我们转换方程为 \(ex-z\equiv0\pmod{p^{r-1}}\)，使用 Coppersmith's technique，选择参数 \(m=7,t=6\) 构造维度为 36 的格，并且选择 \[ X=\left[N^{\frac{r(r-1)}{(r+1)^2}}\right] \] 我们使用多项式 \(f(x1,x2)=x1+ex2\) 构造格，应用 LLL 算法寻找到小根 \((x_1,x_2)\) 满足 \(f(x_1,x_2)\equiv0\pmod{p^{r-1}}\)。

这样我们就可以计算出 \[ g=\gcd(x_1+ex_2,N) \]

脚本

当方程 \(ex-z\equiv0\) 的解已知一个时，我们仅需一元 Copper 即可求解：

N=539779851369541956878655738599584730199799866957191805784596190682932284216781781433367450841202917758999300635019369629627621029957135109806205877317954671312041249493462048283611940752235036153024920172209763260723728345918562258401803973624430150143563078517485996070862532682695228590709019451174548520135142052216785774589096706631010293690859363524584240662502290912412366366114571976050857239915691266377257797199583543940504695517331512813468837128344612227973709974625418257243011036826241599265375741977853552204640800449679679351666009764297016524814036295707311913711955324055690490892097177271718850857268982130811714517356073266905474635370690445031512184247179039751734276906533177939993769044135143389748416635981226449566039039202521305851567296884751935162651063209779647359922622084851547605090230221057349511482738300221222563908357379545905837110168948295030747460300104202323692732549831403834387939156877086852393515817984772384147449841124275061609701453997579569931391166586163299940486204581696722731952467570857217406030804590055255431828403195798003509083922294733709507134156466158642941338493323430671502043066148246348074878064089651235355282144209668143249348243220714471988019011613749340243917652821
e=8179300978753084587812861894047395225516049110376948812109811319430275614612773726672345893359691900281432484382670047044697374818043512731533402576374645405477207239801498428774783768163880078495448747421425078521981578408638790336528372019271073712013371141939808017049399434858687299480461753638164719404612128939787055797762174745092074547412183349192156638711750872083313795551439465507724807626674514935170104573715458782366469587138508845980490673890245713729782917089910271980557159592807350504157192913530007199510144004848020221181558472160543018733124225266127379373751910439604459368078652499029070936707349862139053913745186413782066470461478961703013591655136140060879250067379283913798867648758171004535775565306842444545755351202796833177560656564652632975685912935281581268141803696686952259539945588609591385807620108279333498170028167338690235117003515264281843953984997958878272347778561933726792473981855755454522886321669676790813189668084373153897754540290867346751033567500922477317530445967753955221454744946208555394588111484610700789566547507402309549957740815535069057837915204852490930168843605732632328017129154852857227895362549146737618906180651623216848500491438142456250653458053922622240299736136335179639180898730269690699965799644757774472147210271111150769048976871249731156387939260749192370361488285775377622944817570292095201906142567403539151179209316853493906909989301225903409448461436855145
b=17623328397444755087284107444487160871617682792372566887446834913712379373851213638071138745775127796589871734472781755930251379295485892067473329763997583502625804363418069062645997342172778252731889437
r = 7
idx = (r*(r-1)) / ((r+1)*(r+1))
P.<x> = PolynomialRing(Zmod(N))
f = e * x - b
roots = f.monic().small_roots(beta=idx)
print(roots)

当方程 \(ex-z\equiv0\) 的两个解都未知时，我们需要使用二元 Copper 求解：

import itertools

def small_roots(f, bounds, m=1, d=None):
    if not d:
        d = f.degree()

    R = f.base_ring()
    N = R.cardinality()

    f /= f.coefficients().pop(0)
    f = f.change_ring(ZZ)

    G = Sequence([], f.parent())
    for i in range(m + 1):
        base = N ^ (m - i) * f ^ i
        for shifts in itertools.product(range(d), repeat=f.nvariables()):
            g = base * prod(map(power, f.variables(), shifts))
            G.append(g)

    B, monomials = G.coefficient_matrix()
    monomials = vector(monomials)

    factors = [monomial(*bounds) for monomial in monomials]
    for i, factor in enumerate(factors):
        B.rescale_col(i, factor)
    B = B.dense_matrix().LLL()
    
    print("LLL done")

    B = B.change_ring(QQ)
    for i, factor in enumerate(factors):
        B.rescale_col(i, 1 / factor)

    H = Sequence([], f.parent().change_ring(QQ))
    for h in filter(None, B * monomials):
        H.append(h)
        I = H.ideal()
        if I.dimension() == -1:
            H.pop()
        elif I.dimension() == 0:
            roots = []
            for root in I.variety(ring=ZZ):
                root = tuple(R(root[var]) for var in f.variables())
                roots.append(root)
            return roots

    return []


a = 68823589009570846705623113091430161477799561031575612135516351257937127579444
b = 76549169049080319489163980188997771079750670038002598167961495813084486794567
c = 99215492498952976642031024510129092308042391285395704888838178561670205468882
p = 12775129699668988473759438271274836254349225413222075887429387682336494103348583050672280757042383792640084197832605411237644937815012509935794275313643031


x1 = 664406108637237317109372377546194289077117372665932150107678757303243619963182733408311739663233548725195723884930183442927293177078507
x2 = 269639459994260392015105629865659210391196381140872041084737976688017858324308345528702533444116282512066085580909718347778743354754352


HiddenBits = 64
bound = 2**HiddenBits - 1

P = PolynomialRing(Zmod(p), names=('x', 'y'))
(x, y) = P.gens()

f = y + (x2 << HiddenBits) - a * (x + (x1 << HiddenBits))**2 - b * (x+(x1<<HiddenBits)) - c
bounds = (bound, bound)
roots = small_roots(f, bounds)
print(roots)

The Second Attack

原理

设模数 \(N=p^rq\)，私钥 \((d_1,d_2)\)，如果 \[ |d_1-d_2|<N^{\frac{r(r-1)}{(r+1)^2}} \] 那么我们可以在多项式时间内对 \(N\) 进行分解。

证明如下：

假设 \(e_1d_1-k_1\phi(N)=1\) 、\(e_2d_2-k_2\phi(N)=1\) 且 \(e_1>e_2\)

因此 \(e_1d_1\equiv 1\pmod{\phi(N)}\)，\(e_2d_2\equiv 1\pmod{\phi(N)}\)

前一个式子乘以 \(e_2\)，后一个式子乘以 \(e_1\)，相减得到 \[ e_1e_2(d_1-d_2)\equiv e_2-e_1\pmod{\phi(N)} \] 因为 \(\phi(N)=p^{r-1}(p-1)(q-1)\)，故 \(e_1e_2(d_1-d_2)\equiv e_2-e_1\pmod{p^{r-1}}\)

现在，只需考虑以下这个线性方程 \[ e_1e_2x-(e_2-e_1)\equiv 0\pmod{p^{r-1}} \] 其中 \(d_1-d_2\) 是该方程的根。

进一步假设 \(|d_1-d_2|<N^\delta\)，运用定理一，其中 \(u=r,v=r-1\ 且\ \beta=\frac{1}{r+1}\)，可以推导出如果 \(\delta<uv\beta^2=\frac{r(r-1)}{(r+1)^2}\)，那么 \(x=d_1-d_2\) 可在多项式时间内求解。

也就是说如果 \(|d_1-d_2|<N^{\frac{r(r-1)}{(r+1)^2}}\)，可以在多项式时间内求解方程。

通过计算 \[ \gcd(e_1e_2x-(e_2-e_1),N)=\gcd(p^{r-1}(p-1)(q-1)y,p^rq)=g \]

证明忽略，最后得到结论 \[ \gcd(e_1e_2x-(e_2-e_1),N)=\gcd(p^{r-1}(p-1)(q-1)y,p^rq)=g \] 可以确定 \(p\)，因此将 \(N\) 按如下情况分解：

• 当 \(g=p^{r-1}\) 时，\(p=g^{\frac{1}{r-1}}\)
• 当 \(g=p^r\) 时，\(p=g^{\frac{1}{r}}\)
• 当 \(g=p^{r-1}q\) 时，\(p=\frac{N}{g}\)

因此我们分解了 \(N\)。

实例

尝试求解 \(N=p^2q\)，其中 \[ N = 6093253851486120878859471958399737725885946526553626219\\ e1 = 2749600381847487389715964767235618802529675855606377411\\ e2 = 3575081244952414009316396501512372226545892558898276551 \] 构建方程 \(f(x)=e_1e_2x-(e_2-e_1)\equiv 0\pmod{p^{r-1}}\) （我们可以将这个方程看作 \(g(x)=x-a\equiv 0\pmod{p^{r-1}}\)，其中 \(a\equiv (e_2-e_1)(e_1e_2)^{-1}\pmod{N}\)）

使用 \(m=8\) 和 \(t=6\)，我们构建一个维度为 \(\omega=9\) 的格。

应用 LLL 算法并接触第一个简化多项式，我们得到了解 \(x_0= 1826732340\)

因此 \(\gcd(f(x_0),N)=p=1789386140116417697\)，且最终 \(q=\frac{N}{p^2}=1903010275819064491\)

整个求解程序的运行时间使用现成的计算机是少于 4 秒的。

最后我们可以恢复私钥 \((d_1,d_2)\)。

但是要注意到，如果 \(d_1\approx d_2\approx N^{0.99}\)，且 Sarkar’s method 与存在范围 \(d<N^{0.395}\)，可能无法恢复私钥。

脚本

r = 7
N = 209798341155088334158217087474227805455138848036904381404809759100627849272231840321985747935471287990313456209656625928356468120896887536235496490078123448217785939608443507649096688546074968476040552137270080120417769906047001451239544719039212180059396791491281787790213953488743488306241516010351179070869410418232801398578982244984544906579574766534671056023774009163991804748763929626213884208260660722705479782932001102089367261720194650874553305179520889083170973755913964440175393646890791491057655226024046525748177999422035469428780228224800114202385209306803288475439775037067014297973202621118959024226798935588827359265962780792266516120013602384766460619793738405476219362508944225007365127768741191310079985425349292613888185378948854602285379329682053663283534930182589905986063348509703027498270111412063194971956202729807710253369312175636837558252924035002153389909587349043986253518050303628071319876207392440085675892353421232158925122721273720564784886530611286461575045181073744696415657043278123662980166364494583141297996445429477446442693717498789391918530672770193730629928408766563592081857706608049076318165712479742423149330311238462044666384622153280310696667586565906758451118241914402257039981388209

e1 = 167033559384298522723574512241709447697750495062051373339874928117760768631565225663704494711294488556402223152830158600944819657473430506318731286655519728589208977191031849602792050411662024383502548579402516538753112670329781366297260905517214408459895097308286783418547254449419676568096534767340832822470233461516097690657337287889405321592527860524201824371955082411119548743528220794151774190322092515459637969925138496615421690273925560390321721643556915400569894100488394008220811596560968566833296068500476868375996187754631888256419438775013308064639754700359028260289266420692324376220460340153811660590804281527733243177527178698523018103373311259548716062006020121615186595491453534952848570829485638553678760994354019044715078062414748269425818079274218450448217803229617020494546843594180682307375768323235309661628678546003718924902228908888185484412626429441196588985691713767554591991735686919964937441820738008046218954331990752603146125777571183543616375946363623251491371247594696767767918341279655251868517380267258878990871525012299220182939441091806206624720194246691865367941280852353547267930167542329486261552854451001546455904682702366584763940463481732752992487773878685793275652314513014646439770319249
e2 = 69076592619651589706691933313826601279528159013379300261609967352748175972662567592943146333144902972780621576811778115958019397062270814057821407036352529372113467206560849267275602453288227390740346959857322649956992529510338912182696854496200041245775322561359546062736323363354733510660780489558215103581313753430117471361013972291126160134685745917715386613876414886325025010348396410346222272648657374977901786530969589123771261040601627906959627351426881111464943086191212001374558078570830214670111422731410682212770683631011038163623234630601007231955235905528750031898751733232446644402069580930596404887288935724879795199659228145390574503341087565153744389617539607111733080406125228559950446336384154674927952991964565965760896308198785777527690939982523416579778957846249005801121682470447753074839399698315364445972142571727376297422736232659133510385808957304351692629177239808890209690661982628408419571131470406142532800330250274534615063537773403062635865734585850821677880659194795963303700015814615804751909674946908768425855361277478190640780518117596780808975720826484146074528564147729911624750510271539697935538038871993380673492022099183863825435237650082706168588306816635866830411481021066926833372846305

x = polygen(Zmod(N))
f = e1 * e2 * x - (e2 - e1)
f = f.monic()

# epsilon=0.05
roots = f.small_roots(beta=((r*(r-1)) / ((r+1)*(r+1))))
print(roots)

题目

from Crypto.Util.number import bytes_to_long, getPrime
from secret import msg
from sympy import nextprime
from gmpy2 import invert
from hashlib import md5

flag = 'd3ctf{'+md5(msg).hexdigest()+'}'
p = getPrime(256)
q = getPrime(256)
assert p > q
n = p * q
e = 0x10001
m = bytes_to_long(msg)
c = pow(m, e, n)

N = pow(p, 7) * q
phi = pow(p, 6) * (p - 1) * (q - 1)
d1 = getPrime(2000)
d2 = nextprime(d1 + getPrime(1000))
e1 = invert(d1, phi)
e2 = invert(d2, phi)

print(f'c = {c}')
print(f'N = {N}')
print(f'e1 = {e1}')
print(f'e2 = {e2}')
'''
c = 2420624631315473673388732074340410215657378096737020976722603529598864338532404224879219059105950005655100728361198499550862405660043591919681568611707967
N = 1476751427633071977599571983301151063258376731102955975364111147037204614220376883752032253407881568290520059515340434632858734689439268479399482315506043425541162646523388437842149125178447800616137044219916586942207838674001004007237861470176454543718752182312318068466051713087927370670177514666860822341380494154077020472814706123209865769048722380888175401791873273850281384147394075054950169002165357490796510950852631287689747360436384163758289159710264469722036320819123313773301072777844457895388797742631541101152819089150281489897683508400098693808473542212963868834485233858128220055727804326451310080791
e1 = 425735006018518321920113858371691046233291394270779139216531379266829453665704656868245884309574741300746121946724344532456337490492263690989727904837374279175606623404025598533405400677329916633307585813849635071097268989906426771864410852556381279117588496262787146588414873723983855041415476840445850171457530977221981125006107741100779529209163446405585696682186452013669643507275620439492021019544922913941472624874102604249376990616323884331293660116156782891935217575308895791623826306100692059131945495084654854521834016181452508329430102813663713333608459898915361745215871305547069325129687311358338082029
e2 = 1004512650658647383814190582513307789549094672255033373245432814519573537648997991452158231923692387604945039180687417026069655569594454408690445879849410118502279459189421806132654131287284719070037134752526923855821229397612868419416851456578505341237256609343187666849045678291935806441844686439591365338539029504178066823886051731466788474438373839803448380498800384597878814991008672054436093542513518012957106825842251155935855375353004898840663429274565622024673235081082222394015174831078190299524112112571718817712276118850981261489528540025810396786605197437842655180663611669918785635193552649262904644919
'''

exp 如下

from Crypto.Util.number import *
from hashlib import md5
import gmpy2

r = 7

N = 1476751427633071977599571983301151063258376731102955975364111147037204614220376883752032253407881568290520059515340434632858734689439268479399482315506043425541162646523388437842149125178447800616137044219916586942207838674001004007237861470176454543718752182312318068466051713087927370670177514666860822341380494154077020472814706123209865769048722380888175401791873273850281384147394075054950169002165357490796510950852631287689747360436384163758289159710264469722036320819123313773301072777844457895388797742631541101152819089150281489897683508400098693808473542212963868834485233858128220055727804326451310080791

e1 = 425735006018518321920113858371691046233291394270779139216531379266829453665704656868245884309574741300746121946724344532456337490492263690989727904837374279175606623404025598533405400677329916633307585813849635071097268989906426771864410852556381279117588496262787146588414873723983855041415476840445850171457530977221981125006107741100779529209163446405585696682186452013669643507275620439492021019544922913941472624874102604249376990616323884331293660116156782891935217575308895791623826306100692059131945495084654854521834016181452508329430102813663713333608459898915361745215871305547069325129687311358338082029
e2 = 1004512650658647383814190582513307789549094672255033373245432814519573537648997991452158231923692387604945039180687417026069655569594454408690445879849410118502279459189421806132654131287284719070037134752526923855821229397612868419416851456578505341237256609343187666849045678291935806441844686439591365338539029504178066823886051731466788474438373839803448380498800384597878814991008672054436093542513518012957106825842251155935855375353004898840663429274565622024673235081082222394015174831078190299524112112571718817712276118850981261489528540025810396786605197437842655180663611669918785635193552649262904644919

c = 2420624631315473673388732074340410215657378096737020976722603529598864338532404224879219059105950005655100728361198499550862405660043591919681568611707967

x = polygen(Zmod(N))
f = e1 * e2 * x - (e2 - e1)
f = f.monic()

roots = f.small_roots(beta=((r - 1) / (r + 1) - 0.01), epsilon=0.05)

res = gmpy2.iroot(gmpy2.gcd(int(f(roots[0])), N), r - 1)
assert res[1]
print("p =", res[0])
p = res[0]

q = N // (p ** r)
n = p * q
phi = (p - 1) * (q - 1)
d = inverse(0x10001, phi)
msg = long_to_bytes(power_mod(c, d, n))
print(msg)
flag = 'd3ctf{' + md5(msg).hexdigest() + '}'
print(flag)

The Third Attack

原理

设 \(N_1=p_1^rq_1,N_2=p_2^rq_2\)，其中 \(p_1>p_2\)，如果 \[ |p_1-p_2|<\frac{p_1}{2rq_1q_2} \] 那么我们可以在多项式时间内分解 \(N\)。（连分数定理）

类似于维纳攻击，结论是 \(\displaystyle\frac{q_2}{q_1}\) 是 \(\displaystyle\frac{N_2}{N_1}\) 的一个有理近似，这样同样的我们分解了 \(N_1,N_2\)： \[ p_1=\frac{N_1}{q_1}^{\frac{1}{r}}\\ p_2=\frac{N_2}{q_2}^{\frac{1}{r}} \]

脚本

N1=170987233913769420505896917437304719816691353833034482461
N2=120532911819726882881630714003135237766675602824250965921

def factorization(N1, N2):
    cf = continued_fraction(N2/N1)
    convers = cf.convergents()
    for pqq in convers:
        q2 = pqq.numerator()
        q1 = pqq.denominator()
        if q2 == 0:
            continue
        if q1 != 1 and N1 % q1 == 0:
            return q1, q2
    raise ValueError('Could not factor N!')

factorization(N1, N2)

参考文献

《New attacks on RSA with Moduli \(N = p^rq\)》