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:

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. Let c = 2^256 - 2*q. Crandall's modular reduction algorithm [3] requires multiplications by c. We want c < 2^128, so c can fit into two machine words on 64-bit CPUs, making modular reductions faster.
4. We want both Ep and Eq to have the form of y^2 = x^3 - 3*x + b for some constant b, which allows for more efficient group law formulae [4].
5. We want the curve constant b to be a non-square number. This prevents the point x = 0 from being on the curve, which protects from some power analysis attacks [5].

Results

The attached script tower_cycle_25519.sage implements the search. The first discriminant value that produces a prime q that matches our conditions (1), (2) and (3) is D = -7857907 and the prime is:

q = 0x7fffffffffffffffffffffffffffffffbf7f782cb7656b586eb6d2727927c79f

The resulting curve cycle consists of two curves Ep and Eq, which we will call "Helios" and "Selene". They are the first curves that match our conditions (4) and (5).

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

Helios corresponds to the 11th root of the Hilbert class polynomial over GF(p) and Selene corresponds to the 5th root of the Hilbert class polynomial over GF(q).

We define the base point of each curve as the point with the smallest x-coordinate:

Helios: G = [3, 37760095087190773158272406437720879471285821656958791565335581949097084993268]
Selene: G = [1, 55227837453588766352929163364143300868577356225733378474337919561890377498066]

Security evaluation

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

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)/2, 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 -7857907.

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 ❌

Both Helios and Selene lack twist security, so they should not be used in single-coordinate ladders.

10. Completeness ❌

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

11. Indistinguishability ✔️

Both curves support indistinguishability using Elligator Squared [7].

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] Jacobian coordinates with a=-3 for short Weierstrass curves: https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html

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

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

[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(d, m):
    g = 1
    try:
        r = Mod(-d, m).sqrt(False)
    except (ValueError):
        for (p,e) in factor(m):
            e = e//2
            if e != 0:
                m //= p^(2*e)
                g *= p^e
        try:
            r = Mod(-d, m).sqrt(False)
        except (ValueError):
            return (None, "No sqrt mod m")
    r = int(r)
    if r > m//2:
       r = m - r
    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 (g*r, g*s)
    return (None, "sqrt(s2) not an integer")

# 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
    u = us[0]
    bx = b / u^6
    Eu = EllipticCurve([f(ax), bx])
    assert Eu.is_isomorphic(E)
    return Eu

def find_iso(E):
    # a = -3 for efficient doubling https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html
    iso = get_iso(E, -3)
    if iso is not None:
        return iso
    return E

def find_curves(d, q):
    H = hilbert_class_polynomial(d);
    Sq = H.roots(GF(q))
    Sp = H.roots(GF(p))
    Eqs = []
    Eps = []
    # check the first 12 roots
    num_roots = 12
    for i in range(num_roots):
        j0 = Sq[i][0]
        aq = 27*j0/(4*(1728-j0))
        Eq = EllipticCurve([aq, -aq])
        rq = Eq.count_points()
        if rq == p:
            Eqs.append(Eq)
        Eq = Eq.quadratic_twist()
        rq = Eq.count_points()
        if rq == p:
            Eqs.append(Eq)
    assert len(Eqs) == num_roots
    for i in range(num_roots):
        j0 = Sp[i][0]
        ap = 27*j0/(4*(1728-j0))
        Ep = EllipticCurve([ap, -ap])
        rp = Ep.count_points()
        if rp == q:
            Eps.append(Ep)
        Ep = Ep.quadratic_twist()
        rp = Ep.count_points()
        if rp == q:
            Eps.append(Ep)
    assert len(Eps) == num_roots
    return [Eps, Eqs]

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):
    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 not b.is_square():
        estr += " (b is non-square)"
    return estr

def print_curves(d, q):
    output = "D = %i, q = %s" % (d, hex(q))
    if q % 4 == 3:
        output += " (q = 3 mod 4)"
    if 2^256 - 2*q < 2^128:
        output += " (2^256 - 2*q < 2^128)"
    output += "\n"
    (Eps, Eqs) = find_curves(d, q)
    output += "    Ep:\n" 
    for Ep in Eps:
        Ep = find_iso(Ep)
        output += "        %s\n" % curve_to_str(Ep)
    output += "    Eq:\n"
    for Eq in Eqs:
        Eq = find_iso(Eq)
        output += "        %s\n" % curve_to_str(Eq)
    sys.stdout.write(output)
    sys.stdout.flush()

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):
        print_curves(d, q)

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

processes = cpu_count()
print("p = %s" % hex(p))
print("Searching for curves using %i processes." % (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()

https://github.com/kayabaNerve/fcmp-plus-plus/tree/develop/crypto/helioselene is a Rust implementation of these two curves using constant time, yet generic, arithmetic. It uses the complete Weierstrass addition formulas (and I have yet to implement dedicated doubling). It's presumed not performant enough for production before tailored arithmetic is implemented.