eccons.cc

// (C) 2020 Riad S. Wahby <rsw@cs.stanford.edu>
//
// This file is part of eccons.
//
// Licensed under the Apache License, Version 2.0 (see
// LICENSE or https://www.apache.org/licenses/LICENSE-2.0).
// This file may not be copied, modified, or distributed
// except according to those terms.

// eccons: elliptic curve construction toolkit

// on input q, find (p, D) s.t. F_p has an elliptic curve of order q with discriminant D
//
// This is an implementation of Alg 2.2, steps 1-3a plus the primality testing of 3b,
// from Br\"{o}ker and Stevenhagen, "Constructing Elliptic Curves of Prime Order."
//     https://arxiv.org/abs/0712.2022
//
// I have extended this algorithm to support finding curves of order 4*q, q prime by
// specializing the algorithm of Br\"{o}ker and Stevenhagen, "Efficient CM-Constructions
// of Elliptic Curves over Finite Fields." Math. Comp. vol 76 no 260, Oct. 2007.
//     https://www.jstor.org/stable/40234483
//
// The above specialization leverages the fact that the algorithm for finding curves
// of arbitrary order N decays to almost exactly the algorithm for finding curves of
// prime order q when N=4q. Specifically, to find curves for order 4*q, run the
// prime-order algorithm for order q, but for the primality test in step 3b, rather
// than test p = q + 1 +/- x (for x a solution to x^2 + (-D)y^2 = 4*q), one should
// test p' = 4*q + 1 +/- 2*x. Then if p' is prime, then F_p' has a curve of order N
// with j-invariant a root of H_D.
//
// When constructing elliptic curves with specified embedding degree, this program
// uses the Cocks-Pinch method. See Boneh, Rubin, and Silverberg, "Finding composite
// order ordinary elliptic curves using the Cocks-Pinch method." J. Number Theory
// vol 131 issue 5, May 2011.
//     https://www.sciencedirect.com/science/article/pii/S0022314X10001344
// See also Freeman, "Methods for constructing pairing-friendly elliptic curves,"
// 10th Workshop on Elliptic Curve Cryptography, Sept. 2006.
//     http://cacr.uwaterloo.ca/conferences/2006/ecc2006/freeman.pdf
//
// When constructing Montgomery curves from the Weierstrass form, this program uses
// the work of Okeya, Kurumatani, Sakurai, "Elliptic curves with the Montgomery-form
// and their cryptographic applications," PKC, January 2000.
//     https://link.springer.com/chapter/10.1007/978-3-540-46588-1_17
// The Edwards form follows from the birational equivalence given by Bernstein,
// Birkner, Joye, Lange, and Peters, "Twisted Edwards Curves," AFRICACRYPT, June 2008.
//     https://eprint.iacr.org/2008/013.pdf
// The authors also give the completeness condition for twisted Edwards curves.
//
// Edwards curves constructed as described above are not complete. To construct a
// complete curve, this program uses work from Morain, "Edwards curves and CM curves."
//     https://arxiv.org/abs/0904.2243
// The algorithm used here to trace a descending path into the isogeny volcano is
// slightly adapted from the one given by Fouquet and Morain, "Isogeny Volcanoes
// and the SEA Algorithm," ANTS, July 2002.
//     https://link.springer.com/chapter/10.1007/3-540-45455-1_23
//
// Finally, the alternative form of the Edwards curve that this program computes
// (when possible) is the one suggested by Hisil, Wong, Carter, and Dawson,
// "Twisted Edwards curves revisited," ASIACRYPT, Dec 2008.
//     https://link.springer.com/chapter/10.1007/978-3-540-89255-7_20
// It has a fast point addition algorithm when using extended coordinates, which
// the above paper also introduces.

#include <array>
#include <cstdlib>
#include <cstring>
#include <gmp.h>
#include <gmpxx.h>
#include <iostream>
#include <random>
#include <stack>
#include <tuple>
#include <vector>

unsigned get_two_adicity(const mpz_class &p);

bool maybe_print_entry(const mpz_class &D, mpz_class &p, long two_adicity, bool found) {
    unsigned count = two_adicity ? get_two_adicity(p) : 0;
    if (count < two_adicity) {
        return found;
    }

    if (!found) {       // first one, so print out column header
        gmp_printf("DPvalues = [\n#     =D=        =p=\n");
    }

    gmp_printf("    [%Zd, %Zd],", D.get_mpz_t(), p.get_mpz_t());
    if (two_adicity) {
        printf("    # 2-adicity: %d\n", count);
    } else {
        printf("\n");
    }

    return true;
}

// modular square roots via Tonelli-Shanks
int mpz_sqrtm(mpz_t q, const mpz_t n, const mpz_t p) {
    mpz_class w_c, y_c, y_save_c;
    auto *w = w_c.get_mpz_t();
    auto *y = y_c.get_mpz_t();
    auto *y_save = y_save_c.get_mpz_t();
    unsigned int i = 0, s = 0;

    // if n is a multiple of p, then sqrt is 0
    if (mpz_divisible_p(n, p)) {
        mpz_set_ui(q, 0);
        return 1;
    }

    // if p=3 mod 4, q = n ^ ((p+1) / 4) mod p
    if(mpz_tstbit(p, 1) == 1) {
        mpz_set(q, p);
        mpz_add_ui(q, q, 1);
        mpz_fdiv_q_2exp(q, q, 2);
        mpz_powm(q, n, q, p);
        goto return_principle_root;
    }

    // factor out 2^s from p-1
    mpz_set(q, p);
    mpz_sub_ui(q, q, 1);
    for (s = 0; mpz_tstbit(q, s) == 0; s++);
    mpz_fdiv_q_2exp(q, q, s);

    // Search for a non-residue mod p by picking the first w s.t. (w/p) = -1
    for (mpz_set_ui(w, 2); mpz_legendre(w, p) != -1; mpz_add_ui(w, w, 1));

    // w = w^Q mod p
    mpz_powm(w, w, q, p);
    // y = n^Q mod p
    mpz_powm(y, n, q, p);
    // q = n^((Q+1)/2) mod p
    mpz_add_ui(q, q, 1);
    mpz_fdiv_q_2exp(q, q, 1);
    mpz_powm(q, n, q, p);

    // Tonelli-Shanks main loop
    while (true) {
        // find the order of y mod p
        // if y = 1 mod p, then we've found the sqrt
        mpz_set(y_save, y);
        for (i = 0; (i < s) && (mpz_cmp_ui(y, 1) != 0); i++, mpz_powm_ui(y, y, 2, p));
        if (i == 0) {
            goto return_principle_root;
        } else if (i == s) {
            return 0;
        }

        // otherwise, let t = w^(2^(s-i-1))
        mpz_powm_ui(w, w, 1 << (s-i-1), p);

        // update values for next iteration
        s = i;
        // q = q * w
        mpz_mul(q, q, w);
        mpz_mod(q, q, p);
        // w = w^2
        mpz_mul(w, w, w);
        mpz_mod(w, w, p);
        // y = y w
        mpz_mul(y, y_save, w);
        mpz_mod(y, y, p);
    }

return_principle_root:
    // let q be the principal root
    mpz_fdiv_q_2exp(y, p, 1);
    if (mpz_cmp(q, y) == 1) {
        mpz_sub(q, p, q);
    }
    // check that q is indeed the sqrt
    mpz_mul(y, q, q);
    mpz_mod(y, y, p);
    mpz_mod(w, n, p);
    if (mpz_cmp(y, w) != 0) {
        return 0;
    }
    return 1;
}

// we return x s.t. x^2 + |D|y^2 = 4q
//
// This follows section 1.5.3 of Cohen, "A Graduate Course in Computational Number Theory"
int cornacchia(mpz_class &x, const mpz_class &q, const mpz_class &D, const mpz_class &sqD) {
    // variables
    mpz_class tmp;
    mpz_class a;
    mpz_class b;
    mpz_class l;
    auto *mtmp = tmp.get_mpz_t();
    auto *ml = l.get_mpz_t();
    const auto *mD = D.get_mpz_t();

    // set up temporary variables
    a = 2 * q;                          // r_0 = N = 2 * q
    b = sqD;                            // r_1 = sqrt(D) - NOTE: main() ensures that sqD is sqrt(D) mod 4q
    l = 4 * q;
    mpz_sqrt(ml, ml);                   // floor(sqrt(4*q))

    // Euclidean algorithm
    while (b > l) {
        tmp = a % b;
        a = b;
        b = tmp;
    }

    // primitive solutions have (4q - r_k^2) divisible by D and square
    tmp = (4 * q) - (b * b);        // 4*q - r_k^2
    mpz_abs(ml, mD);
    if (mpz_divisible_p(mtmp, ml)) {
        mpz_divexact(mtmp, mtmp, ml);
        if (mpz_perfect_square_p(mtmp)) {
            x = b;
            return 1;
        }
    }

    return 0;
}

void convert_next_value(char *current, unsigned &len, char *&next, char &sgn) {
    // find first non-numeric character
    for (; (len > 0) && (current[0] >= '0') && (current[0] <= '9'); len--, current++);
    sgn = current[0];
    // the below is safe even when len == 0 because len was calculated by strlen
    current[0] = '\0';
    next = current;
}

bool convert_formula (char *current, mpz_class &q) {
    unsigned len = std::strlen(current);
    char *next = nullptr;
    long tmp_l = 0;
    bool op_is_add = true;
    char sgn = '+';
    mpz_class tmp_z;
    mpz_class one = 1;

    q = 0;
    while (len > 0) {
        // grab either 2^x or x
        tmp_z = 0;
        if ((len > 1) && (std::strncmp("2^", current, 2) == 0)) {
            // found 2^, now get the exponent
            current += 2;
            len -= 2;
            if (len == 0) {
                return false;
            }
            convert_next_value(current, len, next, sgn);
            tmp_l = std::strtol(current, nullptr, 10);
            next[0] = sgn;      // (for error reporting)
            current = next;

            // add or subtract 2^exponent
            mpz_mul_2exp(tmp_z.get_mpz_t(), one.get_mpz_t(), tmp_l);
        } else {
            convert_next_value(current, len, next, sgn);
            tmp_z.set_str(current, 0);
            next[0] = sgn;      // (for error reporting)
            current = next;
        }

        // do the operation
        if (op_is_add) {
            q += tmp_z;
        } else {
            q -= tmp_z;
        }

        // grab the next sign, if there is one
        if (len > 0) {
            if (sgn == '+') {
                op_is_add = true;
            } else if (sgn == '-') {
                op_is_add = false;
            } else {
                return false;
            }
            current += 1;
            len -= 1;
        }
    }
    return true;
}

unsigned get_two_adicity(const mpz_class &p) {
    unsigned count = 0;
    mpz_class pm1 = p - 1;
    auto *pm1t = pm1.get_mpz_t();
    while (mpz_divisible_ui_p(pm1t, 2)) {
        count++;
        mpz_divexact_ui(pm1t, pm1t, 2);
    }

    return count;
}

void print_usage (const char *name) {
    std::cout << "Usage:" << std::endl
              << "  " << name << " <q>          [options]" << std::endl
              << R"EOF(
You can specify the target either as an integer or as a formula, e.g., '2^255-19'.

  option      description                                       default 
  --          --                                                --
  -r <max_r>  sets how large a D to search for                  6

  -D <maxD>   Do not consider discriminants |D| > maxD.         ((r+1)*log q)^2

  -k <deg>    Search for a curve with embedding degree <deg>    unconstrained
              via Cocks-Pinch method. Cannot be used with -4.

  -R <seed>   Use <seed> to search for candidate trace values   random
              when -k is set.

  -4          Search for a curve of order 4*q rather than q.    false
              Cannot be used with -k.

  -E          Construct complete twisted Edwards curves only.   false
              Implies -s; should be used with -4 or -k.

  -2 <min>    Print 2-adicity of each identified p; reject      unconstrained
              any p with less than <min> 2-adicity.

  -j <nproc>  Use at most <nproc> parallel processes.           # of cpus

Example:
  )EOF"
              << name << " 2^251-9 -r 10 -s -4" << std::endl;
}

void check_nargs(int nargs, int expect, const char *name, const char* argname) {
    if (nargs < expect) {
        print_usage(name);
        std::cout << std::endl << argname << " requires an argument!" << std::endl;
        std::exit(-1);
    }
}

int main (int argc, char **argv) {
    if (argc < 2) {
        print_usage(argv[0]);
        std::exit(-1);
    }

    // STEP 0: parse argumnts
    bool times_four = false;
    bool compl_edw = false;
    int arg_curr = 1;
    long dmax_user = 0;
    long k_target = 0;
    long r = 6;
    long two_adicity = 0;
    long rseed_user = 0;
    long max_nproc = 0;
    mpz_class q;
    while (arg_curr < argc) {
        if (std::strncmp("-R", argv[arg_curr], 2) == 0) {
            check_nargs(argc - arg_curr, 2, argv[0], "-R");
            rseed_user = std::strtol(argv[arg_curr+1], nullptr, 10);
            arg_curr += 2;
        } else if (std::strncmp("-D", argv[arg_curr], 2) == 0) {
            check_nargs(argc - arg_curr, 2, argv[0], "-D");
            dmax_user = std::strtol(argv[arg_curr+1], nullptr, 10);
            arg_curr += 2;
        } else if (std::strncmp("-k", argv[arg_curr], 2) == 0) {
            check_nargs(argc - arg_curr, 2, argv[0], "-k");
            k_target = std::strtol(argv[arg_curr+1], nullptr, 10);
            arg_curr += 2;
        } else if (std::strncmp("-4", argv[arg_curr], 2) == 0) {
            times_four = true;
            arg_curr += 1;
        } else if (std::strncmp("-E", argv[arg_curr], 2) == 0) {
            compl_edw = true;
            arg_curr += 1;
        } else if (std::strncmp("-2", argv[arg_curr], 2) == 0) {
            check_nargs(argc - arg_curr, 2, argv[0], "-2");
            two_adicity = std::strtol(argv[arg_curr+1], nullptr, 10);
            arg_curr += 2;
        } else if (std::strncmp("-r", argv[arg_curr], 2) == 0) {
            check_nargs(argc - arg_curr, 2, argv[0], "-r");
            r = std::strtol(argv[arg_curr+1], nullptr, 10);
            arg_curr += 2;
        } else if (std::strncmp("-j", argv[arg_curr], 2) == 0) {
            check_nargs(argc - arg_curr, 2, argv[0], "-j");
            max_nproc = std::strtol(argv[arg_curr+1], nullptr, 10);
            arg_curr += 2;
        } else if (argv[arg_curr][0] == '-') {
            print_usage(argv[0]);
            std::cout << std::endl << "Unrecognized option " << argv[arg_curr] << std::endl;
            std::exit(-1);
        } else {
            if (!convert_formula(argv[arg_curr], q)) {
                print_usage(argv[0]);
                std::cout << std::endl << "Failed to convert '" << argv[arg_curr] << "' to an integer." << std::endl;
                std::exit(-1);
            }
            arg_curr += 1;
        }
    }

    // sanity checks and preparatory stuff
    q = abs(q);
    k_target = std::abs(k_target);
    r = std::abs(r);
    if (rseed_user) {
        rseed_user = std::abs(rseed_user);
    } else if (k_target != 0) {
        std::random_device rd;
        rseed_user = rd();
    }
    dmax_user = std::abs(dmax_user);
    if (q <= 3 || !mpz_probab_prime_p(q.get_mpz_t(), 128)) {
        print_usage(argv[0]);
        std::cout << std::endl << "This program requires q > 3 and prime." << std::endl;
        std::exit(-1);
    } else if (k_target != 0) {
        constexpr const char *const ERR1 = "This program does not support -4 and -k at the same time, but often the\nresults from -k will have cofactor divisible by four (hint: try -R).";
        constexpr const char *const ERR2 = "This program requires q to be congruent to 1 mod k.";
        const char *err = nullptr;
        if (times_four) {
            err = ERR1;
        } else if (q % k_target != 1) {
            err = ERR2;
        }
        if (err != nullptr) {
            print_usage(argv[0]);
            std::cout << std::endl << err << std::endl;
            std::exit(-1);
        }
    }
    mpz_class N = q * (times_four ? 4 : 1);
    gmp_printf("N = %Zd\n", N.get_mpz_t());

    // now begins the algorithm
    using DTableEnt = std::array<mpz_class, 2>;
    using DTable = std::vector<DTableEnt>;
    DTable dtab{};

    unsigned log_q = mpz_sizeinbase(q.get_mpz_t(), 2);
    mpz_class p{3};
    mpz_class tmp_z;
    mpz_class tmp_z2;
    mpz_class tmp_z3;
    mpz_class kq_div;
    if (k_target != 0) {
        kq_div = (q - 1) / k_target;
        printf("# random seed was %ld\n", rseed_user);
    }
    gmp_randclass rand(gmp_randinit_default);
    rand.seed(rseed_user);

    // STEP 1: Create a table of (D, sqrt(D)) candidates
    while (p < (r+1) * log_q) {
        if (mpz_legendre(q.get_mpz_t(), p.get_mpz_t()) == 1) {
            if ((p-1) % 4 == 0) {
                tmp_z2 = p;
            } else {
                tmp_z2 = -p;
            }

            if (mpz_sqrtm(tmp_z.get_mpz_t(), tmp_z2.get_mpz_t(), q.get_mpz_t()) != 1) {
                gmp_printf("# ERROR detected computing sqrt(%Zd) mod %Zd\n", tmp_z2.get_mpz_t(), q.get_mpz_t());
                std::exit(-1);
            } else {
                dtab.emplace_back(DTableEnt{tmp_z2, tmp_z});
            }
        }

        // go to the next prime
        mpz_nextprime(p.get_mpz_t(), p.get_mpz_t());
    }

    // STEP 2: Compute products of the table elements to find D and sqrt(D)
    // max D we will consider is ((r+1)log_q)^2
    mpz_class x;
    mpz_class Dmax;
    if (dmax_user == 0) {
        Dmax = (r+1) * log_q;
        Dmax *= Dmax;
    } else {
        Dmax = dmax_user;
    }

    // walk the D candidates, finding candidates
    using DStackEnt = std::tuple<size_t, DTableEnt>;
    using DStack = std::stack<DStackEnt>;
    DStack dstk{};

    dstk.push(DStackEnt{0, {1, 1}});
    auto last_ent = dtab.size();
    bool found = false;

    while (!dstk.empty()) {
        DStackEnt curr = std::move(dstk.top());
        dstk.pop();
        auto next = std::get<0>(curr);
        auto& Dent = std::get<1>(curr);

        while (true) {
            // if we've run out of dtab, give up
            if (next >= last_ent) {
                break;
            }

            // if we've gotten too big, give up
            tmp_z = Dent[0] * dtab[next][0];
            tmp_z = abs(tmp_z);
            if (tmp_z > Dmax) {
                break;
            }

            // if there are two paths, push one onto the stack
            if (next + 1 < last_ent) {
                dstk.push(DStackEnt{next+1, Dent});
            }

            // extend current path
            Dent[0] *= dtab[next][0];
            Dent[1] = (Dent[1] * dtab[next][1]) % q;
            next++;

            // for order q or 4q, unconstrained k, check that candidate D is -3 mod 8
            if ((k_target == 0) && (Dent[0] % 8 == -3)) {
                // make sure that we have the odd sqrt of D s.t. it's also a sqrt mod 4q
                if (mpz_even_p(Dent[1].get_mpz_t())) {
                    Dent[1] = q - Dent[1];
                }

                // STEP 3: run Cornacchia's algorithm; if q + 1 +/- x is prime, we've found p
                //         If we're looking for a curve of order 4q, check 4q + 1 +/- 2x instead.
                if (cornacchia(x, q, Dent[0], Dent[1])) {
                    for (unsigned i = 0; i < 2; i++) {
                        p = N + 1 + (times_four ? 2 : 1) * (i ? x : -x);
                        if (mpz_probab_prime_p(p.get_mpz_t(), 128) > 0) {
                            found = maybe_print_entry(Dent[0], p, two_adicity, found);
                        }
                    }
                }
            // for constrained k, just need D < 0 and =0 or =1 mod 4
            } else if ((k_target != 0) && (Dent[0] < 0) && ((Dent[0] % 4 == 0) || Dent[0] % 4 == -3)) {
                unsigned maxrep = log_q * log_q;
                for (unsigned i = 0; i < maxrep; i++) {
                    mpz_invert(tmp_z3.get_mpz_t(), Dent[1].get_mpz_t(), q.get_mpz_t());
                    do {
                        // X of order k mod q
                        tmp_z2 = rand.get_z_range(q);
                        mpz_powm_ui(tmp_z.get_mpz_t(), tmp_z2.get_mpz_t(), k_target, q.get_mpz_t());
                    } while (tmp_z == 1);
                    for (unsigned j = 0; j < 2; j++) {
                        // Y = (X-1) * sqrt(-D) mod q
                        tmp_z2 = (tmp_z - 1);
                        if (j) {
                            tmp_z2 *= -tmp_z3;
                        } else {
                            tmp_z2 *= tmp_z3;
                        }
                        tmp_z2 %= q;

                        // p = ((X+1)^2 + DY^2)
                        p = tmp_z + 1;
                        p *= p;
                        tmp_z2 *= tmp_z2;
                        tmp_z2 *= -Dent[0];
                        p += tmp_z2;

                        // make sure p is divisible by 4
                        if (!mpz_divisible_ui_p(p.get_mpz_t(), 4)) {
                            continue;
                        }
                        mpz_divexact_ui(p.get_mpz_t(), p.get_mpz_t(), 4);

                        if (mpz_probab_prime_p(p.get_mpz_t(), 128) > 0) {
                            found = maybe_print_entry(Dent[0], p, two_adicity, found);
                            if (found) {
                                // break two levels
                                i = maxrep;
                                break;
                            }
                        }
                    }
                }
            }
        }
    }

    // STEP 4: Print out info for the user.
    if (!found) {
        std::cout << "# Sorry! I didn't find any candidates. Try running with higher -r."
                  << std::endl;
    } else {
        std::cout << ']'
                  << std::endl
                  << "# change this value to increase or decrease number of parallel threads"
                  << std::endl
                  << "NPROC = ";
        if (max_nproc) {
            std::cout << max_nproc;
        } else {
            std::cout << "mp.cpu_count()";
        }
        std::cout << std::endl
                  << "# when True, generate only complete twisted Edwards curves"
                  << std::endl
                  << "COMPL_EDW = "
                  << (compl_edw ? "True" : "False")
                  << std::endl;
    }
}