Elliptic curve tower-cycle for Curve25519

readme.md

Elliptic curve tower-cycle for Curve25519

Some advanced cryptographic protocols require a "curve cycle", which is a pair of elliptic curves Ep and Eq over prime fields GF(p) and GF(q) such that #Ep = q and #Eq = p, i.e. the scalar field of each curve exactly equals to the base field of the other curve.

Unfortunately, Curve25519 has a composite order (divisible by 8), so it cannot be part of a curve cycle (both p and q need to be primes).

However, if we could find a curve cycle with p = 2^255 - 19, points on Curve25519 could enter the cycle via curve Eq, which has an order of p. We call this arrangement a "tower-cycle":

   -- Ep --
 /          \
|            ^
v            |
 \          /
   -- Eq --
      ^
      |
  Curve25519
    

The search method

To find a curve of exact order p = 2^255 - 19, we will use an algorithm described in the 2007 paper "Constructing elliptic curves of prime order" [1].

The algorithm roughly consists of the following steps:

1. Searching for CM field discriminant D values that produce a trace t such that q = p + 1 ± t is a prime number.
2. Finding a root j0 of the Hilbert class polynomial over GF(q).
3. Calculating a = 27*j0/(4*(1728 - j0))
4. Finally, we have Eq: y^2 = x^3 + a*x - a or its quadratic twist (one of them will have order p).

The Hilbert class polynomial has more than one root, so there will be multiple curves for each prime.

Due to the symmetry of the algorithm between p and q, we can use the same process to find Ep (except we can skip the first step because we already know the value of p).

Search conditions

For performance and security reasons, we place additional conditions on q, Eq and Ep. We split the search conditions into two groups.

Q-conditions

Q-conditions restrict the choice of the prime q.

1. We want q < 2^255 so we can encode compressed points of Eq in 256 bits. This means we will only use q = p + 1 - t in step 1 above.
2. We want q = 3 mod 4 for efficient square root calculation mod q [2].
3. We want 2^256 - 2*q < 2^128, which makes certain modular reduction algorithms more efficient in software (notably the Crandall's reduction algorithm [3]).
4. We want both Ep and Eq to have twist security of at least 100 bits. Although this is not required for the intended use cases of a tower cycle, it is a common search condition recommended by SafeCurves [4].

C-conditions

C-conditions place restrictions on the coefficients a and b of the curve equation y^2 = x^3 + a*x + b. They don't affect the choice of the prime q, but they allow us to select better candidates among all the isogeneous curves of prime order p, q over GF(q), GF(p).

1. We want a = -3, which allows for more efficient group law formulas [5].
2. We want b to be a non-square number, which prevents the point x = 0 from being on the curve. This protects from some power analysis attacks [6] and is also beneficial for avoiding the point at infinity when hashing from Eq to Ep and vice versa.
3. We want the point x = 1 to exist on the curve. This allows us to select the point [1, y] with y even, as the base point for both curves.

Results

The attached script tower_cycle_25519.sage implements the search. The discriminant with the lowest absolute value that produces a prime q that matches our conditions Q1 to Q4 is D = -31750123 and the prime is:

q = 0x7ffffffffffffffffffffffffffffffff735481d1969f317f9850b68df11df53

The resulting curve cycle consists of two curves Ep and Eq, which we will call "Helios" and "Selene". They correspond to the lowest values of j0 that match our conditions C1 to C3.

Helios: y^2 = x^3 - 3*x + 17523451383230374900436292617863907649717438939964238673872692863501483215968 over GF(p)
Selene: y^2 = x^3 - 3*x + 25675911719867737339625140396204798989996478324626569376465022644547366285284 over GF(q)

We define the base point of each curve as the point with x = 1 and an even y coordinate:

Helios: G = [1, 43927350165885181914572701368652294970994947138804342515295004363921039321018]
Selene: G = [1, 25798700515841442074436724357845010259191889815205036617843407906692357567936]

The attached file tower_cycle_25519_output.txt shows the full output of the search script, including all parameters that were not selected and the reason that caused them not to be selected.

Security evaluation

Both Helios and Selene meet some of the SafeCurves criteria [4].

1. Fields ✔️

Both p and q are primes.

2. Equation ✔️

For both Helios and Selene, b != ±2, so the curves are not singular.

3. Base points ✔️

For both curves, the base point G lies on the curve. For Helios, q*G == O and for Selene, p*G == O.

4. Rho ✔️

Both curves provide 127.3 bits of ECDLP security, which is more than Curve25519.

5. Transfers ✔️

Since p != q, both curves are safe from additive transfers. The embedding degree of Helios is (q-1)/2 and the embedding degree of Selene is (p-1)/12, so both are also safe from multiplicative transfers.

6. Discriminants ❌

Due to the way they were constructed, both curves have a CM field discriminant value of -31750123.

7. Rigidity ✔️

The curve generation process is completely explained.

8. Ladders ❌

Both curves have a cofactor of 1, so they don't support the Montgomery ladder.

9. Twists ✔️

Helios has 105.2 bits of twist security and Selene has 110.5 bits of twist security. While this meets the minimum requirements of SafeCurves, we don't recommend using either of these curves with single-coordinate ladder formulas.

10. Completeness ❌

The curves don't support complete formulas that meet the SafeCurves requirements.

11. Indistinguishability ✔️

Both curves support indistinguishability using Elligator Squared [7].

Errata

Due to a bug in the search script, this document previously incorrectly stated that the first discriminant that matches the search conditions Q1 to Q3 was D = -7857907. We thank Veridise for finding the mistake during their security audit.

The search condition Q4 was added later as recommended by the security audit. This caused a different pair of curves to be selected in the end.

References

[1] Constructing elliptic curves of prime order: https://arxiv.org/abs/0712.2022

[2] A simple algorithm for finding square root modulo p: https://arxiv.org/pdf/2008.11814.pdf

[3] Crandall's modular reduction: https://patents.google.com/patent/US5159632A/en

[4] SafeCurves: choosing safe curves for elliptic-curve cryptography https://safecurves.cr.yp.to/index.html

[5] Jacobian coordinates with a=-3 for short Weierstrass curves: https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html

[6] A Refined Power-Analysis Attack on Elliptic Curve Cryptosystems https://www.iacr.org/archive/pkc2003/25670199/25670199.pdf

[7] Elligator Squared: Uniform Points on Elliptic Curves of Prime Order as Uniform Random Strings https://eprint.iacr.org/2014/043.pdf

tower_cycle_25519.sage

#!/usr/bin/env sage
# -*- coding: utf-8 -*-

# This script finds efficient elliptic curve tower-cycles for Curve25519.
# The algorithm is described in this paper: https://arxiv.org/pdf/0712.2022.pdf

import sys
from multiprocessing import Pool, cpu_count
from traceback import print_exc

p = 2^255 - 19

# https://en.wikipedia.org/wiki/Cornacchia%27s_algorithm
# Finds [x,y] that solves x^2+d*y^2=m

def cornacchia_inner(r, m, d):
    sm = isqrt(m)
    t = m
    while r > sm:
        rn = t % r
        t = r
        r = rn
    s2 = m - r^2
    if s2 % d != 0:
        return (None, "s2 not divisible by d")
    s2 //= d
    s = isqrt(s2)
    if s*s == s2:
        return (r, s)
    return (None, "sqrt(s2) not an integer")

def cornacchia(d, m):
    rall = Mod(-d, m).sqrt(extend=False, all=True)
    if not rall:
        return (None, "No sqrt mod m")
    # test all square roots
    for r in rall:
        (x, y) = cornacchia_inner(int(r), m, d)
        if x:
            break
    return (x, y)

# Find an isomorphism with a specific value of the curve parameter "a".
def get_iso(E, ax):
    a = E.a4()
    b = E.a6()
    f = E.base_field()
    u4 = a / ax
    us = u4.nth_root(4,all=True)
    if not us:
        return None
    # sorting for determinism
    us.sort()
    u = us[0]
    bx = b / u^6
    Eu = EllipticCurve([f(ax), bx])
    assert Eu.is_isomorphic(E)
    return Eu

def check_curve(E):
    # a = -3 for efficient doubling https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html
    iso = get_iso(E, -3)
    if iso is None:
        return (E, "a = -3")
    E = iso
    b = E.a6()
    if b.is_square():
        return (E, "b is non-square")
    try:
        E.lift_x(1)
    except (ValueError):
        return (E, "x=1 on curve")
    return (E, None)

def to_signed(x):
    o = int(x.order())
    i = int(x)
    if i > o - 1000:
        return i - o
    return i

def curve_to_str(E, reason):
    a = to_signed(E.a4())
    b = E.a6()
    estr = "        y^2=x^3"
    if a == 1:
        estr += "+"
    elif a == -1:
        estr += "-"
    else:
        estr += "%+i*" % a
    estr += "x+%i" % b
    if reason:
        estr += ", failed condition: %s" % reason
    return estr + "\n"

def find_curves(d, q):
    H = hilbert_class_polynomial(d)
    Fp = GF(p)
    Fq = GF(q)
    Sp = H.roots(GF(p), False)
    Sq = H.roots(GF(q), False)
    # sorting for determinism
    Sp.sort()
    Sq.sort()
    zp = Fp(2)
    zq = Fq(2)
    # smallest non-squares in the fields
    while zp.is_square():
        zp += 1
    while zq.is_square():
        zq += 1
    output = "    Ep:\n" 
    for j0 in Sp:
        ap = 27*j0/(4*(1728-j0))
        Ep = EllipticCurve([ap, -ap])
        rp = Ep.count_points()
        if rp == q:
            (Ep, reason) = check_curve(Ep)
            output += curve_to_str(Ep, reason)
            if not reason:
                break
        # explicit twisting parameter for determinism
        Ep = Ep.quadratic_twist(zp)
        rp = Ep.count_points()
        if rp == q:
            (Ep, reason) = check_curve(Ep)
            output += curve_to_str(Ep, reason)
            if not reason:
                break
    output += "    Eq:\n"
    for j0 in Sq:
        aq = 27*j0/(4*(1728-j0))
        Eq = EllipticCurve([aq, -aq])
        rq = Eq.count_points()
        if rq == p:
            (Eq, reason) = check_curve(Eq)
            output += curve_to_str(Eq, reason)
            if not reason:
                break
        # explicit twisting parameter for determinism
        Eq = Eq.quadratic_twist(zq)
        rq = Eq.count_points()
        if rq == p:
            (Eq, reason) = check_curve(Eq)
            output += curve_to_str(Eq, reason)
            if not reason:
                break
    return output

pi_4 = (pi/4).numerical_approx()

def twist_security(p, q):
    # the number of points of the quadratic twist of Ep
    nt = 2 * (p + 1) - q
    # the largest prime order subgroup
    qt = factor(nt)[-1][0]
    # ECDLP security in bits = log2(sqrt(pi/4 * qt)) = log4(pi/4 * qt)
    return log(pi_4 * qt, 4)

def print_curves(d, q):
    reason = None
    if q % 4 != 3:
        reason = "q = 3 mod 4"
    elif 2^256 - 2*q >= 2^128:
        reason = "2^256 - 2*q < 2^128"
    elif twist_security(p, q) < 100:
        reason = "Ep twist security"
    elif twist_security(q, p) < 100:
        reason = "Eq twist security"
    if reason:
        sys.stdout.write("D = %i, q = %s, failed condition: %s\n" % (d, hex(q), reason))
        sys.stdout.flush()
        return False
    output = "D = %i, q = %s\n" % (d, hex(q))
    output += find_curves(d, q)
    sys.stdout.write(output)
    sys.stdout.flush()
    return True

def get_discriminants(wid, processes):
    # -D must be 3 mod 8
    for d in range(3 + 8 * wid, 1000000000, 8 * processes):
        yield d

def find_primes(wid, processes):
    for d in get_discriminants(wid, processes):
        (x, y) = cornacchia(d, 4*p)
        if x == None:
            continue
        q = p + 1 - x
        if is_prime(q):
            yield (-d, q)

def run_search(wid, processes):
    for (d, q) in find_primes(wid, processes):
        if print_curves(d, q):
            break

def worker(*args):
    try:
        run_search(*args)
    except (KeyboardInterrupt, SystemExit):
        pass
    except:
        print_exc()

processes = 1
print("p = %s" % hex(p))
print("Searching for curves using %i process(es)." % (processes,))

pool = Pool(processes=processes)

try:
    for wid in range(processes):
        pool.apply_async(worker, (wid, processes))

    pool.close()
    pool.join()
except (KeyboardInterrupt, SystemExit):
    pass
finally:
    pool.terminate()

tower_cycle_25519_output.txt

p = 0x7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffed
Searching for curves using 1 process(es).
D = -15203, q = 0x7ffffffffffffffffffffffffffffffecc1495732332411d29df513daa7d2cef, failed condition: 2^256 - 2*q < 2^128
D = -91515, q = 0x7fffffffffffffffffffffffffffffff9f6b957d81aa63f7041f01f802665879, failed condition: q = 3 mod 4
D = -124123, q = 0x7ffffffffffffffffffffffffffffffe9a8cf937fcfedc8f2ef44fdf1e7bc009, failed condition: q = 3 mod 4
D = -136827, q = 0x7ffffffffffffffffffffffffffffffecc1495732332411d29df513daa7d2cef, failed condition: 2^256 - 2*q < 2^128
D = -167995, q = 0x7fffffffffffffffffffffffffffffff286998bd82ef71f5ca89d6f16cb336b3, failed condition: 2^256 - 2*q < 2^128
D = -823635, q = 0x7fffffffffffffffffffffffffffffff9f6b957d81aa63f7041f01f802665879, failed condition: q = 3 mod 4
D = -1099643, q = 0x7fffffffffffffffffffffffffffffff96c281c3d13af039410ee81688034cb9, failed condition: q = 3 mod 4
D = -1117107, q = 0x7ffffffffffffffffffffffffffffffe9a8cf937fcfedc8f2ef44fdf1e7bc009, failed condition: q = 3 mod 4
D = -1231443, q = 0x7ffffffffffffffffffffffffffffffecc1495732332411d29df513daa7d2cef, failed condition: 2^256 - 2*q < 2^128
D = -1571315, q = 0x7fffffffffffffffffffffffffffffffbb2d703f7a1e51cdf484069145e5f0f7, failed condition: Ep twist security
D = -2504947, q = 0x7ffffffffffffffffffffffffffffffebc59fd79f85dd9f7f817bbb209a47933, failed condition: 2^256 - 2*q < 2^128
D = -3103075, q = 0x7ffffffffffffffffffffffffffffffe9a8cf937fcfedc8f2ef44fdf1e7bc009, failed condition: q = 3 mod 4
D = -3574995, q = 0x7ffffffffffffffffffffffffffffffe9a99b589414ec21bd5055787e2ebc05b, failed condition: 2^256 - 2*q < 2^128
D = -4059339, q = 0x7fffffffffffffffffffffffffffffff76ad8ac9b36c172dd4a1a2bfbd9068db, failed condition: 2^256 - 2*q < 2^128
D = -4871323, q = 0x7ffffffffffffffffffffffffffffffe963d636e78d179163e764dc747cfd2dd, failed condition: q = 3 mod 4
D = -6384387, q = 0x7ffffffffffffffffffffffffffffffefcdf3087b9118d9525b17fc4b23800ff, failed condition: 2^256 - 2*q < 2^128
D = -7194667, q = 0x7fffffffffffffffffffffffffffffff8e5e4ad5a38cf9f73dbc544f4eb2221f, failed condition: Eq twist security
D = -7412715, q = 0x7fffffffffffffffffffffffffffffff9f6b957d81aa63f7041f01f802665879, failed condition: q = 3 mod 4
D = -7857907, q = 0x7fffffffffffffffffffffffffffffffbf7f782cb7656b586eb6d2727927c79f, failed condition: Ep twist security
D = -9502827, q = 0x7ffffffffffffffffffffffffffffffeacaa2c51af9f7e5ea312752ff681364b, failed condition: 2^256 - 2*q < 2^128
D = -9896787, q = 0x7fffffffffffffffffffffffffffffff96c281c3d13af039410ee81688034cb9, failed condition: q = 3 mod 4
D = -14141835, q = 0x7fffffffffffffffffffffffffffffffbb2d703f7a1e51cdf484069145e5f0f7, failed condition: Ep twist security
D = -15578291, q = 0x7ffffffffffffffffffffffffffffffef50020d7e67c2d7897c09c8603552277, failed condition: 2^256 - 2*q < 2^128
D = -18522507, q = 0x7ffffffffffffffffffffffffffffffeac73c1c4a258cc42b7d533f0c77630f9, failed condition: q = 3 mod 4
D = -18923731, q = 0x7fffffffffffffffffffffffffffffff51a9d8fe03456636af2e650d284e38f5, failed condition: q = 3 mod 4
D = -22544523, q = 0x7ffffffffffffffffffffffffffffffebc59fd79f85dd9f7f817bbb209a47933, failed condition: 2^256 - 2*q < 2^128
D = -23684667, q = 0x7ffffffffffffffffffffffffffffffeb330afb1a411f124ba702a0f45fdd4eb, failed condition: 2^256 - 2*q < 2^128
D = -25165435, q = 0x7ffffffffffffffffffffffffffffffe98111cd8b9e9f4cd4df91dae02366eeb, failed condition: 2^256 - 2*q < 2^128
D = -25799635, q = 0x7fffffffffffffffffffffffffffffff069b3e0cf7b28e6e5b13282aae4498c9, failed condition: q = 3 mod 4
D = -27927675, q = 0x7ffffffffffffffffffffffffffffffe9a8cf937fcfedc8f2ef44fdf1e7bc009, failed condition: q = 3 mod 4
D = -31617403, q = 0x7fffffffffffffffffffffffffffffffee2dc09ab0631175c86423a9fe395b15, failed condition: q = 3 mod 4
D = -31750123, q = 0x7ffffffffffffffffffffffffffffffff735481d1969f317f9850b68df11df53
    Ep:
        y^2=x^3+1033607210507379184208924050302136195805720573571421383494254635561887108817*x+49627186934599064238114100101926864360189227744248910951774754919461467949413, failed condition: a = -3
        y^2=x^3+52505961432300074133204910020422701614252430600860146337254941544280186197162*x+5390083186358023578580582483921252312382561731960135682473850459676378622787, failed condition: a = -3
        y^2=x^3+8336692105860895464494377149111174415304366447201757372288740038389816954504*x+49098552390429031707615967815798512530835053088026505061147663700794594003866, failed condition: a = -3
        y^2=x^3+43631315466416890222200713544267428145437668999997332009045157888621573601741*x+14264729152241207489584778960076525781197323332822950010683634115334991218208, failed condition: a = -3
        y^2=x^3+10718014349005935679783726389847356234061618319152887801658469495776967913955*x+47178030269652162032001766114496597692573374013667394218070322508179596905994, failed condition: a = -3
        y^2=x^3+15624216191366223226488712870342128124010365689108705277049642762696634386781*x+48694404325044507323446774550294836787822051485591203842789233910296619365599, failed condition: a = -3
        y^2=x^3+12304709481512440019315335133214466853056011064801950722721680258419402622213*x+45591335137145657692470157371129487073578981268018331297007111745537162197736, failed condition: a = -3
        y^2=x^3+12150711221301539429012210309928491746408875114152044221976110498990616453811*x+18586399466903879991473302529259973881998983752424210263648700015988198009410, failed condition: a = -3
        y^2=x^3+20352987735221473516257277849241032471976275353924914706737612573969553035842*x+10864231974202505005298254719103602004094774167061528405285475420113270173111, failed condition: a = -3
        y^2=x^3+34089701995091023433675218039289855372260463600262882177709246863389648519990*x+23806342623567074278110274465054098554374528732557399842019545140566916299959, failed condition: a = -3
        y^2=x^3+32116662612214091576046505133538633668976789884294304346440691661547717112308*x+25779382006444006135738987370805320257658202448525977673288100342408847707641, failed condition: a = -3
        y^2=x^3-3*x+14608398204382016092671352939512451595537051853337674948011728275991088358317, failed condition: x=1 on curve
        y^2=x^3+9867099783571271650429190879713936506846088795282890315267408047308663357579*x+36855290968746022220137457970976415798501274303377441517318319629443822779266, failed condition: a = -3
        y^2=x^3+16204191650857772540614494264928546723067535026255455076116378408821663035906*x+44054600649112112810440523393603487995364696788417205450255348741296390172599, failed condition: a = -3
        y^2=x^3+30879249372722643630711594488646191194612294546316692679105973798574064064416*x+42446228111509339513234706612550240076276605293567868665796169631190311584417, failed condition: a = -3
        y^2=x^3+37582355755967460661262837839111590073155885029476606508354257789755814449893*x+46717421664208900980610252313171002974562873761108840051538689705692873320550, failed condition: a = -3
        y^2=x^3+47791527421039661649560821810851240653278166901108902033388035237247639654872*x+10104517197618436062224670693492713273356825431711379986340756766708925165077, failed condition: a = -3
        y^2=x^3+52726709650279579748837842120410762098846020134909805879396289587535698491933*x+41354679747028143703581203071465534622311777583283809122660019331366930624128, failed condition: a = -3
        y^2=x^3+47595175772201916250776332654020935473451823360651982217736017533930971776129*x+24510906152991353976287786298240193698830359444526116396213403756248179530611, failed condition: a = -3
        y^2=x^3-3*x+2854736054848353016025685075093639650597736603899477707703474608129419074168, failed condition: b is non-square
        y^2=x^3+54678392529107448115216101876365230082269297995195898318224187625924318876862*x+3217652089550649596569390627978723844365694337624383701504604378032245943087, failed condition: a = -3
        y^2=x^3+22022289704527479614372451450826688147915890399430568982714917181852234474880*x+35873754914130618097413041053517265778719101933389713037013874822104330345069, failed condition: a = -3
        y^2=x^3+55919354375465395964783524393576441683832024324427569755314571933379384260809*x+15813521945541613976015744886140097942423744067141698115313760564617444473120, failed condition: a = -3
        y^2=x^3+35451570591353544259084171314595070128353247539569817612243614033661532516949*x+22444474027304553452701321189748883798281744793250464407485177970295032303000, failed condition: a = -3
        y^2=x^3+34459757314974515773727607989556547101414454925608225659563187999008234507867*x+13802164573494362369106598605267392821859322259235604822138456027716948036809, failed condition: a = -3
        y^2=x^3-3*x+27074457538933191770719452708909629724744865620566577762023870439369894546176, failed condition: b is non-square
        y^2=x^3+42720603979093501958130287675444272946196921987672704908562370434262581823427*x+5611435879200570605670653622509539990234578095540052849873788549638734332278, failed condition: a = -3
        y^2=x^3-3*x+17539399157613241668790821773730052815957928319994112209279180327989982238224, failed condition: x=1 on curve
        y^2=x^3+42288458218875609718633299057253771931574167194833889568962102537380736779210*x+15607586399782487993152193447090181995060825137986392450766689466575828040739, failed condition: a = -3
        y^2=x^3+38288906519104567863334240884591260175440728274305192252177957172505819652135*x+19607138099553529848451251619752693751194264058515089767550834831450745167814, failed condition: a = -3
        y^2=x^3+42584302050570547471275687299971492171928616697106886701002945954936213690655*x+15311742568087550240509805204372461754706375635713395318725846049020351129294, failed condition: a = -3
        y^2=x^3-3*x+16418245240054761468360431793194228060262425048634775903619116753376469899862, failed condition: b is non-square
        y^2=x^3-3*x+27932929889965221102504386229895426583953429647519720236629462785080568835827, failed condition: b is non-square
        y^2=x^3+45452524228860808596935499280401001266280149316758665550083800186158426959076*x+12443520389797289114849993223942952660354843016061616469644991817798137860873, failed condition: a = -3
        y^2=x^3+1337866911313614356133230802470690144972085864169451782809069738527110761918*x+47193109328149182862719646084578432766858305419464667757256234095739678724605, failed condition: a = -3
        y^2=x^3+26715197177871855425353739275013832169023226344411763729153445695770755453569*x+31180847440786242286431753229330121757611765988408518290575346308185809366380, failed condition: a = -3
        y^2=x^3+46164764360747834331696853083910876986184642307350355203627214833148223882462*x+11731280257910263380088639420433076940450350025469926816101577170808340937487, failed condition: a = -3
        y^2=x^3-3*x+56715151266164652212524681749070172710143018146820880732158616339591971693524, failed condition: b is non-square
        y^2=x^3-3*x+13125004042814677155980591990878122966786258531669459881929988596999689589844, failed condition: b is non-square
        y^2=x^3+55002070084426747848378313661333392215363871088455637973604243056176826858297*x+2893974534231349863407178843010561711271121244364644046124548947779737961652, failed condition: a = -3
        y^2=x^3+3992806938126571536950788303011743475229559843260920508805985053535441170207*x+53903237680531526174834704201332210451405432489559361510922806950421123649742, failed condition: a = -3
        y^2=x^3-3*x+44092894459232229858783207215038906776675359108435248717845712561513390477054, failed condition: b is non-square
        y^2=x^3+54505904224929198267274812759134798008330274989949121353075844927219714907158*x+3390140393728899444510679745209155918304717342871160666652947076736849912791, failed condition: a = -3
        y^2=x^3+52097805020752008503155336716140645530189051389048558889868450210309397815687*x+46385916783248713669041246305626467171567527550173785038882734349177336034096, failed condition: a = -3
        y^2=x^3-3*x+47633078033327599743976807464873443075855076213564979039560011553108838026418, failed condition: b is non-square
        y^2=x^3-3*x+17523451383230374900436292617863907649717438939964238673872692863501483215968
    Eq:
        y^2=x^3+31526843997166882251594039438938336564635331629232624123151661316169450841186*x+37265471115955430546175147010213077116033876642292474559810757343481000245647, failed condition: a = -3
        y^2=x^3-3*x+15486530054027411065741275638098053507404864026203974095843587485538832614328, failed condition: b is non-square
        y^2=x^3+12734266304742372807687783311640788963746262521196466938881198097941319035336*x+45161778313915724904097709192703164962877043414034226570123611476626002359691, failed condition: a = -3
        y^2=x^3-3*x+37780485019954152974853906487973746184891557711560259611446846225047339702911, failed condition: b is non-square
        y^2=x^3-3*x+40377735979706220099746269888849985732466022427879643779527920385401382769984, failed condition: x=1 on curve
        y^2=x^3-3*x+45952012646381141163754966798223495239768269078591201635143546047706821570234, failed condition: b is non-square
        y^2=x^3-3*x+13356085066323316900905067377209210237441418757059735731047660514044026713231, failed condition: x=1 on curve
        y^2=x^3+27293028946492398781767729665545155524080903051685227691770188461386927940131*x+30603015672165698930017762838798798402542402883545465817234621113180393454896, failed condition: a = -3
        y^2=x^3-3*x+25675911719867737339625140396204798989996478324626569376465022644547366285284

Yet another update: as recommended by Veridise, we added twist security as one of the search conditions. So the final curve cycle that was selected has D = -31750123. Note that we don't expect to use the curves for anything that actually needs twist security.