一道“连分数+copperSmith+ECC“的题
By 天璇Merak——Rain
题目来源：X-NUCA‘2020’
题目
#!/usr/bin/env sage
from secret import FLAG
assert FLAG.startswith(b"X-NUCA{") and FLAG.endswith(b"}")
def key_gen(bits):
while True:
p = random_prime(2**bits)
q = random_prime(2**bits)
if p % 4 == 3 and q % 4 == 3:
break
if p < q:
p, q = q, p
N = p * q
while True:
x = getrandbits(bits // 2)
y = getrandbits(bits // 2)
if gcd(x, y) == 1 and (x * y) < (int(sqrt(2 * N)) // 12):
e = randint( int(((p + 1) * (q + 1) * 3 * (p + q) - (p - q) * int(N**0.21)) * y / (x * 3 * (p + q))), int(((p + 1) * (q + 1) * 3 * (p + q) + (p - q) * int(N**0.21)) * y / (x * 3 * (p + q))) )
if gcd(e, (p + 1) * (q + 1)) == 1:
k = inverse_mod(e, (p + 1) * (q + 1))
break
return (N, e, k)
if __name__ == "__main__":
bits = 1024
N, e, _ = key_gen(bits)
pt = (int.from_bytes(FLAG[:32], 'big'), int.from_bytes(FLAG[32:], 'big'))
ct = (0, 1)
d = (((pt[1])**2 - 1) * inverse_mod(((pt[1])**2 + 1) * (pt[0])**2, N)) % N
# 2000 years later...:)
for _ in range(e):
ct = ( int((ct[0] * pt[1] + ct[1] * pt[0]) * inverse_mod(1 + d * ct[0] * pt[0] * ct[1] * pt[1], N) % N), int((ct[1] * pt[1] + d * ct[0] * pt[0]) * inverse_mod(1 - d * ct[0] * pt[0] * ct[1] * pt[1], N) % N) )
f = open("output.txt", "wb")
f.write(str((e, N)).encode() + b'\n')
f.write(str(ct).encode())
f.close()
output.txt
(288969517294013178236187423377607850772706067194956328319540958788120421760563745859661120809993097599452236235703456953461446476016483100948287481999230043898368061651387268308645842547879026821842863879967704742559469599469159759360184157244907772315674219466971226019794131421405331578417729612598931842872757269134756215101980595515566901217084629217607502582265295755863799167702741408881294579819035951888562668951997777236828957162036234849207438819692480197365737237130918496390340939168630111890207700776894851839829623749822549994705192645373973493114436603297829506747411555800330860323339168875710029679, 6321130275268755691320586594611921079666212146561948694592313061609721619539590734495630218941969050343046016393977582794839173726817429324685098585960482266998399162720208269336303520478867387042992449850962809825380612709067651432344409349798118550026702892042869238047094344883994914342037831757447770321791092478847580639207346027164495372017699282907858775577530313354865815011726710796887715414931577176850854690237886239119894136091932619828539390021389626283175740389396541552356118540397518601098858527880603493380691706649684470530670258670128352699647582718206243920566184954440517665446820063779925391893)
(5899152272551058285195694254667877221970753694584926104666866605696215068207480540407327508300257676391022109169902014292744666257465490629821382573289737174334198164333033128913955350103258256280828114875165476209826215601196920761915628274301746678705023551051091500407363159529055081261677043206130866838451325794109635288399010815200512702451748093168790121961904783034526572263126354004237323724559882241164587153748688219172626902108911587291552030335170336301818195688699255375043513696525422124055880380071075595317183172843771015029292369558240259547938684717895057447152729328016698107789678823563841271755, 253027286530960212859400305369275200777004645361154014614791278682230897619117833798134983197915876185668102195590667437488411251835330785944874517235915807926715611143830896296709467978143690346677123639363900536537534596995622179904587739684155397043547262126131676366948937690378306959846311626889534352806134472610026603322329394769864728875293696851590640974817297099985799243285824842399573006841275494668451690794643886677303573329060084436896592291515021246248961538322485059619863786362159459122242131918702862396595818404578595841492379025543989260901540257216728185425462070297720884398220421012139424567)题目分析
密钥
n = p*q，p，q分别是1024bit左右的数
e是有特殊性质的数
k = inverse_mod(e, (p + 1) * (q + 1))，是不是有些眼熟：
RSA加密算法中的私钥 d = inverse_mod(e, (p-1)*(q-1))
所以猜测，k就是解密的私钥。
加密算法
回到主函数中，关注 “#2000 years later……）”后面的内容也有些眼熟：
for _ in range(e):
ct = ( int((ct[0] * pt[1] + ct[1] * pt[0]) * inverse_mod(1 + d * ct[0] * pt[0] * ct[1] * pt[1], N) % N), int((ct[1] * pt[1] + d * ct[0] * pt[0]) * inverse_mod(1 - d * ct[0] * pt[0] * ct[1] * pt[1], N) % N) )ECC的乘法运算代码大概就长这样，但是还需要判断具体是那种ECC算法。
把循环体中的式子用表达式写出来：
（ct[0], ct[1]）= (x1, y1) (pt[0], pt[1]) = (x1, x2)
好巧不巧，有个Edward ECC算法中定义的加法刚好是这样的（其实是Lord告诉我的），而且这个ECC算法的零点是（0，1），也就是题目中最开始的ct，Edward ECC的阶也刚好是
（p+1）（q+1），也就解释了为什么密钥生成中 k = inverse_mod(e, (p + 1) (q + 1))了。
解密思路
下面的乘法运算都是Edward ECC定义下的乘法运算
$$
因为\:k=e^{-1}\:mod\:(p+1)(q+1),\quad ct=ept
$$
$$
所以\:kct = ke*pt = pt
$$
因此，求出私钥k，计算 k*ct ，结果就是明文。
# ct = 密文，N，e为公钥，均在output文件里
d = (((ct[1]) ** 2 - 1) * inverse(((ct[1]) ** 2 + 1) * (ct[0]) ** 2, N)) % N
zero = (0, 1)
def Edward_add(a, b):
ct = a
pt = b
ct = (int((ct[0] * pt[1] + ct[1] * pt[0]) * inverse(1 + d * ct[0] * pt[0] * ct[1] * pt[1], N) % N),
int((ct[1] * pt[1] + d * ct[0] * pt[0]) * inverse(1 - d * ct[0] * pt[0] * ct[1] * pt[1], N) % N))
return ct
def Edward_mul(n, p):
r = zero
tmp = p
while 0 < n:
if n & 1 == 1:
r = Edward_add(r, tmp)
n, tmp = n >> 1, Edward_add(tmp, tmp)
return r
phi = (p+1)*(q+1)
k = invert(e, phi)
tmp = Edward_mul(k, ct) #tmp就是明文
print(long_to_bytes(tmp[0]), long_to_bytes(tmp[1]))分解N
所以又回到老问题：分解n，求出p，q
因为题目中的e有特殊的要求，所以以此为突破点，求解p，q
if p < q:
p, q = q, p
N = p * q
while True:
x = getrandbits(bits // 2)
y = getrandbits(bits // 2)
if gcd(x, y) == 1 and (x * y) < (int(sqrt(2 * N)) // 12):
e = randint( int(((p + 1) * (q + 1) * 3 * (p + q) - (p - q) * int(N**0.21)) * y / (x * 3 * (p + q))), int(((p + 1) * (q + 1) * 3 * (p + q) + (p - q) * int(N**0.21)) * y / (x * 3 * (p + q))) )把e写出来：
$$
\frac{[3(p+1)(q+1)(p+q)-(p-q)N^{0.21}]y}{3x(p+q)}<e<\frac{[3(p+1)(q+1)(p+q)+(p-q)N^{0.21}]y}{3x(p+q)}
$$
$$
其中\quad (p+1)*(q+1)=\phi(n)\quad (Edward\:ECC下的\phi)
$$
$$
\Rightarrow |ex-y\phi (n)|<\frac{(p-q)n^{0.21}y}{3(p+q)}
$$
有了上面的结论，加上q < p < 2q，以及 x*y 不会很大，就可以判断出：
$$
\frac{y}{x}存在于\frac {e}{N}的连分数序列中
$$
通过遍历e/N的连分数序列，就可以确定x，y，然后：
$$
看下e的取值下限：\quad \frac{3y(p+1)(q+1)(p+q)}{3x(p+q)}-\frac{(p-q)N^{0.21}y}{3x*(p+q)}
$$
$$
前一项的位数是2048bit（分子分母的x、y、(p+q)约掉了），后一项是430bit（N^{0.21}）
$$
$$
所以e的位数是2048bit
$$
$$
观察式子\quad |ex-y\phi (n)|<\frac{(p-q)n^{0.21}y}{3(p+q)}
$$
$$
左边的位数是1048+512=2560bit，右边是430+512=942bit
$$
$$
\Rightarrow \:ex\approx y\phi(n)
$$
$$
ex\approx y*(p+1)(q+1)
$$
$$
ex-y(pq+1)\approx y(p+q)
$$
$$
[ex-y(n+1)]//y\approx (p+q)
$$
就可以计算出 (p+q) 的大概值，然后：
$$
(p+q)^2-4pq=(p-q)^2
$$
$$
[\sqrt {(p-q)^2}+(p+q)]//2\approx p
$$
这样就可以得到p的高位了
然后确定一下p的低位中有哪些位是不确定的：
从前面的推导有：
$$
e*x-y\phi (n)=r,\quad r是942bit
$$
$$
\frac{e*x}{y}-(n+1)=(p+q)+\frac{r}{y}
$$
$$
\frac{r}{y}是942-512 = 430bit
$$
所以p+q的低430bit是不确定的。
def rational_to_quotients(x, y):
"""
求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):
"""求连分数序列的第k个收敛子，convergent元组中的第二个数就是k
"""
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 x_y(e, n): # 连分数， 返回 p+q 的大概值, 高位是确定的
quotients = rational_to_quotients(e, n)
convergents = convergents_from_quotients(quotients)
#print(convergents)
for (y, x) in convergents:
if(y > 0):
if (1023 < int((e*x-y*(n+1))//y).bit_length() < 1026) :
print(x,y,sep=' ***\n')
return (e*x-y*(n+1))//y
pq = int(x_y(e, N))
p_h = (pq +int(iroot((pq**2)-4*N, 2)[0])) // 2 #得到p的高位是确定的 Copper Smith求解p、q：
前面退出p+q的低430bit不确定，但这也是估计值，为了准确，这里采取440bit：
pbar = 132160284144608950019816194803720605665582054407890340625286428343034451279699999656554400403442321672129341860427814515935184696844617907072796285688260865300923112869612920717393389100884250652210157065314926375313572891619847585468382475525427237283595480199209171431928054228220016002142262735297735560971
kbits = 440
# pbar p的大概值，高位已知； kbits: p未确定的位数值
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)
print("q: ", q)