from itertools import combinations

def cornacchia(d, m):
    assert(gcd(d,m)==1)
    try:
        rs = [ ZZ(x) for x in Zmod(m)(-d).sqrt(all=True) ]
    except:
        return None
    for r in rs:
        r0 = m; r1 = r
        while r0^2 >= m:
            r0, r1 = r1, r0 % r1
        if is_square((m - r0^2)/d):
            return r0, sqrt((m - r0^2)/d)
    return None


def find_curve_of_order(n):
    r = 0
    logn = int(n.log(2).n())
    while True:
        p = [ legendre_symbol(-1, p)*p for p in prime_range(max(3,r*logn), (r+1)*logn) if legendre_symbol(n, p) == 1 ]
        S = [ p for p in p ]
        for i in range(1, len(p)):
            for L in combinations(S, i):
                D = prod(L)
                if D % 8 != 5:# or D.log(2).n() > (r * n.log(2).n())^2:
                    continue
                solution = cornacchia(-D, 4*n)
                if solution is None:
                    continue
                (x, y) = solution
                if is_prime(n+1+x):
                    p = n+1+x
                elif is_prime(n+1-x):
                    p = n+1-x
                else:
                    continue
                P = hilbert_class_polynomial(D)
                roots = P.roots(ring=GF(p))
                if len(roots) > 0:
                    E = EllipticCurve_from_j(roots[0][0]).change_ring(GF(p))
                    if E.order() == n:
                        return E
                    else:
                        return E.quadratic_twist()
        r += 1



sage: time E = find_curve_of_order(2**127-1)
CPU times: user 583 ms, sys: 7 µs, total: 583 ms
Wall time: 583 ms
sage: E
Elliptic Curve defined by y^2 = x^3 + 29825378467770542984135446927558639468*x + 54083907264064377833823915102218955679 over Finite Field of size 170141183460469231709232497652167941199
sage: E.order()
170141183460469231731687303715884105727
sage: 2**127-1
170141183460469231731687303715884105727