Summary of Crypto in CTF(RSA)

RSA - How to Use openssl

Given flag.enc, pubkey.pem/pub.key

1	openssl rsa -pubin -text -modulus -in warmup -in pubkey.pem

Then we get (e, n), after getting d:

from Crypto.Util.number import bytes_to_long, long_to_bytes

c = bytes_to_long(open('flag.enc','rb').read())
print(long_to_bytes(pow(c, d, n)))

RSA - Given (e, d, (n))

Theorem

[Initialize] Set $k\leftarrow d\cdot e-1$
[Try a random g] Choose random small primes g = 2, 3, 5, 7, …
[Next k] If $2\mid k,k\leftarrow\frac{k}{2}$ and $x\leftarrow g^k$ mod n, else go to step 2 (g $\leftarrow$ next_prime(g))
[Finish?] If $y\leftarrow gcd(x-1,n)\ mod\ n>1,p\leftarrow y,q\leftarrow n//y$, quit. Else go to step 3.

exp

from gmpy2 import next_prime, gcd


def Factorize(n, e, d):
    g = 2
    while True:
        k = e * d - 1
        while not k & 1:
            k //= 2
            p = int(gcd(pow(g, k, n) - 1, n)) % n
            if p > 1:
                return (p, n // p)
        g = int(next_prime(g))


if __name__ == "__main__":
    n = 
    e = 
    d = 
    print(Factorize(n, e, d))

Note

Given only (e, d) - [NCTF2019] babyRSA

RSA - Given (e, dp, n)

Theorem

$Known(e,dp,n),assume\ dp\equiv d(mod\ (p-1))$

$\because dp\cdot e\equiv d\cdot e\equiv 1(mod\ (p-1))$

$\therefore dp\cdot e-1=k\cdot (p-1)$

$\because dp\cdot e\equiv d\cdot e\equiv 1(mod\ (p-1))$

$\therefore dp\cdot e-1=k\cdot (p-1)$

$\therefore (dp\cdot e-1)\cdot d\cdot e=k’\cdot (p-1),\quad k’=k\cdot d\cdot e$

$\Leftrightarrow d*e=-k’\cdot (p-1)+dp\cdot e\cdot d\cdot e\equiv 1(mod\ \varphi(n))\Leftrightarrow -k’\cdot (p-1)+dp\cdot e\equiv 1(mod\ \varphi(n))$

$\therefore k_{1}\cdot (p-1)+dp\cdot e-1=k_{2}\cdot (p-1)\cdot (q-1)\Leftrightarrow (p-1)\cdot (k_{2}\cdot (q-1)-k_{1})+1=dp\cdot e$

$\because dp<p-1\quad \therefore (k_{2}\cdot (q-1)-k_{1})\in (0, e)$

$\therefore$遍历i in range(1, e)，当同时满足

$(dp\cdot e-1)\%i==0$和$n\%((dp\cdot e-1)//i+1)==0时$，N成功分解.

exp

from Crypto.Util.number import long_to_bytes
import gmpy2


def decrypt(e, dp, n):
    _ = dp * e - 1
    for i in range(1, e):
        if not _ % i:
            p = _ // i + 1
            if not n % p:
                q = n // p
                phi = (p - 1) * (q - 1)
                d = int(gmpy2.invert(e, phi))
                return d


if __name__ == "__main__":
    e = int(input("e = "))
    dp = int(input("dp = "))
    n = int(input("n = "))
    c = int(input("c = "))
    d = decrypt(e, dp, n)
    print(long_to_bytes(pow(c, d, n)))

RSA - Given (p, q, dp, dq)

Theorem

$Known(p,q,dp,dq),assume\ dp\equiv d(mod\ (p-1)),dq\equiv d(mod\ (q-1))$

$d\equiv dp(mod\ (p-1)),d\equiv dq(mod\ (q-1))\rightarrow ext_CRT\ (\because gcd(p-1,q-1)$不一定等于1)

$\because d=k_{1}(p-1)+dp=k_{2}(q-1)+dq\Leftrightarrow k_{1}(p-1)=(dq-dp)+k_{2}(q-1)$

$\Leftrightarrow k_{1}\frac{p-1}{gcd(p-1,q-1)}=\frac{dq-dp}{gcd(p-1,q-1)}+k_{2}\frac{q-1}{gcd(p-1,q-1)}$

$\Rightarrow k_{1}\frac{p-1}{gcd(p-1,q-1)}\equiv\frac{dq-dp}{gcd(p-1,q-1)}(mod\ \frac{q-1}{gcd(p-1,q-1)})$

$\Leftrightarrow k_{1}\equiv invert(\frac{p-1}{gcd(p-1,q-1)},\frac{q-1}{gcd(p-1,q-1)})\cdot \frac{dq-dp}{gcd(p-1,q-1)}(mod\ \frac{q-1}{gcd(p-1,q-1)})$

将$k_{1}=k_{3}\frac{q-1}{gcd(p-1,q-1)}+invert(\frac{p-1}{gcd(p-1,q-1)},\frac{q-1}{gcd(p-1,q-1)})\cdot \frac{dq-dp}{gcd(p-1,q-1)}$代入$d=k_{1}(p-1)+dp$

$d=k_{3}\frac{(p-1)(q-1)}{gcd(p-1,q-1)}+invert(\frac{p-1}{gcd(p-1,q-1)},\frac{q-1}{gcd(p-1,q-1)})\cdot \frac{(dq-dp)(p-1)}{gcd(p-1,q-1)}+dp$

$\Rightarrow d\equiv invert(\frac{p-1}{gcd(p-1,q-1)},\frac{q-1}{gcd(p-1,q-1)})\cdot \frac{(dq-dp)(p-1)}{gcd(p-1,q-1)}+dp(mod\ \frac{(p-1)(q-1)}{gcd(p-1,q-1)})$

exp

import gmpy2
from Crypto.Util.number import long_to_bytes


def decrypt(p, q, dp, dq):
    phi = (p - 1) * (q - 1)
    dd = int(gmpy2.gcd(p - 1, q - 1))
    modules = phi // dd
    d = int(gmpy2.invert((p - 1) // dd, (q - 1) // dd)) * \
        (dq - dp) * (p - 1) // dd + dp
    d %= modules
    while True:
        if gmpy2.gcd(d, phi) == 1:
            return d
        d += modules


if __name__ == "__main__":
    p = int(input("p = "))
    q = int(input("q = "))
    dp = int(input("dp = "))
    dq = int(input("dq = "))
    c = int(input("c = "))
    d = decrypt(p, q, dp, dq)
    plaintext = long_to_bytes(pow(c, d, p * q))
    print(plaintext)

sage

1	d = crt(dp, dq, p - 1, q - 1)

RSA - Wiener’s Attack

Theorem

$Let \ N=pq\ with\ q<p<2q. Let\ d<\frac{1}{3}N^{\frac{1}{4}}.\\\\ Known (N,e),we\ can\ use\ Wiener's\ attack.\\\\$

$Since\ ed\equiv 1(mod\ \varphi(N)),\exists k\in N,s.t.\ ed=k\varphi(N)+1\Rightarrow\left|\frac{e}{\varphi(N)}-\frac{k}{d}\right|=1.$
From q<p<2q, it follows $\sqrt{N}<p<\sqrt{2N}$, hence
$N-\varphi(N)=p+q-1=p+\frac{N}{p}-1<\frac{3\sqrt{2}}{2}\sqrt{N}<3\sqrt{N}$
$\therefore\left|\frac{e}{N}-\frac{k}{d}\right|=\left|\frac{ed-k\varphi(N)+k\varphi(N)-kN}{Nd}\right|=\left|\frac{1+k(\varphi(N)-N)}{Nd}\right|\leq\left|\frac{3k\sqrt{N}}{Nd}\right|=\left|\frac{3k}{\sqrt{N}d}\right|$
$\because k\varphi(N)=ed-1<ed,\ e<\varphi(N)\quad \therefore k\varphi(N)<d\varphi(N)\rightarrow k<d<\frac{1}{3}N^{\frac{1}{4}}$
$\therefore\left|\frac{e}{N}-\frac{k}{d}\right|\leq\left|\frac{3k}{\sqrt{N}d}\right|<\left|\frac{3k}{\sqrt{N}d}\right|<\frac{1}{3d\cdot d}<\frac{1}{2d^2}$

Therefore,we can get $\frac{k}{d}$ from the convergents of $\frac{e}{N}$.

exp

import gmpy2
from Crypto.Util.number import long_to_bytes


def rational_to_quotients(x, y):
    a = x // y
    quotients = [a]
    while a * y != x:
        x, y = y, x - a * y
        a = x // y
        quotients.append(a)
    return quotients


def convergents_from_quotients(quotients):
    convergents = [(quotients[0], 1)]
    for i in range(2, len(quotients) + 1):
        quotients_partion = quotients[0:i]
        denom = quotients_partion[-1]  # 分母
        num = 1
        for _ in range(-2, -len(quotients_partion), -1):
            num, denom = denom, quotients_partion[_] * denom + num
        num += denom * quotients_partion[0]
        convergents.append((num, denom))
    return convergents


def WienerAttack(e, n):
    quotients = rational_to_quotients(e, n)
    convergents = convergents_from_quotients(quotients)
    for (k, d) in convergents:
        if k and not (e * d - 1) % k:
            phi = (e * d - 1) // k
            # check if (x^2 - coef * x + n = 0) has integer roots
            coef = n - phi + 1
            delta = coef * coef - 4 * n
            if delta > 0 and gmpy2.iroot(delta, 2)[1] == True:
                print('d = ' + str(d))
                return d


if __name__ == "__main__":
    e = int(input("e = "))
    n = int(input("n = "))
    c = int(input("c = "))
    d = WienerAttack(e, n)
    plaintext = long_to_bytes(pow(c, d, n))
    print(plaintext)

RSA - Partial_m

Theorem

求解一元同余方程的问题可以这样来表达：给定一个未知分解的合数N，b是N的一个因子（注意b是合数），且对一个已知的$\beta(0<\beta<1)$，满足$b\geq N^{\beta}$。给定一个度数为$\delta$，首1的单变量多项式$f_{b}(x)\in Z[X]$，在$\delta$和N的比特长度的多项式函数时间内找到所有$b|f_{b}(x)$的满足$|x_{0}|\leq X$的根$x_{0}\in Z$，X是$x_{0}$的上界。目标是通过尽可能宽松的界X获得更大的解出x的可能性。

定理1 (Coppersmith 96)：$N,b,\beta,f,\delta$定义同上。此时我们可以在$O(c\delta^{5}log^{9}N)$内得到满足条件：

$|x_{0}|\leq \frac{1}{2}N^{\frac{\beta^{2}}{\delta}-\epsilon}$

的，一元同余方程f(x)=0 mod b的解 ($\epsilon$多取$\frac{\beta}{7}$, 证明略)

定理2 (CopperSmith 96[2])：定理1在N=b时的特殊情况，可以在(log N,$\delta$)的多项式时间内求出所有满足：

$|x_{0}|\leq N^{\frac{1}{\delta}}$

的模等式$f_{N}(x)=0\ mod\ N$的解

$\therefore$ e足够小，且部分明文泄露时，可以采用coppersmith单变量模等式的攻击，如下：

$c:=m^{e}\ mod\ n=(mbar+x_{0})^{e}\ mod\ n$ , (mbar = (m >> kbits) << kbits)

定理2知，$|x_{0}|\leq N^{\frac{1}{e}}$时，可以在logN和e的多项式时间内求出x

exp

n = 
e = 
c = 
mbar = 
kbits = 

beta = 1

nbits = n.nbits()
print("upper {} bits of {} bits is given".format(nbits - kbits, nbits))

PR.<x> = PolynomialRing(Zmod(n))
f = (mbar + x)^e - c

x0 = f.small_roots(X=2^kbits, beta=1)[0]  # find root < 2^kbits with factor = n
print("m:", mbar + x0)

Note

m高位泄露是coppersmith的原版结论，实际上可以扩展到高位/低位/高低位泄露的一般情况

$m=M+x_{0}\cdot 2^{k}+M’$

$f_{n}(x):=(M+x_{0}\cdot 2^{k}+M’)^{e}-c$

化为首1即可sage求解（只要未知明文$|x_{0}|\leq N^{\frac{1}{e}}$）

RSA - Partial_p

Theorem

定理3 (Coppersmith: kp高位比特)：$p>q,k\in N^{*},q\nmid k$，如果对于kp的一个估计值$\widetilde{p}$满足$|kp-\widetilde{p}|\leq 2N^{\frac{1}{4}}$

那么在多项式时间内可以得到N的分解

Proof：

$f_{p}(x):=x+\widetilde{p}$的一个解$x_{0}=(kp-\widetilde{p})\ mod\ p$ 满足$x_{0}\leq 2N^{\frac{1}{4}}$。该多项式首1，且$p\geq N^{\frac{1}{2}}$，由定理1知该方程有解

往一般形式上推广，发现p的低位比特泄露也能使用coppersmith攻击，化为首1即可，可类比partial_m的扩展形式

exp

n = 
e = 
c = 
pbar = 
kbits = 

print("upper %d bits (of %d bits) is given" % (pbar.nbits()-kbits, pbar.nbits()))

PR.<x> = PolynomialRing(Zmod(n))
f = x + pbar

x0 = f.small_roots(X=2^kbits, beta=0.4)[0]  # find root < 2^kbits with factor >= n^0.4
p = x0 + pbar
print("p:", p)
q = n // int(p)
d = inverse_mod(e, (p-1)*(q-1))
print("m:", pow(c, d, n))

Note

sage的small_root传参X不能过大，需自行判断阈值并调整（如果X过大，即使存在X内的解，也无法求出）

比如p的低位泄露时因为不确定缺失高位的具体比特数，所以要在$2^{\frac{n.nbits()}{2}-kbits}$附近作X的阈值估计

无法确定拿到的p是否大于q，所以对$\beta=0.5$进行调整至0.4

RSA - Partial_d

Theorem

Given the $\frac{n.nbits()}{4}$ least significant bits of d，we can compute all of d in polynomial time in n and e.

Proof：

$s \leftarrow p+q$，$ed=1+k\phi (n)=1+k(n-s+1)$

$\therefore ed_{0}\equiv 1+k(n-s+1)\ mod\ 2^{\frac{n}{4}}\quad (1)$

$\because p^{2}-sp+n\equiv 0\ mod\ 2^{\frac{n}{4}}\quad (2)$

$\therefore p\cdot(1),k\cdot(2)\rightarrow ed_{0}p\equiv p+kpn-kp^{2}-kn+kp\ mod\ 2^{\frac{n}{4}}\quad (2)$

$\therefore$解一元同余方程$ed_{0}x-kx(n-x+1)+kn\equiv x\ (mod\ 2^{d_{0}.nbits()})$即可得到p的部分低位$p_{0}$

问题转化为partial_p

exp

def partial_p(p0, kbits, n):
    PR.<x> = PolynomialRing(Zmod(n))
    nbits = n.nbits()
    f = 2^kbits*x + p0
    f = f.monic()
    roots = f.small_roots(X=2^(nbits//2-kbits), beta=0.4)  # find root < 2^(nbits//2-kbits) with factor >= n^0.4
    if roots:
        x0 = roots[0]
        p = gcd(2^kbits*x0 + p0, n)
        return ZZ(p)

    
def find_p(d0, kbits, e, n):
    X = var('X')
    for k in range(1, e+1):
        results = solve_mod([e*d0*X - k*X*(n-X+1) + k*n == X], 2^kbits)
        for x in results:
            p0 = ZZ(x[0])
            p = partial_p(p0, kbits, n)
            if p and p != 1:
                return p


if __name__ == '__main__':
    n = 
    e = 
    c = 
    d0 = 
    beta = 0.5
    
    nbits = n.nbits()
    kbits = d0.nbits()
    print("lower %d bits (of %d bits) is given" % (kbits, nbits))

    p = int(find_p(d0, kbits, e, n))
    print("found p: %d" % p)
    q = n//int(p)
    print("d:", inverse_mod(e, (p-1)*(q-1)))

Theorem

当padding长度$m\in (0,\lfloor\frac{n.nbits()}{e^{2}}\rfloor]$时，$g_{1}=x^{e}-c_{1},g_{2} = (x+y)^{e}-c_{2}$，可以先用结式求解上述二元高次方程组

定理4：如果$(x_{0},y_{0})$是方程组 $\left\{\begin{aligned}f(x,y)=0\\\\g(x,y)=0\end{aligned}\right.$ 的解，则$y_{0}$是$Res_{x}(f,g)$的一个根

sage里可以用resultant来实现结式计算，$g_{1}.resultant(g_{2})$默认将结式以第二个变量来表示，也可以指定$g_{1}.resultant(g_{2},y)$来将结式以x表示

求出y后，采用Related Message Attack，简要概括就是求解$gcd(g_{1},g_{2})$，如果结果是线性的，则攻击成功($x=-gcd(g_{1}, g_{2})[0]$)

exp

def short_pad_attack(c1, c2, e, n):
    PRxy.<x,y> = PolynomialRing(Zmod(n))
    PRx.<xn> = PolynomialRing(Zmod(n))
    PRZZ.<xz,yz> = PolynomialRing(Zmod(n))

    g1 = x^e - c1
    g2 = (x+y)^e - c2

    q1 = g1.change_ring(PRZZ)
    q2 = g2.change_ring(PRZZ)

    h = q2.resultant(q1)
    h = h.univariate_polynomial()
    h = h.change_ring(PRx).subs(y=xn)
    h = h.monic()

    kbits = n.nbits()//(2*e*e)
    diff = h.small_roots(X=2^kbits, beta=0.4)[0]  # find root < 2^kbits with factor >= n^0.4

    return diff

def related_message_attack(c1, c2, diff, e, n):
    PRx.<x> = PolynomialRing(Zmod(n))
    g1 = x^e - c1
    g2 = (x+diff)^e - c2

    def gcd(g1, g2):
        while g2:
            g1, g2 = g2, g1 % g2
        return g1.monic()

    return -gcd(g1, g2)[0]


if __name__ == '__main__':
    n = 
    e = 

    c1 =
    c2 = 

    diff = short_pad_attack(c1, c2, e, n)
    print("difference of two messages is %d" % diff)

    m1 = related_message_attack(c1, c2, diff, e, n)
    print("m1:", m1)
    print("m2:", m1 + diff)

Note

exp的change_ring缺失即报错，但同环下转换的具体原因，估计要翻sage手册orz

RSA - Boneh Durfee Attack

Theorem

比维纳攻击强度更高的小私钥攻击

引用ctf-wiki上的简略推导过程：

后续应该用到了coppersmith对于多元模等式的攻击（坑有点儿深，具体证明略）

exp

import time

############################################
# Config
##########################################

"""
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
    c = 0xe3505f41ec936cf6bd8ae344bfec85746dc7d87a5943b3a7136482dd7b980f68f52c887585d1c7ca099310c4da2f70d4d5345d3641428797030177da6cc0d41e7b28d0abce694157c611697df8d0add3d900c00f778ac3428f341f47ecc4d868c6c5de0724b0c3403296d84f26736aa66f7905d498fa1862ca59e97f8f866c

    # the hypothesis on the private exponent (the theoretical maximum is 0.292)
    delta = .28 # 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)
        print("plaintext:", pow(c, d, N))
    else:
        print("=== no solution was found ===")

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

if __name__ == "__main__":
    example()

RSA - $gcd(e,\varphi(n)\neq 1)$

Theorem

$gcd(e,\varphi(n))\neq 1$时，$e’\leftarrow e//gcd(e,\varphi(n)),d\leftarrow invert(e’,\varphi(n))$，$c^{d}\equiv m^{gcd(e,\varphi(n))}\ mod\ n$

已知$m^{gcd(e,\varphi(n))}$后攻击方向不定，视题目而定

Example

De1CTF2019 - Baby RSA

最后步骤的加密：

1 2	assert(c1==pow(flag,e1,pq1)) assert(c2==pow(flag,e2,pq2))

p, q1, q2 均为素数，$f_{1}=gcd(e_{1},\varphi(pq_{1}))\neq 1,f_{2}=gcd(e_{2},\varphi(pq_{2}))\neq 1$

但发现$gcd(e_{1},\varphi(pq_{1}))=gcd(e_{2},\varphi(pq_{2}))=14$，于是可以用ext_gcd求解

def GCRT(mi, ai):
    # mi,ai分别表示模数和取模后的值,都为列表结构
    curm, cura = mi[0], ai[0]
    for (m, a) in zip(mi[1:], ai[1:]):
        d = gmpy2.gcd(curm, m)
        c = a - cura
        assert (c % d == 0) #不成立则不存在解
        K = c // d * gmpy2.invert(curm // d, m // d)
        cura += curm * K
        curm = curm * m // d
        cura %= curm
    return (cura % curm, curm) #(解,最小公倍数)


f3, lcm = GCRT([n1, n2], [f1, f2])
assert(f3 % n1 == f1)
assert(f3 % n2 == f2)
assert(lcm == q1*q2*p)

又发现gcd=14的因子7来自p，所以令$n_{3}=q_{1}\cdot q_{2},c_{3}=f_{3}\%n_{3}$

1 2	d3 = gmpy2.invert(7, phi3) m3 = pow(c3, d3, n3)

$m_{3}\equiv m^{2}\ mod\ n_{3}$

作小公钥指数的爆破即可.

但m如果太大时，就要用到一些有限域上开二次根的算法，在本文后面会有所提及.

RSA - Polynomial based RSA

Theorem

设p为素数，$F_{p}$为p元有限域，$F_{p}$上多项式形式的RSA是基于$F_{p}$上多项式的一个已知定理，可类比欧拉定理

设$F_{p}[x]$是$F_{p}$上的一元多项式环，P和Q为$F_{p}(x)$上的既约多项式，$N=P\cdot Q$，$F_{p}[x]$的模N商环下的多项式M满足gcd(M,N)=1

则$M^{\varphi(N)}=1$（对于多项式下的$\varphi(N)$，不同于ZZ上RSA的$(p-1)\cdot (q-1)$，但基本定义相同，即度小于degree(N)的所有满足gcd(*,N)的多项式*）

$\because$ P, Q均为既约多项式，
$\therefore\varphi(P)=p^{P.degree()}-1,\varphi(Q)=p^{Q.degree()}-1,\varphi(N)=(p^{P.degree()}-1)\cdot (p^{Q.degree()}-1)$

加密过程中，明文m的多项式表示方法大致有以下两种：

对应ascii码作多项式系数

R.<y> = PolynomialRing(GF(p))

def gen_irreducable_poly(deg):
    while True:
        out = R.random_element(degree=deg)
        if out.is_irreducible():
            return out
       
P = gen_irreducable_poly(ZZ.random_element(length, 2*length))
Q = gen_irreducable_poly(ZZ.random_element(length, 2*length))

e = 65537

N = P*Q
S.<x> = R.quotient(N)

flag = 'xxxx'
flag = list(bytearray(flag.encode()))
m = S(flag)
c = m^e

在$GF(p^{nbits})$上fetch_int

P.<x>=PolynomialRing(GF(p))
## 略去P,Q,N生成代码
R.<y> = GF(p^2049)
m_int = Integer(bytes_to_long(flag))
m_poly = P(R.fetch_int(m_int))
c_poly = pow(m_poly, e, n)
c_int = R(c_poly).integer_representation()

exp

p = 
R.<y> = PolynomialRing(GF(p))
N = 
S.<x> = R.quotient(N)
c = 

p, q = N.factor()
p, q = p[0], q[0]
phi = (p**P.degree()-1)*(p**Q.degree()-1)
e = 
d = inverse_mod(e, phi)

m = c^d
print("".join([chr(c) for c in m.list()]))

RSA - Large Common Factor of p-1 and q-1

Theorem

Given $N=pq,g|(p-1),g|(q-1)$. (We call it common prime RSA)

If g is large enough, we have an efficient factoring method for N. (魔改rho)

[可行分析]

$p-1=ga,q-1=gb,N-1=gh(h=gab+a+b)$
$\therefore x^{N-1}\ (mod\ p)$至多有a个值
$\therefore$用$x\rightarrow x^{N-1}+3$代替rho中的迭代函数，能在$O(\sqrt a)=O(N^{1/4-\gamma /2})$下分解N.

exp

from Crypto.Util.number import *


def f(x, n):
    return (pow(x, n - 1, n) + 3) % n


def rho(n):
    i = 1
    while True:
        a = getRandomRange(2, n)
        b = f(a, n)
        j = 1
        while True:
            p = GCD(abs(a - b), n)
            print('{} in {} circle'.format(j, i))
            if p == n:
                break
            elif p > 1:
                return (p, n // p)
            else:
                a = f(a, n)
                b = f(f(b, n), n)
            j += 1
        i += 1


def main():
    n = 2062899536811871554818178359324161185631864322612928428737135410318396234838612759353304630466467882907238599097155276236974469077407927587176395102072905755064272779705304103155044175760772870357960042973800466612294120431806053764425995115279842108284227580763784948966673592802930476849233889221985236020626654664840038198318033843640935727802767372332931755961464416884826343968074190784107977317020369409515252779597705798495979144688260140280406347997234735504377161004933444363855269475346955181189441066456092467816293921375731674289071178645028262851698651731976240122083422857066357348367076175249101149547
    print(rho(n))


if __name__ == '__main__':
    main()

common prime RSA在g过小时，g也能轻易通过分解N-1来解出（rho等算法分解出的小因子尝试即可）

关于这类的题型我还没在CTF赛题中见到过，下面给出生成的脚本

from Crypto.Util.number import *


def gen_prime(nbits, gamma):
    g = 2 * getPrime(int(nbits * gamma))
    while True:
        a = getRandomNBitInteger(int((0.5 - gamma) * nbits - 1))
        p = g * a + 1
        if isPrime(p):
            b = getRandomNBitInteger(int((0.5 - gamma) * nbits - 1))
            q = g * b + 1
            while not isPrime(q) or GCD(a, b) != 1:
                b = getRandomNBitInteger(int((0.5 - gamma) * nbits - 1))
                q = g * b + 1
            return p, q


def gen_key(nbits, gamma):
    p, q = gen_prime(nbits, gamma)
    n = p * q
    lcm = (p * q) // GCD(p, q)
    e = getPrime(16)
    while GCD(e, lcm) != 1:
        e = getPrime(16)
    d = inverse(e, lcm)
    return (n, e), (p, q, d)


def main():
    (n, e), (p, q, d) = gen_key(2048, 0.485)
    print('n =', n)
    print('p =', p)
    print('q =', q)
    print('e =', e)


if __name__ == '__main__':
    main()

RSA - Small CRT Private Exponents

Theorem

平衡的RSA-CRT密码系统下，q<p<2q. 且满足$ed_{p}\equiv 1(mod\ p-1),ed_{q}\equiv 1(mod\ q-1)$

则$d_{p},d_{q}<min\\{\frac{1}{4}(\frac{N}{e})^{\frac{2}{5}},\frac{1}{3}N^{\frac{1}{4}}\\}$成立时，可有效破解RSA-CRT算法.

$\because ed_{p}=1+k(p-1),ed_{q}=1+l(q-1)\\\\ \therefore e^{2}d_{p}d_{q}+kl(1-N)+e(d_{p}(l-1)+d_{q}(k-1))=k+l-1.\\\\ 令x=d_{p}(l-1)+d_{q}(k-1),y=kl,z=k+l-1,\omega =d_{p}d_{q}\\\\ 有ex+(1-N)y+e^{2}\omega = z,矩阵表示下则为$

构造一组基如下：

我们有如下结论：

记$v_{0}=(4ex,4N^{\frac{1}{2}}y,e^{\frac{3}{5}}N^{\frac{2}{5}}z)$. 则可以看出$v_{0}$一定是L约化基的三个行向量的线性组合.

LLL计算后的约化基$\\{b_{1},b_{2},b_{3}\\}$，选取一个适当的上界M，令$|a_{i}|\leq M$, 则验证线性组合$v=a_{1}b_{1}+a_{2}b_{2}+a_{3}b_{3}$对应的$x,y,z,\omega$，看得到的数据能否正确分解N即可.

M视具体情况而定

exp

暂无（因为生成满足要求的$d_{p},d_{q}$以及使对应的$d>N^{0.292}$有点困难）

RSA - Parity Oracle

Theorem

假设存在一个oracle，能对给定密文进行解密并给出对应明文的奇偶信息，则我们只需要log(N)次就能解密任意密文.

$c\equiv m^{e}(mod\ N)$, 向oracle机发送$2^{e}c\equiv (2m)^{e}(mod\ N)$,则返回奇时，说明$2m>N$，反之$2m<N$

继续往下分析，此时发送$(2^{2})^{e}c(mod\ N)$, 如果先前返回奇，则此时等同于在校验$2(2m-N)$与N的大小关系，即返回奇说明$m\in (\frac{3N}{4},N)$，反之$m\in (\frac{N}{2},\frac{3N}{4})$，先前返回偶的这一轮次不再分析

接下来轮次也类似，二分log(N)次即可.

exp

import decimal
from pwn import *
from Crypto.Util.number import long_to_bytes


def oracle(c):
    io = remote('111.198.29.45', '42484')
    io.recvuntil('You can input ciphertext(hexdecimal) now\n')
    c = hex(c)[2:]
    if len(c) & 1:
        c = '0' + c
    io.sendline(c)
    res = io.recvline(keepends=False)
    io.close()
    assert(res == b'odd' or res == b'even')
    return res == b'odd'


def partial(c, e, n):
    nbits = n.bit_length()
    decimal.getcontext().prec = nbits
    low = decimal.Decimal(0)
    high = decimal.Decimal(n)
    for i in range(nbits):
        c = (c * pow(2, e, n)) % n
        if not oracle(c):
            high = (low + high) / 2
        else:
            low = (low + high) / 2
        print(i, '/', nbits)
    return int(high)


def main():
    #c = int(input("c = "))
    #e = int(input("e = "))
    #n = int(input("n = "))
    e = 0x10001
    n = 0x0b765daa79117afe1a77da7ff8122872bbcbddb322bb078fe0786dc40c9033fadd639adc48c3f2627fb7cb59bb0658707fe516967464439bdec2d6479fa3745f57c0a5ca255812f0884978b2a8aaeb750e0228cbe28a1e5a63bf0309b32a577eecea66f7610a9a4e720649129e9dc2115db9d4f34dc17f8b0806213c035e22f2c5054ae584b440def00afbccd458d020cae5fd1138be6507bc0b1a10da7e75def484c5fc1fcb13d11be691670cf38b487de9c4bde6c2c689be5adab08b486599b619a0790c0b2d70c9c461346966bcbae53c5007d0146fc520fa6e3106fbfc89905220778870a7119831c17f98628563ca020652d18d72203529a784ca73716db
    c = 0x4f377296a19b3a25078d614e1c92ff632d3e3ded772c4445b75e468a9405de05d15c77532964120ae11f8655b68a630607df0568a7439bc694486ae50b5c0c8507e5eecdea4654eeff3e75fb8396e505a36b0af40bd5011990663a7655b91c9e6ed2d770525e4698dec9455db17db38fa4b99b53438b9e09000187949327980ca903d0eef114afc42b771657ea5458a4cb399212e943d139b7ceb6d5721f546b75cd53d65e025f4df7eb8637152ecbb6725962c7f66b714556d754f41555c691a34a798515f1e2a69c129047cb29a9eef466c206a7f4dbc2cea1a46a39ad3349a7db56c1c997dc181b1afcb76fa1bbbf118a4ab5c515e274ab2250dba1872be0
    m = partial(c, e, n)
    print(long_to_bytes(m))


if __name__ == "__main__":
    main()

RSA - William’s p+1 and Pollard’s p-1

Theorem

如果要分解的N包含的素数因子p，使得p+1或p-1是平滑的（即只有适当上界的小因子），则William’s p+1和Pollard’s p-1算法能很好实现分解.

William’s p+1 Algorithm

Step 1：选择整数A(A>2)来表征卢卡斯序列$V_{j}=AV_{j-1}-V_{j-2}$，并初始化$V_{0}=2,V_{1}=A$.

Step 2：计算序列到$j=m!$.

ps：这步有很大的操作空间（在More里会涉及相关简化计算和证明）

Step 3：计算$gcd(V_{j}-2,N)$，如果结果等于1或N，则m++并返回Step2，否则进入Step4.

Step 4：检查p+1是否是m!的因子，如果是，则说明找到的p是p+1平滑的.

因为我们事先不知道p，所以也就无从计算雅各比符号(D/p)，但平方非剩余概率趋于1/2，考虑选取多个A多线程即可（n个A就能使互斥事件的概率降到大约$\frac{1}{2^{n}}$）
Pollard’s p-1 Algorithm

p-1是B-smooth时，很大概率上$(p-1)|B!$（不考虑p-1有很多重复的趋近于B的质因子）

所以计算$a^{m!}-1$，以及校验gcd即可，$m_{max}$可以是B也可以适当上调.

exp

import threading
from gmpy2 import iroot
from functools import reduce
from Crypto.Util.number import *
import time

factor_state = False


def pow_with_sqrt(a, c, b, n, mod):  # (a+c*sqrt(b))^n
    states = [(a % mod, c % mod)]
    for i in range(int(math.log(n, 2))):
        new_state = ((pow(states[i][0], 2, mod) + b * pow(states[i][1], 2, mod)) %
                     mod, (2 * states[i][0] * states[i][1]) % mod)
        states.append(new_state)
    n_bin = bin(n)[2:][::-1]
    res_x = 1
    res_y = 0
    for i in range(len(n_bin)):
        if n_bin[i] == '1':
            res_x, res_y = (res_x * states[i][0] + res_y * states[i][1] *
                            b) % mod, (res_x * states[i][1] + res_y * states[i][0]) % mod
    return res_x, res_y


def william_factor(n, index):
    start = time.clock()
    global factor_state
    A = getRandomRange(3, n)
    #m = 2
    #next_pos = 1
    m = 3000
    next_pos = reduce(lambda x, y : x * y, [i + 1 for i in range(m - 1)])
    B = A**2 - 4
    #res_x1, res_y1 = A, -1
    #res_x2, res_y2 = A, 1
    res_x1, res_y1 = pow_with_sqrt(A, -1, B, next_pos, n)
    res_x2, res_y2 = pow_with_sqrt(A, 1, B, next_pos, n)
    while True:
        if factor_state == True:
            return
        print((index, m))
        next_pos *= m #next_pos = m!
        C = inverse(pow(2, next_pos, n), n)
        res_x1, res_y1 = pow_with_sqrt(res_x1, res_y1, B, m, n)
        res_x2, res_y2 = pow_with_sqrt(res_x2, res_y2, B, m, n)
        res_x, res_y = (res_x1 + res_x2) % n, (res_y1 + res_y2) % n
        assert(iroot(B, 2)[1] == False and res_y != n)
        if iroot(B, 2)[1] == True:
            res_x = (res_x + res_y * iroot(B, 2)[0]) % n
        # Vi = C((A-sqrt(B))^(m!)+(A+sqrt(B))^(m!))
        Vi = (C * res_x) % n
        p = GCD(Vi - 2, n)
        assert(p != n) #p=n说明lucas序列下标过大
        if p != 1:
            factor_state = True
            print('p =', p)
            end = time.clock()
            print('cost {}s'.format(end - start))
            return
        m += 1


def main():
    n = 7941371739956577280160664419383740967516918938781306610817149744988379280561359039016508679365806108722198157199058807892703837558280678711420411242914059658055366348123106473335186505617418956630780649894945233345985279471106888635177256011468979083320605103256178446993230320443790240285158260236926519042413378204298514714890725325831769281505530787739922007367026883959544239568886349070557272869042275528961483412544495589811933856131557221673534170105409
    #n = 112729
    threads = []
    for i in range(3):
        t = threading.Thread(target=william_factor, args=(n, i))
        threads.append(t)
    for t in threads:
        t.start()
    for t in threads:
        t.join()


if __name__ == '__main__':
    main()

关于William’s p+1 Algorithm的证明：

Lucas Functions：设$\alpha,\beta$是$x^{2}-Px+Q=0$的根，且有

则定义$\Delta=(\alpha-\beta)^{2}=P^{2}-4Q$，知$P=\alpha + \beta,Q=\alpha\beta$.

有如下定理：

p为奇素数，且$p\nmid Q$，勒让德符号$(\frac{\Delta}{p})=\epsilon$，则有

$U_{(p-\epsilon)m}(P,Q)\equiv 0(mod\ p)\\\\ V_{(p-\epsilon)m}(P,Q)\equiv 2Q^{\frac{m(1-\epsilon)}{2}}(mod\ p)$

假设我们有p+1|R，Q=1，$(\frac{\Delta}{p})=\epsilon=-1$，则由上述定理可知：

$U_{R}(P,1)\equiv 0(mod\ p)\\\\ V_{R}(P,1)-2\equiv 0(mod\ p)$

因为易推出$V_{n+1}=PV_{n}-QV_{n-1}=PV_{n}-V_{n-1}$，则只要在该lucas序列下求出$V_{R}$即可.

我上面的exp里取的R是$m!(m\in N^{+})$，但m足够大时（如m=1000时，$R≈2^{8000}$），逐项递推变得低效不可解.

因此简化思想如下：

上述Lucas序列对应通项为$\frac{(P-\sqrt{P^{2}-4})^{n}+(P+\sqrt{P^{2}-4})^{n}}{2^{n}}$，先fuzz出一个m（m过大时，$gcd(V_R-2,N)=N$），再从$V_{m!}$开始p+1 factor即可，这里还应记住保留前状态$(P-\sqrt{P^{2}-4})^{(m-1)!}$，以便在求解$(P-\sqrt{P^{2}-4})^{m!}$时只要做前状态的m次幂。

这里也记录一下Coinc1dens师傅的出题笔记：

prime = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997]

def Lucas_pow(P, R):
    A, B = P, 2
    R = bin(R)[3:]
    for i in R:
        if int(i) == 1:
            A, B = (P * A**2 - A * B - P) % n, (A**2 - 2) % n
        else:
            A, B = (A**2 - 2) % n, (A * B - P) % n
    return gcd(A - 2, n)

def Williams_p_1():
    R = 1
    B = iroot(n, 2)[0]
    B = log(B)
    for pi in prime:
        for i in range(B//log(pi)):
            R *= pi
    while 1:
        P = random.randint(2, B)
        p = Lucas_pow(P, R)
        if p > 1 and p < n:
            return p

可以看到上述算法的R选取和先前选的m!不同，且采用$V_{f}\rightarrow V_{2f}$的跨度进行递推.

实验发现在p+1的小素数因子上限B已知时，采取上述方法会明显快一些（当然也可以求出该方法中的R，再通项直接求R，大约5s就能出）.

Link

https://math.boisestate.edu/~liljanab/BOISECRYPTFall09/Sundling.pdf

https://www.ams.org/journals/mcom/1982-39-159/S0025-5718-1982-0658227-7/S0025-5718-1982-0658227-7.pdf

RSA - Square Root in $F_{p}$

Theorem

Given $b^{2}\equiv a(mod\ p)$

Case 1：$p\equiv 3(mod\ 4)\rightarrow b:=\pm a^{\frac{p+1}{4}}(mod\ p)$
Case 2：Cipolla’s algorithm（一般性）
- Step 1：在$F_{p}$上找到一个t，使得$u:=t^{2}-a$是一个模p的平方非剩余（因为平方非剩余的概率趋于1/2，所以随机选取并用欧拉判别即可）
- Step 2：定义$F_{p^{2}}=F_{p}(\sqrt{t^{2}-a})=\\{x+y\sqrt{t^{2}-a}:x,y\in F_{p}\\}$，关于$F_{p^{2}}$是域的证明此处略. （当然，$\omega:=\sqrt{t^{2}-a}$在$F_{p}$上不存在，因此类似虚数存在）.先来证明两个引理：(1) $(x+y\omega)^{p}=\sum_{i=0}^{p}C_{p}^{i}x^{p-i}(y\omega)^{i}$，因为u是平方非剩余，所以$\omega^{p-1}=(\omega^{2})^{\frac{p-1}{2}}=u^{\frac{p-1}{2}}=-1$. 且p均整除二项式展开中间系数，上式可转化为$x^{p}+y^{p}\omega^{p}(费马小定理)=x+y\omega^{p}=x-y\omega$.；(2) $((t+\omega)^{\frac{p+1}{2}})^{2}=(t+\omega)^{p+1}=(t+\omega)(t-\omega)=t^{2}-\omega^{2}=a$.
  
  所以计算$(t+\omega)^{\frac{p+1}{2}}$即为二次根b.
Case 3：Atkin’s algorithm（$p\equiv 5(mod\ 8)$）

$c\leftarrow (2a)^{\frac{p-5}{8}}$

$i\leftarrow 2ac^{2}$

$b\leftarrow ac(i-1)$

exp

import math
from Crypto.Util.number import *


def euler_judge(x, p):
    if pow(x, (p - 1) // 2, p) == 1:
        return True
    else:
        return False


def pow_in_Fp2(a, b, n, p):  # (a+sqrt(b))^n in Fp^2
    states = [(a % p, 1)]
    for i in range(int(math.log(n, 2))):
        new_state = ((pow(states[i][0], 2, p) + b * pow(states[i][1], 2, p)) %
                     p, (2 * states[i][0] * states[i][1]) % p)
        states.append(new_state)
    n_bin = bin(n)[2:][::-1]
    res_x = 1
    res_y = 0
    for i in range(len(n_bin)):
        if n_bin[i] == '1':
            res_x, res_y = (res_x * states[i][0] + res_y * states[i][1]
                            * b) % p, (res_x * states[i][1] + res_y * states[i][0]) % p
    return res_x, res_y


def square_root(a, p):
    assert(euler_judge(a, p) == True)
    if p % 4 == 3:
        b = pow(a, (p + 1) // 4, p)
    elif p % 8 == 5:
        c = pow(2 * a, (p - 5) // 8, p)
        i = (2 * a * c**2) % p
        b = (a * c * (i - 1)) % p
    else:
        t = getRandomRange(1, p)
        u = t**2 - a
        while euler_judge(u, p):
            t = getRandomRange(1, p)
            u = t**2 - a
        b, check = pow_in_Fp2(t, u, (p + 1) // 2, p)
        assert(check == 0)
    return b, -b % p


def main():
    root = square_root(19, 65537)
    print(root)


if __name__ == "__main__":
    main()

有限域下的n-th root其实在sympy中也有函数处理

1
2
3

from sympy.ntheory.residue_ntheory import nthroot_mod

x = nthroot_mod(a, 2, p, all_roots=False) #求二次根

但p和n过大时，nthroot_mod在有限域上求解也会变得极其困难，因此下面再做一个该类型的拓展.

RSA - $e\mid (p-1),e\mid (q-1)$

Theorem

前面提到过$gcd(e,\varphi(n))\neq 1$的情况，但不针对$gcd(e,\varphi(n))=e$，于是这里就对$e\mid (p-1),e\mid (q-1)$的特殊情况进行讨论.

解题思路即求解m mod p和m mod q，再通过CRT还原m mod n. 主要难点则是在$GF(p)$上求e次根.

在有限域上求r-th root有两个常见算法（Adleman-Manders-Miller algorithm和Cipolla-Lehmer algorithm），但这里介绍的是Namhun Koo提出的一种更具一般性的开根算法，且在s足够小的时候更高效（$r^{s}\mid (p-1),r^{s}\nmid (p-1)$）.

New r-th root extraction formula over $F_{q}$ for $q\equiv lr^{s}+1(mod\ r^{s+1})$ and r is a prime

[Theorem 1]

设c是在$F_q$上的一个r次剩余，则$\exists b\in F_{q},s.t.\ ord_{q}(c^{r-1}b^{r})=r^{t}.(0\leq t<s)$

Proof：

$\because gcd(r,l)=1\quad \therefore\exists(\alpha,\beta),s.t.\ r\beta +r-1=l\alpha$.

$\alpha =\frac{r\beta +r-1}{l}$，令$\zeta=(c^{\alpha})^{\frac{q-1}{r^{s}}}=(c^{\alpha})^{\frac{q-1}{r^{s}}}c^{r\beta +r-1-l\alpha}=c^{r-1}(c^{\beta+\alpha\frac{q-1-lr^{s}}{r^{s+1}}})^{r}=c^{r-1}b^{r}$.

$\because c$是$F_{q}$上的r次剩余，且指数含$\frac{q-1}{r^{s}}$，则$ord_{q}(\zeta)=r^{t}.(t<s)$

则由Theorem 1，我们令$\xi$为$F_{q}$上的一$r^{s}$阶单位根，要取符合要求的$\xi$，我们可以使得$\xi=d^{\frac{q-1}{r^{s}}}$，其中d为模q的r次非剩余（概率约为$\frac{r-1}{r}$）.

则存在唯一的$(i,j)$，使得$\xi^{r^{s-t}}=\zeta^{i},\zeta=(\xi^{r^{s-t}})^j$，即$ij\equiv 1(mod\ r^{t})$.

[Theorem 2]

令$u\equiv j(r^{t}-1)r^{s-t-1}\equiv -jr^{s-t-1}(mod\ r^{s-1})$，则$F_{q}$下c的一个r次根即为$cb\xi^{u}$，b在Theorem 1中已给出定义.

Proof：

$(cb\xi^{u})^{r}=cc^{r-1}b^{r}\xi^{ur}=c\zeta\xi^{ur}$.

$\because u=j(r^{t}-1)r^{s-t-1}+kr^{s-1},\xi^{r^{s}}=1.$

$\therefore \zeta\xi^{ur}=(\xi^{r^{s-t}})^{j}\xi^{ur}=\xi^{jr^{s-t}+ur}=\xi^{jr^{s-t}+kr^{s}+j(r^{t}-1)r^{s-t}}=\xi^{jr^{s}}=1.$

$\therefore (cb\xi^{u})^{r}=c.$

[Remark 1]

Theorem里的$r\beta +r-1=l\alpha$确定的$\alpha$ mod r和$\beta$ mod l是唯一的. 且上式可以扩展到$\frac{q-1}{r}$的模域上（$r\beta + r-1\equiv l\alpha(mod\ \frac{q-1}{r})$），因为c是$F_{q}$上的r次剩余，即$c^{\frac{q-1}{r}}=1$.

[Remark 2]

$cb=c^{\beta+1+\alpha\frac{q-1-lr^{s}}{r^{s+1}}}=c^{\frac{(\beta+1)r^{s+1}+\alpha(q-1)-\alpha lr^{s}}{r^{s+1}}}=c^{\frac{r^{s}((\beta+1)r-\alpha l)+\alpha(q-1)}{r^{s+1}}}=c^{\frac{1+\frac{\alpha(q-1)}{r^{s}}}{r}}.$

[Remark 3]

Theorem 1结尾提到的$(i,j)$在$r^{t}$较大时求解困难，因此本方法在r和t过大时并不高效.

[利用]

s = 1时（即$q\equiv lr+1(mod\ r^{2})$时）

$\because t<s\quad\therefore t=0,u=0.$

$cb\xi^{u}=cb$即为一个满足要求的r次根.

由Remark 3知，$cb=c^{\frac{1+\frac{\alpha(q-1)}{r^{s}}}{r}}=c^{\frac{1+\frac{\alpha(q-1)}{r}}{r}}$，至于$\alpha$可以通过扩展欧几里得算法很轻易得到.
s > 1时，t可能落在{0, …, s - 1}上，因此分析时要分的情况会多不少（在这里就不占篇幅去写了，paper链接会在Link中给出）
上述算法给出的只是$G_{q}$上c的其中一个r次根，但要找出所有的r次根也很简单，只要借助r次单位根即可

$G_{q}$上的r次单位根($r\mid (q-1)$)为{$x^{i(\frac{q-1}{r})}$}($i\in[0,q)$), 我们取一个不为1的r次单位根与我们先前找到的特解生成一个r阶循环群即找到通解.

exp

import time
from binascii import unhexlify
from Crypto.Util.number import *


def rthroot(c, r, q):
    c %= q
    assert(isPrime(r) and (q - 1) % r == 0 and (q - 1) % (r**2) != 0)
    l = ((q - 1) % (r**2)) // r
    alpha = (-inverse(l, r)) % r
    root = pow(c, ((1 + alpha * (q - 1) // r) // r), q)
    return root


def allroot(r, q, root):
    all_root = set()
    all_root.add(root)
    while len(all_root) < r:
        new_root = root
        unity = pow(getRandomRange(2, q), (q - 1) // r, q)
        for i in range(r - 1):
            new_root = (new_root * unity) % q
            all_root.add(new_root)
    return all_root


'''
def crt(ai, mi):
    a1, m1 = ai[0], mi[0]
    a2, m2 = ai[1], mi[1]
    return (a1 * inverse(m2, m1) * m2 + a2 * inverse(m1, m2) * m1) % (m1 * m2)
'''


def decrypt(proot, qroot, p, q):
    count = 0
    total = len(proot) * len(qroot)
    step = total // 100
    t1 = inverse(q, p)
    t2 = inverse(p, q)
    for i in proot:
        for j in qroot:
            count += 1
            if count % step == 0:
                print('\r{}%'.format(count // step), end='')
            #root = (i, j)
            #m = crt(root, (p, q))
            m = (i * t1 * q + j * t2 * p) % (p * q)
            if unhexlify(hex(m)[2:10]) == b'NCTF':
                print('\n', long_to_bytes(m))
                return


def main():
    p = 199138677823743837339927520157607820029746574557746549094921488292877226509198315016018919385259781238148402833316033634968163276198999279327827901879426429664674358844084491830543271625147280950273934405879341438429171453002453838897458102128836690385604150324972907981960626767679153125735677417397078196059
    q = 112213695905472142415221444515326532320352429478341683352811183503269676555434601229013679319423878238944956830244386653674413411658696751173844443394608246716053086226910581400528167848306119179879115809778793093611381764939789057524575349501163689452810148280625226541609383166347879832134495444706697124741
    e = 0x1337
    c = 10562302690541901187975815594605242014385201583329309191736952454310803387032252007244962585846519762051885640856082157060593829013572592812958261432327975138581784360302599265408134332094134880789013207382277849503344042487389850373487656200657856862096900860792273206447552132458430989534820256156021128891296387414689693952047302604774923411425863612316726417214819110981605912408620996068520823370069362751149060142640529571400977787330956486849449005402750224992048562898004309319577192693315658275912449198365737965570035264841782399978307388920681068646219895287752359564029778568376881425070363592696751183359
    print('[+] Calculating e-th root...')
    start = time.clock()
    proot = rthroot(c, e, p)
    qroot = rthroot(c, e, q)
    end = time.clock()
    print('[*] Cost {}s'.format(end - start))
    print('[+] Calculating all e-th roots...')
    start = time.clock()
    all_proot = allroot(e, p, proot)
    all_qroot = allroot(e, q, qroot)
    end = time.clock()
    print('[*] Cost {}s'.format(end - start))
    print('[+] CRT cracking...')
    start = time.clock()
    decrypt(all_proot, all_qroot, p, q)
    end = time.clock()
    print('[*] Cost {}s'.format(end - start))


if __name__ == '__main__':
    main()

exp对应2019-NCTF-easyRSA.

做完以后去找了官方wp，发现预期解是Adleman-Manders-Miller algorithm，这里给出链接

http://www.soreatu.com/ctf/writeups/Writeup%20for%20Crypto%20problems%20in%20NCTF%202019.html#easyrsa

[注]

exp第一次调试的时候发现速度很慢-.-后来发现是大量crt的调用使得计算了接近$e^{2}$次的inverse(p, q)和inverse(q,p)，因此将这两个crt中用到的变量锁定即可.

Link

https://eprint.iacr.org/2013/117.pdf

https://arxiv.org/pdf/1111.4877.pdf

RSA - How to Use openssl

RSA - Given (e, d, (n))

Theorem

exp

Note

RSA - Given (e, dp, n)

Theorem

exp

RSA - Given (p, q, dp, dq)

Theorem

exp

sage

RSA - Wiener’s Attack

Theorem

exp

RSA - Partial_m

Theorem

exp

Note

RSA - Partial_p

Theorem

exp

Note

RSA - Partial_d

Theorem

exp

RSA - Short Pad Attack & Related Message Attack

Theorem

exp

Note

RSA - Boneh Durfee Attack

Theorem

exp

RSA - $gcd(e,\varphi(n)\neq 1)$

Theorem

Example

RSA - Polynomial based RSA

Theorem

exp

RSA - Large Common Factor of p-1 and q-1

Theorem

exp

More

RSA - Small CRT Private Exponents

Theorem

exp

RSA - Parity Oracle

Theorem

exp

RSA - William’s p+1 and Pollard’s p-1

Theorem

exp

More

Link

RSA - Square Root in $F_{p}$

Theorem

exp

More

RSA - $e\mid (p-1),e\mid (q-1)$

Theorem

exp

More

Link