my solver for DDLP in HackTM 2023 Finals

solve.sage

primes = list(prime_range(1000000))

# 4p = 1 + Dy^2, p = p1 * p2 * ... * pm + 1
# 4 * p1 * ... * pm + 3 = D * y^2 (p1 should be 2)
# 8 * p2 * ... * pm = D * y^2 - 3 = (sqrt(D) * y - sqrt(3)) * (sqrt(D) * y + sqrt(3))
# when D = 27, r.h.s. = 3 * (3*y - 1) * (3*y + 1)
cnt = 0
y3_min = int(sqrt(2 ** 255 * 4 // 3))
y3_max = int(sqrt(2 ** 256 * 4 // 3))
while True:
    cnt += 1
    if cnt % 10000 == 0:
        print(cnt)
    # tmp is a candidate for 3*y - 1, which must be smooth
    tmp = 4
    while tmp < y3_min:
        tmp *= choice(primes)
    if not y3_min <= tmp <= y3_max:
        continue
    if (tmp + 1) % 3 != 0:
        continue
    y = (tmp + 1) // 3
    factors = list(map(lambda x: x[0], factor(3*y + 1)))
    max_factor = factors[-1]
    if max_factor < 2**40:  # if 3*y + 1 is also smooth, p - 1 should be smooth
        tmp_p = 3 * (3 * y - 1) * (3 * y + 1) // 8 * 2 + 1
        if is_prime(tmp_p):
            break
p = tmp_p
D = 27
j = -12288000
a = -3 * j * pow(j - 1728, -1, p) % p
b = 2 * j * pow(j - 1728, -1, p) % p
E = EllipticCurve(GF(p), [a, b])
if E.order() != p:
    E = E.quadratic_twist()