ECDLP 考点

ECDLP (Elliptic Curve Discrete Logarithm Problem)，即椭圆曲线上的离散对数问题。

1.Edwards Curves DLP

题目代码:

题目来源：Crypto CTF 2021 Rohald

#!/usr/bin/env sage

from Crypto.Util.number import *
from secret import flag, Curve

def ison(C, P):
	c, d, p = C
	u, v = P
	return (u**2 + v**2 - c**2 * (1 + d * u**2*v**2)) % p == 0

def teal(C, P, Q):
	c, d, p = C
	u1, v1 = P
	u2, v2 = Q
	assert ison(C, P) and ison(C, Q)
	u3 = (u1 * v2 + v1 * u2) * inverse(c * (1 + d * u1 * u2 * v1 * v2), p) % p
	v3 = (v1 * v2 - u1 * u2) * inverse(c * (1 - d * u1 * u2 * v1 * v2), p) % p
	return (int(u3), int(v3))

def peam(C, P, m):
	assert ison(C, P)
	c, d, p = C
	B = bin(m)[2:]
	l = len(B)
	u, v = P
	PP = (-u, v)
	O = teal(C, P, PP)
	Q = O
	if m == 0:
		return O
	elif m == 1:
		return P
	else:
		for _ in range(l-1):
			P = teal(C, P, P)
		m = m - 2**(l-1)
		Q, P = P, (u, v)
		return teal(C, Q, peam(C, P, m))

c, d, p = Curve

flag = flag.lstrip(b'CCTF{').rstrip(b'}')
l = len(flag)
lflag, rflag = flag[:l // 2], flag[l // 2:]

s, t = bytes_to_long(lflag), bytes_to_long(rflag)
assert s < p and t < p

P = (398011447251267732058427934569710020713094, 548950454294712661054528329798266699762662)
Q = (139255151342889674616838168412769112246165, 649791718379009629228240558980851356197207)

print(f'ison(C, P) = {ison(Curve, P)}')
print(f'ison(C, Q) = {ison(Curve, Q)}')

print(f'P = {P}')
print(f'Q = {Q}')

print(f's * P = {peam(Curve, P, s)}')
print(f't * Q = {peam(Curve, Q, t)}')

# ison(C, P) = True
# ison(C, Q) = True
# P = (398011447251267732058427934569710020713094, 548950454294712661054528329798266699762662)
# Q = (139255151342889674616838168412769112246165, 649791718379009629228240558980851356197207)
# s * P = (730393937659426993430595540476247076383331, 461597565155009635099537158476419433012710)
# t * Q = (500532897653416664117493978883484252869079, 620853965501593867437705135137758828401933)

题目分析:

由下面代码可知曲线方程为 Edwards Curves ：$u^2+v^2=c^2(1+d{u}^2{v}^2)(modp)$ ，在有限域 $p$ 上的。

def ison(C, P):
	c, d, p = C
	u, v = P
	return (u**2 + v**2 - c**2 * (1 + d * u**2*v**2)) % p == 0

已知量有曲线上的四个点 $P,Q,sP,tQ$ ；未知量有 $c,d,p$ ，已知公式：$sP=s\ast P，tQ=t\ast Q$

我们需要计算出 $s,t$ 的值即得到 flag 值了，这里就是解椭圆曲线上的离散对数问题。

• Step1. 求 p值：对满足曲线上的任意点$(x,y)$ 有

$x^2+y^2\equiv c^2(1+d{x}^2{y}^2)(modp)=>x^2+y^2-c^2(1+d{x}^2{y}^2)=kp$

所以对于曲线上的 $P,Q$ 两点有：

$X_1=x_1^2+y_1^2-c^2(1+d{x}_1^2{y}_1^2)=k_{1}p$

$X_2=x_2^2+y_2^2-c^2(1+d{x}_2^2{y}_2^2)=k_{2}p$

两式相减：

$X_1-X_2=(x_1^2-x_2^2+y_1^2-y_2^2)-c^{2}d(x_1^2{y}_1^2-x_2^2{y}_2^2)=(k_1-k_2)p$

这里由于 $c,d$ 的值未知，所以要想办法消去 $c,d$:

令：$A_1=(x_1^2-x_2^2+y_1^2-y_2^2)，A_2=(x_1^2{y}_1^2-x_2^2{y}_2^2)$

则：$\frac{A_1}{A_2}\equiv c^{2}dmodp$

同理对于$sP,tQ$ 两点有：

令：$A_3=(x_3^2-x_4^2+y_3^2-y_4^2)，A_2=(x_3^2{y}_3^2-x_4^2{y}_4^2)$

则：$\frac{A_3}{A_4}\equiv c^{2}dmodp$

联立上面两式有：$\frac{A_1}{A_2}-\frac{A_3}{A_4}=k_{3}p$ ， $A_1{A}_4-A_2{A}_3=A_2{A}_4{k}_{3}p$

同理：对于曲线上的 $P(x_1,y_1),tQ(x_4,y_4)$ 两点和$Q(x_2,y_2),sP(x_3,y_3)$ 两点可以构建等式

令：$A_5=(x_1^2-x_4^2+y_1^2-y_4^2)，A_6=(x_1^2{y}_1^2-x_4^2{y}_4^2)$

$A_7=(x_2^2-x_3^2+y_2^2-y_3^2)，A_8=(x_2^2{y}_2^2-x_3^2{y}_3^2)$

有：$\frac{A_5}{A_6}-\frac{A_7}{A_8}=k_{4}p$，$A_5{A}_8-A_6{A}_7=A_6{A}_8{k}_{3}p$

然后求两式的最大公因数即可得到 $p$ 值 $p=gcd(A_1{A}_4-A_2{A}_3,A_5{A}_8-A_6{A}_7)$。

from gmpy2 import gcd

P = (398011447251267732058427934569710020713094, 548950454294712661054528329798266699762662)
Q = (139255151342889674616838168412769112246165, 649791718379009629228240558980851356197207)
sP = (730393937659426993430595540476247076383331, 461597565155009635099537158476419433012710)
tQ = (500532897653416664117493978883484252869079, 620853965501593867437705135137758828401933)

# 获取 p 的倍数
def get_P(P, Q):
    x1, y1 = P
    x2, y2 = Q
    A1 = (x1 ** 2 - x2 ** 2 + y1 ** 2 - y2 ** 2)
    A2 = (x1 ** 2 * y1 ** 2 - x2 ** 2 * y2 ** 2)
    return A1, A2

A1, A2 = get_P(P, Q)
A3, A4 = get_P(sP, tQ)
A5, A6 = get_P(P, tQ)
A7, A8 = get_P(Q, sP)

p = gcd(A1 * A4 - A2 * A3, A5 * A8 - A6 * A7)
print(p)

# 18983346087633426019400112058348303476269889

因为椭圆曲线使用的有限域为素数，这里检验得到的 p 值并非素数，我们使用 sage 的 ecm.factor() 函数对 p 进行因式分解，进而得到精确的p值。

ecm.factor(18983346087633426019400112058348303476269889)
Out: [3, 7, 903968861315877429495243431349919213155709]

is_prime(903968861315877429495243431349919213155709)
True

• Step2. 求 $c,d$ 的值。

1、$c^{2}d=\frac{(x_1^2-x_2^2+y_1^2-y_2^2)}{(x_1^2{y}_1^2-x_2^2{y}_2^2)}(modp)$

2、由曲线方程：$x^2+y^2=c^2(1+d{x}^2{y}^2)(modp)$ 可得，$c^2=x^2+y^2-c^{2}d{x}^2{y}^2(modp)$

3、求得 $c$ 的值后，由：$d=\frac{c^{2}d}{c^2}$ 即可得到 d 的值。

实现代码如下：

X1, X2 = get_P(P, Q)
ccd = X1 * pow(X2, -1, p) % p
print(ccd)

x, y = P
# cc = c^2
cc = (x ** 2 + y ** 2 - ccd * x ** 2 * y ** 2) % p
d = ccd * pow(cc, -1, p)

# 校验 c,d,p 的值是否满足曲线
def ison(cc, d, P):
    u, v = P
    return (u ** 2 + v ** 2 - cc * (1 + d * u ** 2 * v ** 2)) % p == 0

assert ison(cc, d, P) == True
print("ok")

• Step3. 因为我们用 sage 实现椭圆曲线的函数 EllipticCurve() 是Weierstrass curves 形式，所以这里我们需要先将 Edwards Curves 转化为 Weierstrass curves ，然后再在该曲线上求解离散对数问题。

根据 Edwards Curves(爱德华兹曲线) 理论知识可知：

这里曲线的形式为：$E_{c,d}：x^2+y^2=c^2(1+d{x}^2{y}^2)$

• 将 $E_{c,d}：x^2+y^2=c^2(1+d{x}^2{y}^2)$ 映射为： $E_d：x^2+y^2=1+d{x}^2{y}^2$

相关论文 Faster addition and doubling on elliptic curves 理论在第 5 页

映射的实现基于有限域p，令 $d=d_1\ast c_1^4$，$x=\frac{x_1}{c_1}，y=\frac{y_1}{c_1}$

则有：$x_1^2+y_1^2=c_1^2(1+d_1{x}_1^2{y}_1^2)⟼x^2+y^2=1+d{x}^2{y}^2$

• 映射为 Weierstrass 曲线：$E_d：x^2+y^2=1+d{x}^2{y}^2⟼E_{a2,a4}:v^2=u^3+2(d+1)u^2+(d-1{)}^{2}u)$

$$(x,y)⟼(u,v)=(\frac{A}{x^2},\frac{-2A}{x^3})$$

其中： $A=2y-(2dy+d+1)x^2+2$

相关论文 Edwards Elliptic Curves 理论在第 22—23 页

实现代码：

from gmpy2 import gcd
from sympy import isprime
from Crypto.Util.number import *

P = (398011447251267732058427934569710020713094, 548950454294712661054528329798266699762662)
Q = (139255151342889674616838168412769112246165, 649791718379009629228240558980851356197207)
sP = (730393937659426993430595540476247076383331, 461597565155009635099537158476419433012710)
tQ = (500532897653416664117493978883484252869079, 620853965501593867437705135137758828401933)


# 获取 p 的倍数
def get_P(P, Q):
    x1, y1 = P
    x2, y2 = Q
    A1 = (x1 ** 2 - x2 ** 2 + y1 ** 2 - y2 ** 2)
    A2 = (x1 ** 2 * y1 ** 2 - x2 ** 2 * y2 ** 2)
    return A1, A2


A1, A2 = get_P(P, Q)
A3, A4 = get_P(sP, tQ)
A5, A6 = get_P(P, tQ)
A7, A8 = get_P(Q, sP)

p = gcd(A1 * A4 - A2 * A3, A5 * A8 - A6 * A7)
p = ecm.factor(p)[2]
print("p =", p)
assert is_prime(p)

X1, X2 = get_P(P, Q)
ccd = X1 * pow(X2, -1, p) % p
print("ccd =", ccd)

x, y = P
# cc = c^2
cc = (x ** 2 + y ** 2 - ccd * x ** 2 * y ** 2) % p
d = ccd * pow(cc, -1, p)
print("d =", d)

# 对 c的平方进行有限域上的开平方根
# c = GF(p)(cc).sqrt()
c = cc.nth_root(2)
print("c =", c)


# 校验是否满足爱德华兹曲线方程
def ison(c, d, p, P):
    u, v = P
    return (u ** 2 + v ** 2 - c ** 2 * (1 + d * u ** 2 * v ** 2)) % p == 0


assert ison(c, d, p, P) == True

# 转化为 Edwards 曲线的一般形式
d = GF(p)(d) * GF(p)(c) ** 4


def phi1(P):
    return (P[0] / c, P[1] / c)


P = phi1(P)
Q = phi1(Q)
sP = phi1(sP)
tQ = phi1(tQ)

assert (P[0] ^ 2 + P[1] ^ 2 - (1 + d * P[0] ^ 2 * P[1] ^ 2)) % p == 0
assert (Q[0] ^ 2 + Q[1] ^ 2 - (1 + d * Q[0] ^ 2 * Q[1] ^ 2)) % p == 0

# 映射为 Weierstrass 曲线的一般形式
a2 = 2 * (d + 1)
a4 = (d - 1) ** 2
E = EllipticCurve(GF(p), [0, a2, 0, a4, 0])


def phi2(P):
    x, y = P
    A = 2 * y - (2 * d * y + d + 1) * x ** 2 + 2
    u = A / (x ** 2)
    v = (-2 * A) / (x ** 3)
    return u, v


P = E(phi2(P))
Q = E(phi2(Q))
sP = E(phi2(sP))
tQ = E(phi2(tQ))

# 计算ECDLP离散对数
s = discrete_log(sP, P, operation="+")
t = discrete_log(tQ, Q, operation="+")
print("s =", s)
print("t =", t)
flag = long_to_bytes(s) + long_to_bytes(t)
print(flag)

flag{nOt_50_3a5Y_Edw4rDs_3LlipT!c_CURv3}

2.Twisted Edwards Curves DLP

题目代码:

题目来源：NCTF 2022 superecc

from Crypto.Util.number import *
from secrets import INF, flag

assert flag[:5] == b'nctf{'


class super_ecc:
    def __init__(self):
        self.a = 73101304688827564515346974949973801514688319206271902046500036921488731301311
        self.c = 78293161515104296317366169782119919020288033620228629011270781387408756505563
        self.d = 37207943854782934242920295594440274620695938273948375125575487686242348905415
        self.p = 101194790049284589034264952247851014979689350430642214419992564316981817280629

    def add(self, P, Q):
        (x1, y1) = P
        (x2, y2) = Q
        x3 = (x1 * y2 + y1 * x2) * inverse(self.c * (1 + self.d * x1 * x2 * y1 * y2), self.p) % self.p
        y3 = (y1 * y2 - self.a * x1 * x2) * inverse(self.c * (1 - self.d * x1 * x2 * y1 * y2), self.p) % self.p
        return (x3, y3)

    def mul(self, x, P):
        Q = INF
        x = x % self.p
        while x > 0:
            if x % 2 == 1:
                Q = self.add(Q, P)
            P = self.add(P, P)
            x = x >> 1
        return Q


flag = bytes_to_long(flag[5:-1])
ecc = super_ecc()
G = (30539694658216287049186009602647603628954716157157860526895528661673536165645,
     64972626416868540980868991814580825204126662282378873382506584276702563849986)
S = ecc.mul(flag, G)
print(S)
# (98194560294138607903211673286210561363390596541458961277934545796708736630623,
# 58504021112693314176230785309962217759011765879155504422231569879170659690008)

题目分析:

def add(self, P, Q):
    (x1, y1) = P
    (x2, y2) = Q
    x3 = (x1 * y2 + y1 * x2) * inverse(self.c * (1 + self.d * x1 * x2 * y1 * y2), self.p) % self.p
    y3 = (y1 * y2 - self.a * x1 * x2) * inverse(self.c * (1 - self.d * x1 * x2 * y1 * y2), self.p) % self.p
    return (x3, y3)

分析加密代码可知曲线上的点加为：

$$(x_1,y_1)+(x_2,y_2)=(\frac{x_1{y}_2+x_2{y}_1}{c(1+d{x}_1{x}_2{y}_1{y}_2)},\frac{y_1{y}_2-a{x}_1{x}_2}{c(1-d{x}_1{x}_2{y}_1{y}_2)})$$

题目给的已知量有：曲线信息 $a,c,d,p$ 和基点 $G$ ，考察的是扭曲的爱德华兹曲线，曲线方程为：$E_{a,c,d}：a{x}^2+y^2=c^2(1+d{x}^2{y}^2)$

等价变形可得：$a(\frac{x}{c}{)}^2+(\frac{y}{c}{)}^2=1+d{c}^4(\frac{x}{c}{)}^2(\frac{y}{c}{)}^2$

令 $x_1=\frac{x}{c}，y_1=\frac{y}{c},d_1=d{c}^4$

则原曲线方程可以映射为一般形式： $a{x}_1^2+y_1^2=1+d{x}_1^2{y}_1^2$

根据论文 Twisted Edwards Curves 理论在第 4 页：

• 映射为 Montgomery 曲线

$E_{a,d}⟼M_{A,B}：B{v}^2=u^3+A{u}^2+u$

其中：

$A=\frac{2(a+d)}{a-d}，B=\frac{4}{a-d}$

$$u=\frac{1+y}{1-y}，v=\frac{1+y}{x(1-y)}=\frac{u}{x}$$

• 再将Montgomery 曲线映射为 Weierstrass 曲线即可实现在 sage 上解离散对数问题。

$M_{A,B}⟼E_{a,b}:v^2=u^3+au+b$

其中：$(x,y)⟼(u,v)=(\frac{x}{B}+\frac{A}{3B},\frac{y}{B})$

$a=\frac{3-A^2}{3B^2}$ ，$b=\frac{2A^3-9A}{27B^3}$

然后根据代码：

S = ecc.mul(flag, G)

可知，解有限域上的离散对数问题即可得到 flag。

这里解离散对数我们发现直接使用 discrete_log(S,G) 函数并不能直接得到 flag值，分析椭圆曲线的基点 G 的阶，发现这里 G.order() 是一个光滑数，则可以利用Pohlig-Hellman算法 实现降低离散对数计算中阶的大小，进而降低计算复杂度。

n = G.order()
print(n)
print(ecm.factor(n))

# 101194790049284589034264952247851014979635478581961323247248565451562534055248
# [2, 2, 2, 2, 113, 28921, 749971, 87374629, 20375070355228635421, 1449497564615145263845097285948144399]

曲线上的 Pohlig-Hellman 算法实现文章：Pohlig-Hellman 算法(针对阶是光滑且仅有小素因子)

实现代码:

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


G = (30539694658216287049186009602647603628954716157157860526895528661673536165645,
     64972626416868540980868991814580825204126662282378873382506584276702563849986)
a = 73101304688827564515346974949973801514688319206271902046500036921488731301311
c = 78293161515104296317366169782119919020288033620228629011270781387408756505563
d = 37207943854782934242920295594440274620695938273948375125575487686242348905415
p = 101194790049284589034264952247851014979689350430642214419992564316981817280629

assert (a * G[0] ** 2 + G[1] ** 2) % p == c ** 2 * (1 + d * G[0] ** 2 * G[1] ** 2) % p
print("符合曲线方程 ax^2+y^2 = c^2(1+dx^2y^2)")


def TEdwards_to_TEdwards(G, a, c, d):
    x, y = G
    x1 = x * pow(c, -1, p)
    y1 = y * pow(c, -1, p)
    d = d * c ** 4
    G = (x1, y1)
    return (G, a, d)


def TEdwards_to_Montgomery(G, a, d):
    x, y = G
    A = 2 * (a + d) * pow(a - d, -1, p)
    B = 4 * pow(a - d, -1, p)
    u = (1 + y) * pow(1 - y, -1, p)
    v = u * pow(x, -1, p)
    G = (u, v)
    return (G, A, B)


def Montgomery_to_Weierstrass(G, A, B):
    x, y = G
    u = x * pow(B, -1, p) + A * pow(3 * B, -1, p)
    v = y * pow(B, -1, p)
    a = (3 - A ** 2) * pow(3 * B ** 2, -1, p)
    b = (2 * A ** 3 - 9 * A) * pow(27 * B ** 3, -1, p)
    G = (u, v)
    return (G, a, b)


G, a1, d1 = TEdwards_to_TEdwards(G, a, c, d)
G, A1, B1 = TEdwards_to_Montgomery(G, a1, d1)
G, a1, b1 = Montgomery_to_Weierstrass(G, A1, B1)

S = (98194560294138607903211673286210561363390596541458961277934545796708736630623,
     58504021112693314176230785309962217759011765879155504422231569879170659690008)

S, a2, d2 = TEdwards_to_TEdwards(S, a, c, d)
S, A2, B2 = TEdwards_to_Montgomery(S, a2, d2)
S, a2, b2 = Montgomery_to_Weierstrass(S, A2, B2)

assert G[1] ** 2 % p == G[0] ** 3 + a1 * G[0] + b1 % p
assert S[1] ** 2 % p == (S[0] ** 3 + a2 * S[0] + b2) % p
print("符合曲线方程 y^2 = x^3+ax+b")

E = EllipticCurve(GF(p), [a1, b1])

S = E(S)
G = E(G)

# 曲线基点的阶
n = G.order()
# 将椭圆曲线的阶分解为若干小素数
factors, exponents = zip(*factor(n))
print(factors)
primes = [factors[i] ^ exponents[i] for i in range(len(factors))][:-2]

dlogs = []
for fac in primes:
    t = int(int(n) // int(fac))
    dlog = discrete_log(t * S, t * G, operation="+")
    dlogs += [dlog]
    print("factor: " + str(fac) + ", Discrete Log: " + str(dlog))  # calculates discrete logarithm for each prime order

flag = crt(dlogs,primes)
print(long_to_bytes(flag))

输出结果：

符合曲线方程 ax^2+y^2 = c^2(1+dx^2y^2)
符合曲线方程 y^2 = x^3+ax+b
(2, 113, 28921, 749971, 87374629, 20375070355228635421, 1449497564615145263845097285948144399)
factor: 16, Discrete Log: 10
factor: 113, Discrete Log: 13
factor: 28921, Discrete Log: 1443
factor: 749971, Discrete Log: 406288
factor: 87374629, Discrete Log: 26189651
b'Tw1stzzzz'

对于上面代码中这一步 primes = [factors[i] ^ exponents[i] for i in range(len(factors))][:-2]

我们删掉了最后两位素数，是因为最后两个素数分别为：20375070355228635421和1449497564615145263845097285948144399 二者的乘积计算结果要大于flag的值，所以对 flag 进行的模运算结果仍然为 flag值，所以可以省略掉。

3.Hessian Curves DLP

题目代码

题目来源： CryptoCTF 2023 Barak

#!/usr/bin/env sage
from Crypto.Util.number import *
from flag import flag

def on_barak(P, E):
	c, d, p = E
	x, y = P
	return (x**3 + y**3 + c - d*x*y) % p == 0

def add_barak(P, Q, E):
	if P == (0, 0):
		return Q
	if Q == (0, 0):
		return P
	assert on_barak(P, E) and on_barak(Q, E)
	x1, y1 = P
	x2, y2 = Q
	if P == Q:
		x3 = y1 * (c - x1**3) * inverse(x1**3 - y1**3, p) % p
		y3 = x1 * (y1**3 - c) * inverse(x1**3 - y1**3, p) % p
	else:

		x3 = (y1**2*x2 - y2**2*x1) * inverse(x2*y2 - x1*y1, p) % p
		y3 = (x1**2*y2 - x2**2*y1) * inverse(x2*y2 - x1*y1, p) % p
	return (x3, y3)

def mul_barak(m, P, E):
	if P == (0, 0):
		return P
	R = (0, 0)
	while m != 0:
		if m & 1:
			R = add_barak(R, P, E)
		m = m >> 1
		if m != 0:
			P = add_barak(P, P, E)
	return R

def rand_barak(E):
	c, d, p = E
	while True:
		y = randint(1, p - 1)
		K = Zmod(p)
		P.<x> = PolynomialRing(K) 
		f = x**3 - d*x*y + c + y^3
		R = f.roots()
		try:
			r = R[0][0]
			return (r, y)
		except:
			continue

p = 73997272456239171124655017039956026551127725934222347
d = 68212800478915688445169020404812347140341674954375635
c = 1
E = (c, d, p)

P = rand_barak(E)

FLAG = flag.lstrip(b'CCTF{').rstrip(b'}')
m = bytes_to_long(FLAG) 
assert m < p
Q = mul_barak(m, P, E)
print(f'P = {P}')
print(f'Q = {Q}')
#P = (71451574057642615329217496196104648829170714086074852, 69505051165402823276701818777271117086632959198597714)
#Q = (40867727924496334272422180051448163594354522440089644, 56052452825146620306694006054673427761687498088402245)

题目分析

根据下列代码可知考察的是 Hessian 曲线：$H_{c,d}：x^3+y^3+c=dxy$

def on_barak(P, E):
	c, d, p = E
	x, y = P
	return (x**3 + y**3 + c - d*x*y) % p == 0

由 Q = mul_barak(m, P, E) 可知加密形式为：$Q=mP$ ，需要求解离散对数问题获取 m 的值。

已知量有曲线参数 $c,d,p$ ，以及曲线上的点 $P,Q$ ；

这里因为 $p-1$ 是一个光滑数，所以可能存在攻击离散对数的漏洞，如果曲线的阶也是一个光滑数，那么我们就可以利用 Pohlig-Hellman 算法来快速计算出离散对数的解。

1、我们先将Hessian 曲线映射为等价的 Weierstrass 曲线：

$$H：x^3+y^3+1=3Dxy⟼E_{a,b}：v^2=u^3-27D(D^3+8)\ast u+54(D^6-20D^3-8)$$

（其中 $3D=d$ ，令 $D=d/3$）

坐标的映射为：

$$(x,y)⟼(u,v)=(-9D^2+Bx,3B(y-1))$$

其中 $B=\frac{12(D^3-1)}{Dx+y+1}$

参考wiki介绍：Hessian form of an elliptic curve

2、然后用 sage上的函数 E.order() 计算曲线的阶，发现阶光滑且因式分解后都是小素数，所以我们可以利用椭圆曲线上的 Pohlig-Hellman 算法来进行求解计算离散对数。

# 曲线基点的阶
n = E.order()
print("曲线的阶：",n)
# 将椭圆曲线的阶分解为若干小素数
factors, exponents = zip(*factor(n))
print("阶进行因式分解：",factors)
print("素数指数：",exponents)
primes = [factors[i] ^ exponents[i] for i in range(len(factors))]

dlogs = []
for fac in primes:
    t = int(int(n) // int(fac))
    dlog = discrete_log(t * Q, t * P, operation="+")
    dlogs += [dlog]
    print("factor: " + str(fac) + ", Discrete Log: " + str(dlog))  # calculates discrete logarithm for each prime order

flag = crt(dlogs, primes)
assert flag*P == Q
print(long_to_bytes(flag))

输出结果：

曲线的阶： 73997272456239171124655016995459084401465136460086688
阶进行因式分解： (2, 3, 17, 2341, 23497, 500369, 5867327, 33510311, 13824276503, 67342255597)
素数指数： (5, 3, 1, 1, 1, 1, 1, 1, 1, 1)
factor: 32, Discrete Log: 1
factor: 27, Discrete Log: 1
factor: 17, Discrete Log: 13
factor: 2341, Discrete Log: 1489
factor: 23497, Discrete Log: 18879
factor: 500369, Discrete Log: 339241
factor: 5867327, Discrete Log: 4831116
factor: 33510311, Discrete Log: 31953114
factor: 13824276503, Discrete Log: 10804270204
factor: 67342255597, Discrete Log: 57315765050
b'\xa1UD8\x11T4\x92\xfa\x07Z-\x8e\xc1\xec;\xcc,\xab\xa0\x8b\x01'

我们已知： m < p，这里输出的 flag 是乱码，说明并不是真正的 m 值，但这里的flag值又满足 flag*P == Q 。

注： 当使用 Pohlig-Hellman 算法计算 $Q=mP$ 时，如果椭圆曲线是 Hessian 曲线并且其阶数比对应的 Weierstrass 曲线的阶数大，可能会导致计算出的 $m$ 值比真实值要大，存在偏差。

与普通的 Weierstrass 曲线相比，Hessian 曲线的计算公式可能导致在计算过程中产生额外的阶数倍数。这是因为 Hessian 曲线的点加法规则可能导致一些点的倍数比在对应的 Weierstrass 曲线上计算的倍数要大。

所以我们可以利用 $Q=(m\pm P.order())P$ 来找到准确的 $m$ 值

当我们将 $m$ 增加了 $P.order()$时，实际上相当于在群中循环地向前移动了一个完整的周期。在椭圆曲线密码学中，通常 $P.order()$ 是一个大质数，代表了点 $P$ 生成的子群的大小。如果 $Q=mP$，那么 $Q=(m\pm P.order())P$。因为在椭圆曲线上，$P.order()$ 个 $P$ 相加后会回到无穷远点（标记为 $O$），所以 $Q$ 的值与 $mP$ 的值相同。因此，$Q=mP$ 的性质不会因为增加 $P.order()$ 而改变

for i in range(flag,1, -P.order()):
    flag = long_to_bytes(i)
    print(flag)

输出结果：

b'\xa1UD8\x11T4\x92\xfa\x07Z-\x8e\xc1\xec;\xcc,\xab\xa0\x8b\x01'
b'\x99\x17\xa4W\xb9\x8d\xd4\x8ef\xb34\xf6&\x95\x94\xbd\x1e\x8d~\x94\xdd\xc5'
b'\x90\xda\x04wa\xc7t\x89\xd3_\x0f\xbe\xbei=>p\xeeQ\x890\x89'
b'\x88\x9cd\x97\n\x01\x14\x85@\n\xea\x87V<\xe5\xbf\xc3O$}\x83M'
b'\x80^\xc4\xb6\xb2:\xb4\x80\xac\xb6\xc5O\xee\x10\x8eA\x15\xaf\xf7q\xd6\x11'
b'x!$\xd6ZtT|\x19b\xa0\x18\x85\xe46\xc2h\x10\xcaf(\xd5'
b'o\xe3\x84\xf6\x02\xad\xf4w\x86\x0ez\xe1\x1d\xb7\xdfC\xbaq\x9dZ{\x99'
b'g\xa5\xe5\x15\xaa\xe7\x94r\xf2\xbaU\xa9\xb5\x8b\x87\xc5\x0c\xd2pN\xce]'
b'_hE5S!4n_f0rM_0F_3CC!!'
b'W*\xa5T\xfbZ\xd4i\xcc\x12\x0b:\xe52\xd8\xc7\xb1\x94\x167s\xe5'
b'N\xed\x05t\xa3\x94te8\xbd\xe6\x03}\x06\x81I\x03\xf4\xe9+\xc6\xa9'
b'F\xafe\x94K\xce\x14`\xa5i\xc0\xcc\x14\xda)\xcaVU\xbc \x19m'
b'>q\xc5\xb3\xf4\x07\xb4\\\x12\x15\x9b\x94\xac\xad\xd2K\xa8\xb6\x8f\x14l1'
b'64%\xd3\x9cATW~\xc1v]D\x81z\xcc\xfb\x17b\x08\xbe\xf5'
b'-\xf6\x85\xf3Dz\xf4R\xebmQ%\xdcU#NMx4\xfd\x11\xb9'

实现代码

sage 代码：

from Crypto.Util.number import *

p = 73997272456239171124655017039956026551127725934222347
d = 68212800478915688445169020404812347140341674954375635
c = 1
P = (71451574057642615329217496196104648829170714086074852, 69505051165402823276701818777271117086632959198597714)
Q = (40867727924496334272422180051448163594354522440089644, 56052452825146620306694006054673427761687498088402245)
D = d * pow(3, -1, p)


def on_Hessian(P):
    x, y = P
    if (x ** 3 + y ** 3 + c) % p == (d * x * y) % p:
        return True


def Hessian_to_Weierstrass(P):
    x, y = P
    A = (6 * (D ** 3 - 1) * (y + 9 * D ** 3 - 3 * D * x - 36)) * pow(
        (x + 9 * D ** 2) ** 3 + (3 * D ** 3 - D * x - 12) ** 3, -1, p)
    B = 12 * (D ** 3 - 1) * pow((D * x + y + 1), -1, p)
    u = -9 * D ** 2 + B * x
    v = 3 * B * (y - 1)
    return (u, v)


# 在海森曲线上
assert on_Hessian(P) == True
assert on_Hessian(Q) == True

# 映射为Weierstrass曲线上对应的坐标
P = Hessian_to_Weierstrass(P)
Q = Hessian_to_Weierstrass(Q)

a = -27 * D * (D ** 3 + 8)
b = 54 * (D ** 6 - 20 * D ** 3 - 8)

# 在魏尔斯特拉斯曲线上
E = EllipticCurve(GF(p), [a, b])
P = E(P)
Q = E(Q)

# 曲线基点的阶
n = E.order()
print("曲线的阶：", n)
# 将椭圆曲线的阶分解为若干小素数
factors, exponents = zip(*factor(n))
print("阶进行因式分解：", factors)
print("素数指数：", exponents)
primes = [factors[i] ^ exponents[i] for i in range(len(factors))]

dlogs = []
for fac in primes:
    t = int(int(n) // int(fac))
    dlog = discrete_log(t * Q, t * P, operation="+")
    dlogs += [dlog]
    print("factor: " + str(fac) + ", Discrete Log: " + str(dlog))  # calculates discrete logarithm for each prime order

flag = crt(dlogs, primes)
assert flag * P == Q
print(long_to_bytes(flag))

# 爆破准确的 flag 值
for i in range(flag, 1, -P.order()):
    flag = long_to_bytes(i)
    print(flag)

CCTF{_hE5S!4n_f0rM_0F_3CC!!}

方法二

使用sage自带的 EllipticCurve_from_cubic(F,P=None,morphism=True) 函数

• F– 具有有理系数的三个变量的齐次立方，作为多项式环元素，定义平滑的平面三次曲线$C$。
• P– 3元组(x,y,z)定义一个投影点$C$， 或者None。
• morphism– 布尔值（默认值True：）。如果True 则返回一个双有理同构$C$ 到 Weierstrass 椭圆曲线$E$，否则只返回$E$。

实现代码：

from Crypto.Util.number import *

p = 73997272456239171124655017039956026551127725934222347
d = 68212800478915688445169020404812347140341674954375635
c = 1
P = (71451574057642615329217496196104648829170714086074852, 69505051165402823276701818777271117086632959198597714)
Q = (40867727924496334272422180051448163594354522440089644, 56052452825146620306694006054673427761687498088402245)

# 曲线方程x^3+y^3+1=dxy
R.<x,y,z> = Zmod(p)[]   # 模p下的多项式环
cubic = x^3+y^3+c*z^3-d*x*y*z   # 构造有理系数的三元齐次立方
E = EllipticCurve_from_cubic(cubic,morphism=False) # 只返回一个Weierstrass形式的方程
EF = EllipticCurve_from_cubic(cubic,morphism=True) # 返回一个双有理同构$C$ 到 Weierstrass 椭圆曲线E

[In]: E
[Out]: Elliptic Curve defined by y^2 + 41848246294228984089267712788268151466593112763848031*x*y = x^3 + 8186198475255690096811875068967290988967924587416530*x^2 + 67041837351181533250843830910913901867849396421579685*x + 40673641660869730385035838198800394713965979400200155 over Ring of integers modulo 73997272456239171124655017039956026551127725934222347


[In]: EF
[Out]:
Scheme morphism:
  From: Projective Plane Curve over Ring of integers modulo 73997272456239171124655017039956026551127725934222347 defined by x^3 + y^3 + 5784471977323482679485996635143679410786050979846712*x*y*z + z^3
  To:   Elliptic Curve defined by y^2 + 41848246294228984089267712788268151466593112763848031*x*y = x^3 + 8186198475255690096811875068967290988967924587416530*x^2 + 67041837351181533250843830910913901867849396421579685*x + 40673641660869730385035838198800394713965979400200155 over Ring of integers modulo 73997272456239171124655017039956026551127725934222347
  Defn: Defined on coordinates by sending (x : y : z) to
        (9167218108590654299487009181696317855454416112239889*y + 3279765642505352830157211697766991850038465138705271*z : 5784471977323482679485996635143679410786050979846712*y : 64249990329941531194735759683164021211133831031453155*x + 22688267644108754864364901498104863155005679028019334*y + 42580013058247064141816874924090310799821854717107287*z)

P = EF(P) # 映射坐标
Q = EF(Q)
m = discrete_log(Q,P,operation="+")
print(m)

m = 1780694557271320552511299360138314441283923223949197
assert  Q == m*P
print(long_to_bytes(m))
# 爆破准确的m值
for i in range(m,p-1, P.order()):
    flag = long_to_bytes(i)
    print(flag)

4.Huff Curves DLP

题目代码

题目来源：[DASCTF 2023] CB curve

from Crypto.Util.number import *
from random import randint
from secret import flag,order


class CB_curve:
    def __init__(self):
        self.p = 1141741939958844590498346884870015122543626602665954681008204697160652371664923
        self.a = 727131475903635498678013730344448225340496007388151739960305539398192321065043
        self.b = 840714623434321649308065401328602364673881568379142278640950034404861312007307

    def add(self, P, Q):
        if P == -1:
            return Q
        (x1, y1) = P
        (x2, y2) = Q
        x3 =  (x1+x2)*(1+self.a*y1*y2)*inverse((1+self.b*x1*x2)*(1-self.a*y1*y2),self.p)% self.p
        y3 =  (y1+y2)*(1+self.b*x1*x2)*inverse((1-self.b*x1*x2)*(1+self.a*y1*y2),self.p)% self.p
        return (x3, y3)

    def mul(self, x, P):
        Q = -1
        while x > 0:
            if x & 1:
                Q = self.add(Q, P)
            P = self.add(P, P)
            x = x >> 1
        return Q
    
    def negG(self,G):
        return self.mul(order-1,G)

ecc = CB_curve()
G = (586066762126624229327260483658353973556531595840920560414263113786807168248797, 66727759687879628160487324122999265926655929132333860726404158613654375336028)
P = (ecc.mul(bytes_to_long(flag),G)[0],randint(1,ecc.p))
Q = (460843895959181097343292934009653542386784127282375019764638432240505304648101, 739422832583403823403837831802136107593509589942947902014204968923412689379907)

e = randint(1,p)
pl = [ecc.add(P,ecc.mul(10-i,ecc.negG(Q)))[0] + e for i in range(10)]
ph = [ecc.add(P,ecc.mul(10-i,Q))[0] + e for i in range(10)]

print(pl)
print(ph)

题目分析

分析代码可知曲线加法满足 Huff 曲线方程： $x(a{y}^2-1)=y(b{x}^2-1)$

点加公式为：

$$(x_1,y_1)+(x_2,y_2)=(\frac{(x_1+x_2)(1+a{y}_1{y}_2)}{(1+b{x}_1{x}_2)(1-a{y}_1{y}_2)},\frac{(y_1+y_2)(1+b{x}_1{x}_2)}{(1-b{x}_1{x}_2)(1+a{y}_1{y}_2)})$$

将 Huff 曲线映射为 Weierstrass 曲线公式为：

$x(a{y}^2-1)=y(b{x}^2-1)⟼v^2=u^3+a_2{u}^2+a_{4}u$

其中：

$$u=\frac{bx-ay}{y-x},v=\frac{b-a}{y-x}$$$$a_2=a+b,a_4=ab$$

由函数 negG() 可知，$-G=(G.order-1)\ast G$

def negG(self, G):
    return self.mul(order - 1, G)

已知量有Huff 曲线的参数信息 $p,a,b$ 和 变量 $pl，ph$ 的值；未知量有 $e,order,P$

其中 $P.x=[flag\ast G].x$ ，所以我们需要先求出点 $P$ 的坐标值。

P = (ecc.mul(bytes_to_long(flag), G)[0], randint(1, ecc.p))
pl = [ecc.add(P, ecc.mul(10 - i, ecc.negG(Q)))[0] + e for i in range(10)]
ph = [ecc.add(P, ecc.mul(10 - i, Q))[0] + e for i in range(10)]

根据 $pl,ph$ 的公式可知，这里进行的点加的结果只取 $x$ 坐标，所以我们下面公式只需要写关于 $x$ 坐标的：

令 $x_Q=[(10-i)\ast Q].x$ ，$y_Q=[(10-i)\ast Q].y$

$pl=\frac{(x_P-x_Q)\ast (1-a{y}_P{y}_Q)}{(1-b{x}_P{x}_Q)\ast (1+a{y}_P{y}_Q)}+e$ ，$ph=\frac{(x_P+x_Q)\ast (1+a{y}_P{y}_Q)}{(1+b{x}_P{x}_Q)\ast (1-a{y}_P{y}_Q)}+e$ ，$pl$ 这里的是 $-Q$ 所以公式中是负数。

有： $(pl-e)\ast (ph-e)=\frac{x_p^2-x_Q^2}{1-b^2{x}_P^2{x}_Q^2}$

因为上述等式提供的数据有多组，等式中的未知量为 $e,x_p$，所以可以二元构造多项式为：

$(pl-e)\ast (ph-e)\ast (1-b^2{x}_P^2{x}_Q^2)=x_p^2-x_Q^2$

我们可以使用 格劳布纳基础 求解多项式方程组 ，在SageMath中可以使用groebner_basis()函数来计算一个理想的格劳布纳基础。

class CB_curve:
    def __init__(self):
        self.p = 1141741939958844590498346884870015122543626602665954681008204697160652371664923
        self.a = 727131475903635498678013730344448225340496007388151739960305539398192321065043
        self.b = 840714623434321649308065401328602364673881568379142278640950034404861312007307

    def add(self, P, Q):
        if P == -1:
            return Q
        (x1, y1) = P
        (x2, y2) = Q
        x3 = (x1 + x2) * (1 + self.a * y1 * y2) * inverse((1 + self.b * x1 * x2) * (1 - self.a * y1 * y2),
                                                          self.p) % self.p
        y3 = (y1 + y2) * (1 + self.b * x1 * x2) * inverse((1 - self.b * x1 * x2) * (1 + self.a * y1 * y2),
                                                          self.p) % self.p
        return (x3, y3)

    def mul(self, x, P):
        Q = -1
        while x > 0:
            if x & 1:
                Q = self.add(Q, P)
            P = self.add(P, P)
            x = x >> 1
        return Q

    def negG(self, G):
        return self.mul(order - 1, G)


ecc = CB_curve()

PR.<e,xp> = Zmod(p)[]  # 定义模p下的多项式环
f = []
for i in range(10):
    xp = xp
    xq = ecc.mul(10-i,Q)[0]
    f.append(xp^2-xq^2 - (pl[i]-e)*(ph[i] - e)*(1-b^2*xp^2*xq^2))
I= Ideal(f).groebner_basis()
print("多项式方程组的解：",I)

PR.<xp> = Zmod(p)[]
f2 = xp^2 + 219493165434454878473973957507132663767650700404392831423708684433961924200902
xp = f2.roots()
print("P点的x坐标",xp)

多项式方程组的解： [xp^2 + 219493165434454878473973957507132663767650700404392831423708684433961924200902, e + 716700711017198421972376297958894204723153539777056104579499803899129208364755]

P点的x坐标 [(902098404790942932698136414095200590033317733070113807365359132832241907267881, 1), (239643535167901657800210470774814532510308869595840873642845564328410464397042, 1)]

求出 P 点的坐标后，由 P = (ecc.mul(bytes_to_long(flag), G)[0], randint(1, ecc.p)) 可知，ecc.mul(bytes_to_long(flag), G)[0] 与 P 有相同的 x 坐标，那么我们可以在曲线上找到他们共有的 y坐标值。所以可以直接看成是 $P=flag\ast G$ ，$G$为生成点。

计算 G.order() 发现是光滑阶且因式分解后是小素数，那么我们可以使用 Pohlig-Hellman 算法进行离散对数后的求解。

最后校验：assert flag * G == P 即可

实现代码

sagemath实现：

from Crypto.Util.number import *

p = 1141741939958844590498346884870015122543626602665954681008204697160652371664923
a = 727131475903635498678013730344448225340496007388151739960305539398192321065043
b = 840714623434321649308065401328602364673881568379142278640950034404861312007307

G = (586066762126624229327260483658353973556531595840920560414263113786807168248797,
     66727759687879628160487324122999265926655929132333860726404158613654375336028)
Q = (460843895959181097343292934009653542386784127282375019764638432240505304648101,
     739422832583403823403837831802136107593509589942947902014204968923412689379907)

pl = [908996880816674413953945844149350915331956247471480600840221415119794882139724,
 971918808384910355828135603762747020183688585728289421786279444571287619529246,
 1285550352531583269956802123237391199017403081800977678246201935580429758051904,
 1551774945769448705387900437472951015954157193946719575845523359198154668857591,
 676185408751480221545400062950292727848016906516506232986883519673765317932582,
 1250300209784131850574858927023046353058343552115735540789593580037130054384362,
 1298409778422699298367007023890818793557023853717180295526932023194697263501748,
 1332552452292482549702793642987623159617988974910321945878093492007278710993114,
 1030239404875082841481045525469865919289388171602293245905162820968158543176773,
 1154148024180033719999293176590867264297899817449945744942661351655533433871621]
ph = [584297112520340495757457954416165393828472756298945167299482077258411155766756,
 886432149227960827335266910774569034430464592640209168563805700117347063152246,
 613528590036968449893421430816319461615130635882647544978722093413694101540550,
 576162106332135829961234799085370038425761945928004579456101802617485243023987,
 627570890346195626159365118862437334953500165050236216404858019114288681512171,
 1015503424232985454098149884321288932492551183126601131968495641510550575005042,
 1532737675157046782602115678180407262847166210963507805526455422934164759886583,
 1540047002602145805476906585925538790245968214992837106009502002588479779602195,
 505097517314409449404205152068185149808364887623922221197462411159844816865696,
 873498218680784138428154510303205366133389839886911286745954821800632158315951]

# 满足 Huff曲线
assert G[0] * (a * G[1] ** 2 - 1) % p == G[1] * (b * G[0] ** 2 - 1) % p
assert Q[0] * (a * Q[1] ** 2 - 1) % p == Q[1] * (b * Q[0] ** 2 - 1) % p
print("1. 满足 Huff曲线")


class CB_curve:
    def __init__(self):
        self.p = 1141741939958844590498346884870015122543626602665954681008204697160652371664923
        self.a = 727131475903635498678013730344448225340496007388151739960305539398192321065043
        self.b = 840714623434321649308065401328602364673881568379142278640950034404861312007307

    def add(self, P, Q):
        if P == -1:
            return Q
        (x1, y1) = P
        (x2, y2) = Q
        x3 = (x1 + x2) * (1 + self.a * y1 * y2) * inverse((1 + self.b * x1 * x2) * (1 - self.a * y1 * y2),
                                                          self.p) % self.p
        y3 = (y1 + y2) * (1 + self.b * x1 * x2) * inverse((1 - self.b * x1 * x2) * (1 + self.a * y1 * y2),
                                                          self.p) % self.p
        return (x3, y3)

    def mul(self, x, P):
        Q = -1
        while x > 0:
            if x & 1:
                Q = self.add(Q, P)
            P = self.add(P, P)
            x = x >> 1
        return Q

    def negG(self, G):
        return self.mul(order - 1, G)



def Huff_to_E(P):
    x, y = P
    u = (b * x - a * y) * pow(y - x, -1, p)
    v = (b - a) * pow(y - x, -1, p)
    return (u, v)


ecc = CB_curve()

PR.<e,xp> = Zmod(p)[]  # 定义模p下的多项式环
f = []
for i in range(10):
    xp = xp
    xq = ecc.mul(10-i,Q)[0]
    f.append(xp^2-xq^2 - (pl[i]-e)*(ph[i] - e)*(1-b^2*xp^2*xq^2))
I= Ideal(f).groebner_basis()
print("多项式方程组的解：",I)

PR.<xp> = Zmod(p)[]
f2 = xp^2 + 219493165434454878473973957507132663767650700404392831423708684433961924200902
xp = f2.roots()
print("P点的x坐标",xp)

xp = 902098404790942932698136414095200590033317733070113807365359132832241907267881

# 求解 P 点对应的 y 坐标
PR.<yp> = Zmod(p)[]
f3 = xp*(a*yp^2-1) - yp*(b*xp^2-1)
yp = f3.roots()
print("P点的y坐标",yp)
yp = 672929595307990944197873882889709005621738844588134711458648048321447534353147


# 映射为 Weierstrass 上的坐标
G = Huff_to_E(G)
Q = Huff_to_E(Q)
P = (xp,yp)
P = Huff_to_E(P)

a2 = a+b
a4 = a*b
E = EllipticCurve(GF(p),[0,a2,0,a4,0])
G = E(G)
Q = E(Q)
print("2.成功映射为 Weierstrass曲线")


P = E(P)
print("P点在E上，且值为：",P)



print("\n")
# 曲线基点的阶
n = G.order()
print("曲线的阶：", n)
# 将椭圆曲线的阶分解为若干小素数
factors, exponents = zip(*factor(n))
print("阶进行因式分解：", factors)
print("素数指数：", exponents)
primes = [factors[i] ^ exponents[i] for i in range(len(factors))][:-1]

dlogs = []
for fac in primes:
    t = int(int(n) // int(fac))
    dlog = discrete_log(t * P, t * G, operation="+")
    dlogs += [dlog]
    print("factor: " + str(fac) + ", Discrete Log: " + str(dlog))  # calculates discrete logarithm for each prime order

flag = crt(dlogs, primes)
assert flag * G == P
print(long_to_bytes(flag))

b'DASCTF{goodathuff}'

5. 奇异椭圆曲线 DLP

奇异椭圆曲线知识点学习：奇异椭圆曲线

题目代码

题目来源：[XYCTF 2024] easy_ecc

from Crypto.Util.number import *
from hashlib import sha256
from secret import flag, secret,SECRET

assert flag[6:-1] == sha256(long_to_bytes(secret)).hexdigest().encode()


class ECC_easy:
    def __init__(self):
        self.a = 1365855822212045061018261334821659180641576788523935479
        self.b = 17329427219955161804703200105884322768704934833367341
        self.p = 1365855822212045061018261334821659180641576788523935481

    def add(self, P, Q):
        mul_inv = lambda x: pow(x, -1, self.p)
        x1, y1 = P
        x2, y2 = Q
        if P!=Q:
            l=(y2-y1)*inverse(x2-x1,self.p)%self.p
        else:l=(3*x1**2+2*self.a*x1+1)*inverse(2*self.b*y1,self.p)%self.p
        temp1 = (self.b*l**2-self.a-x1-x2)%self.p
        temp2 = ((2*x1+x2+self.a)*l-self.b*l**3-y1)%self.p
        x3 = temp1
        y3 = temp2
        return x3, y3

    def mul(self, x, P):
        Q = SECRET
        x = x % self.p
        while x > 0:
            if x & 1:
                Q = self.add(Q, P)
            P = self.add(P, P)
            x >>= 1
        return Q

    def ispoint(self, x, y):
        return (self.a * x ** 2 + x ** 3+x) % self.p == (self.b * y ** 2) % self.p


ecc = ECC_easy()
LLLL = (1060114032187482137663886206406014543797784561116139791,752764811411303365258802649951280929945966659818544966)
assert ecc.ispoint(LLLL[0], LLLL[1])
END = ecc.mul(secret, LLLL)
print(END)

# (695174082657148306737473938393010922439779304870471540,414626357054958506867453055549756701310099524292082869)

题目分析

我们需要求的未知量为：secret

已知量为：曲线参数 a、b、p，LLLL这个变量为一个基点值，变量END 值为标量乘法计算后的结果，有END = mul(secret, LLLL)

由题目给的如下函数可知，曲线方程为蒙哥马利曲线方程

def ispoint(self, x, y):
    return (self.a * x ** 2 + x ** 3 + x) % self.p == (self.b * y ** 2) % self.p

我们先将蒙哥马利曲线方程转化为：魏尔斯特拉斯曲线方程，这样方便在 sagemath中操作

def Montgomery_to_Weierstrass(G, A, B):
    x, y = G
    u = (x * pow(B, -1, p) + A * pow(3 * B, -1, p)) % p
    v = y * pow(B, -1, p) % p
    a = (3 - A ** 2) * pow(3 * B ** 2, -1, p) % p
    b = (2 * A ** 3 - 9 * A) * pow(27 * B ** 3, -1, p) % p
    G = (u, v)
    return (G, a, b)

对转化后的值，使用 sagemath 的E = EllipticCurve(GF(p),[a1,b1]) 函数构造椭圆曲线，这时候提示报错信息：

ArithmeticError: invariants (0, 0, 0, 1098066776930223927329092382214459309226361965213, 1263248714105743124314734095577181018742053879965591734) define a singular curve

这是因为 sage 中的 EllipticCurve() 函数不能使用奇异椭圆曲线，我们需要通过公式将其转化为非奇异曲线。

参考密码学书籍：Elliptic Curves: Number Theory and Cryptography 第二章 2.10 Singular Curves ,P59-P62

给出了正式的描述（和证明）。简而言之，该算法是：

• 找到奇点
• 平移曲线，使点位于 (0, 0)
• 将曲线上的点映射到有限域中的元素（映射取决于曲线是尖点还是节点）
• 计算有限域中的离散对数

1、定理 2.30

设曲线 $E：y^2=x^3$ ，设 $E_{ns}(K)$ 是这条曲线上坐标在 $K$ 中的非奇异点， including point at infinity

$\infty =(0:1:0)$ ，映射关系如下：$(x,y)\mapsto \frac{x}{y}$

• 通过使用这个定理，如果 $G$ 是奇异曲线上的一点，且 $P=dG$ ，那么求 $d$ 有如下方法：

$g=G.x\times G.y^{-1}(modp)$， $y=P.x\times P.y^{-1}(modp)$

$d=g^{-1}\times y(modp)$

2、定理 2.31

设曲线 $E：y^2=x^2(x+a)$ 且 $0\ne a\in K$ 。设 $E_{ns}(K)$ 是这条曲线上坐标在 $K$ 中的非奇异点。

设 ${\alpha }^2=a$ ，映射关系如下： $(x,y)\mapsto \frac{y+\alpha x}{y-\alpha x}$

根据上面的定理 2.31 和相关知识点可得：

PR.<x,y> = GF(p)[]
f = x ^ 3 + a1 * x + b1
# 多项式方程构造曲线
C = Curve(-y ^ 2 + f)
print(C)
# 计算椭圆曲线上的奇点
singular_point = C.singular_points()
print("\n奇点为：", singular_point)

x0 = singular_point[0][0]

# 将奇点平移到(0,0)
f_ = f.subs(x=x + x0)
# 其余所有在曲线上的点也需要平移
P_ = (P[0] - x0, P[1])
G_ = (G[0] - x0, G[1])
print("平移后的曲线方程为：", f_.factor())

1、找到奇点

2、将奇点： (821598549758978983064709974196127522156033448658296567, 0) 通过$(x-821598549758978983064709974196127522156033448658296567,y-0)$ 平移后到点 $(0,0)$。

（注意：这里后面解题过程中记得将曲线上的已知点也要平移）

3、将曲线上的点映射到有限域中的元素

可知： 平移后的曲线方程为： $(x+1098939827064891888175868587766723385826523557450954220)\ast x^2$

根据上面条件可知：$y^2=(x+1098939827064891888175868587766723385826523557450954220)\ast x^2$ ，

令 $a=1098939827064891888175868587766723385826523557450954220$ ，在有限域 GF(p) 上 ${\alpha }^2=a$ ，则：

$\alpha =361925395984543395963949683281348788604680612880798808$

根据映射关系： $(x,y)\mapsto \frac{y+\alpha x}{y-\alpha x}$ ，得 $P,G$ 在有限域上的值，然后求离散对数即可。

实现代码

from Crypto.Util.number import *
from hashlib import sha256

a = 1365855822212045061018261334821659180641576788523935479
b = 17329427219955161804703200105884322768704934833367341
p = 1365855822212045061018261334821659180641576788523935481
G = (1060114032187482137663886206406014543797784561116139791, 752764811411303365258802649951280929945966659818544966)
P = (695174082657148306737473938393010922439779304870471540, 414626357054958506867453055549756701310099524292082869)


def Montgomery_to_Weierstrass(G, A, B):
    x, y = G
    u = (x * pow(B, -1, p) + A * pow(3 * B, -1, p)) % p
    v = y * pow(B, -1, p) % p
    a = (3 - A ** 2) * pow(3 * B ** 2, -1, p) % p
    b = (2 * A ** 3 - 9 * A) * pow(27 * B ** 3, -1, p) % p
    G = (u, v)
    return (G, a, b)


G, a1, b1 = Montgomery_to_Weierstrass(G, a, b)
P, a2, b2 = Montgomery_to_Weierstrass(P, a, b)
print("G =", G)
print("P =", P)
print("a1 =", a1)
print("b1 =", b1)

G = (804603363012007759329983017410816754946539644939, 668360700828957783980888938878566241242911721895008218)
P = (933414165833077907509715600260551365988944141925739220, 121346737700219338084994830488363509910434835223666824)
a1 = 1098066776930223927329092382214459309226361965213
b1 = 1263248714105743124314734095577181018742053879965591734

# 判断曲线是否为奇异曲线
if (4 * a1 ^ 3 + 27 * b1 ^ 2) % p == 0:
    print("[*] 奇异椭圆曲线!!!")
else:
    print("[*] 非奇异椭圆曲线!!!")

PR.<x,y> = GF(p)[]
f = x ^ 3 + a1 * x + b1
# 多项式方程构造曲线
C = Curve(-y ^ 2 + f)
print(C)
# 计算椭圆曲线上的奇点
singular_point = C.singular_points()
print("\n奇点为：", singular_point)

x0 = singular_point[0][0]

# 将奇点平移到(0,0)
f_ = f.subs(x=x + x0)
# 其余所有在曲线上的点也需要平移
P_ = (P[0] - x0, P[1])
G_ = (G[0] - x0, G[1])
print("平移后的曲线方程为：", f_.factor())

a = 1098939827064891888175868587766723385826523557450954220
alpha = GF(p)(a).square_root()
print("alpha =", alpha)

# 将曲线上的点映射为有限域上的点
P_ = (P_[1] + alpha * P_[0]) * (P_[1] - alpha * P_[0]) ^ (-1) % p
G_ = (G_[1] + alpha * G_[0]) * (G_[1] - alpha * G_[0]) ^ (-1) % p

secret = discrete_log(P_, G_)
print("secret =", secret)
flag = b"XYCTF{"+ sha256(long_to_bytes(secret)).hexdigest().encode() +b"}"
print(flag)

b'XYCTF{ec9a6e17537e81b7f593f65f7e2ca5d575e6b34c504c24e4afb40c1e9dc4be0d}'

6.Smart’s Attack 算法

适用于解决椭圆曲线上的离散对数问题，当 $E.order()=p$ ，（这里的 $p$ 为有限域） 可以使用该算法快速求解离散对数 $P=k\ast G$ 的问题。

相关论文：https://wstein.org/edu/2010/414/projects/novotney.pdf

例题： Exceptional Curves

实现EXP:

def SmartAttack(P, G, p):
    """
    用于求解当 E(Fp) = p 时,P=kG ,私钥 k 的值
    :param G: 基点 
    :param P: 公钥
    :param p: 模数
    :return: 私钥 k
    """
    E = G.curve()
    Eqp = EllipticCurve(Qp(p, 2), [ ZZ(t) + randint(0,p)*p for t in E.a_invariants() ])

    P_Qps = Eqp.lift_x(ZZ(G.xy()[0]), all=True)
    for P_Qp in P_Qps:
        if GF(p)(P_Qp.xy()[1]) == G.xy()[1]:
            break

    Q_Qps = Eqp.lift_x(ZZ(P.xy()[0]), all=True)
    for Q_Qp in Q_Qps:
        if GF(p)(Q_Qp.xy()[1]) == P.xy()[1]:
            break

    p_times_P = p*P_Qp
    p_times_Q = p*Q_Qp

    x_P,y_P = p_times_P.xy()
    x_Q,y_Q = p_times_Q.xy()

    phi_P = -(x_P/y_P)
    phi_Q = -(x_Q/y_Q)
    k = phi_Q/phi_P
    return ZZ(k)