// Schoof-Elkies-Atkin Algorithm!
// Mike Scott August 1999   mike@compapp.dcu.ie
// Counts points on GF(p) Elliptic Curve, y^2=x^3+Ax+B a prerequisite 
// for implemention of  Elliptic Curve Cryptography
// This program is intended as an aid to Cryptographers in generating 
// suitable curves (ideally with a prime number of points) with respect 
// to prime moduli from 160 to 512 bits in length. This can be done with
// relatively modest computing facilities - the average home computer and
// a little patience will suffice.
//
// An "ideal" curve is defined as one with with prime number of points.
//
// First the credits
//
// Basic algorithm is due to Schoof
// 1. "Elliptic Curves Over Finite Fields and the Computation of Square Roots 
//     mod p", Rene Schoof, Math. Comp., Vol. 44 pp 483-494
// 
// Elkies-Atkin ideas are described in
//
// 2. "Counting points on Elliptic Curves over Finite Fields", Rene Schoof, 
//    Jl. de Theorie des Nombres de Bordeaux 7 (1995) pp 219-254
//
// The particular variation implemented here is due to Mueller. See his thesis
//
// 3. "Ein Algorithmus zur Bestimmung der Punktanzahl elliptischer Kurven
//    uber endlichen Korpern der Charakteristik grosser drei",  
//    available from Volker Mueller's home page 
//    www.informatik.th-darmstadt.de/TI/Mitarbeiter/vmueller.html
//
// Other useful English-language publications are available from this site.
// Strongly recommended is the recent book 
//
// 4. "Elliptic Curves in Cryptography"
//     by Blake, Seroussi and Smart, London Mathematical Society Lecture Note 
//     Series 265, Cambridge University Press. ISBN 0 521 65374 6 
//
// Another useful reference is
// 5.  Elliptic Curve Public Key Cryptosystems", Menezes, 
//     Kluwer Academic Publishers, Chapter 7
//
// The Kangaroo continuation is due to Pollard:-
//
// 6. "Monte Carlo Methods for Index Computation"
//    by J.M. Pollard in Math. Comp. Vol. 32 1978 pp 918-924
//
// Fast FFT methods are largely as described by Shoup
//
// 7. "A New Polynomial Factorisation Algorithm and its implementation",
//     Victor Shoup, Jl. Symbolic Computation, 1996 
//
// A potentially more effective way of using Atkin primes is described in 
//
// 8. "Remarks on the Schoof-Elkies-Atkin Algorithm", L. Dewaghe, Math. Comp.
//    Vol. 67, 223, July 1998, pp 1247-1252
//
// Thanks are due to Richard Crandall for his encouragment, and the idea of 
// using projective co-ordinates.
//
// NOTE: Only for use on curves over large prime modulus P. 
// For smaller P use schoof.exe utility available from the same source.
//
// This first version does not process so-called Atkin primes
// Schoof's original algorithm is used for primes 3, 5 and 7 (this facilitates
// "early abort" using the -s option).
//
// After that only Elkies primes are used.  It is therefore not as fast as it 
// could be, particularly for smaller values of the prime modulus P. 
// However when the asymptotics kick-in, it starts to get competitive (when 
// you need it to). Since the average Cryptographer will only wish to 
// generate a few curves for practical use, this is deemed to be adequate. 
// The final continuation uses Pollard's Lambda ("kangaroo") algorithm.
//
// It is envisaged that from time-to-time this program will be modified
// and hopefully improved - that is speeded up.
// In particular it is planned to exploit Atkin Primes to allow two herds of 
// kangaroos complete the job more efficiently
//
// Asyptotically the algorithm should take time O(log(p)^5)
// However the Kangaroo continuation favours "smaller" curves, while
// the asymptotic behaviour is more accurate for "bigger" curves
//
// Timings in minutes:- random curves 180MHz Pentium Pro 
// (ignoring time to generate Modular Polynomials)
//
//               C1      C2      C3     Ave  Asymptotic multiplier wrt 160 bits
// 160-bit       2.5     3.0     2.0    2.5         1
// 192-bit       5.5     5.5     3.5    4.8         2.5
// 224-bit       9       7.5    10      8.8         5.4
// 256-bit      13.5    21.5    23     19.3        10.5
// 384-bit      86     108     120    105          60
// 512-bit     600     357     398    452         336
//
// As can be seen the asymptotic behaviour of the program would appear to 
// be about right. The wide variations in timing for the same size of curve
// is typical - it depends on how "lucky" you are finding Elkies Primes
//
// ****************
// Download Instructions
//
// To access the Windows 'NT/95/98 executables directly, point your
// browser at ftp://ftp.compapp.dcu.ie/pub/crypto, and download
//
// mueller.exe
// modpol.exe
// process.exe
// sea.exe
//
// The main program source file for these programs may be found in the
// same place, with .cpp extensions.
//
// To obtain the full source code first look at
// the README file on the ftp site, and then download
//
// ftp://ftp.compapp.dcu.ie/pub/crypto/miracl.zip
// To recompile, see the file sea.txt
//
// For more information:-
//
// http://indigo.ie/~mscott
//
// ****************
// Instructions for use
//
// First run the utility "mueller" to build up a collection of Modular 
// Polynomials. Each modular polynomial is associated with a small odd prime. 
// This needs to be done once only - ever, but you can from time 
// to time augment your collection of Polynomials by running it again. 
// Its quite time consuming, but in less than an hour you should have enough 
// to get started. The more you have, the bigger the prime modulus that you 
// can use.
//
// Then run the utility "process" to process the raw polynomial file with 
// respect to your chosen prime modulus P. This need to be done just once for
// every prime modulus of interest to you. This takes only a few minutes at 
// most.
//
// An alternative is to use instead the "modpol" application, which is a
// composite of "mueller" and "process". It directly generates the Modular 
// Polynomials reduced wrt a pre-specified prime modulus, suitable for 
// immediate use by "sea". If working with limited computing resources such 
// that sufficient generic Modular Polynomials cannot be generated by 
// "mueller", this may be your only alternative.
//
// Finally run this program "sea" specifying the A and B parameters of the 
// particular curve. This program can also search through many curves for 
// a curve ideal for cryptographic use (with a prime number of points).
//
// For example try:-
//
// mueller 0 120 -o mueller.raw
// process -f 65112*2#144-1 -i mueller.raw -o test160.pol
// sea -3 49 -i test160.pol
//
// Here the "mueller" program generates modular polynomials for all odd primes
// in the range 0 - 120 into the file mueller.raw. The "process" application
// reduces these 'raw' polynomials wrt the 160 bit prime modulus 
// P = 65112*2^144-1 to a file test160.pol. The "sea" application uses this 
// file to count the points on the curve Y^2 = X^3 - 3X + 49 mod P 
//
// Alternatively:-
//
// modpol -f 65112*2#144-1 0 120 -o test160.pol
// sea -3 49 -i test160.pol
//
// The number of Modular Polynomials required depends on the size of the 
// prime modulus P. It is also random in the sense that it depends on the 
// probability of each small prime being an "Elkies" prime wrt the given curve. 
// In the vast majority of cases the range suggested to "mueller" or "modpol" 
// should be 0 to bits(P), where bits(P) is the number of bits in P. However 
// you might get away with a much smaller value if you are lucky with your 
// "Elkies" primes. If modular polynomials could not be generated for all
// primes in the range, due to the use of the -s2, -s3 or -s6 flag in 
// "mueller" or "modpol" (see comments in mueller.cpp), then a somewhat 
// larger range might be needed.
//
// When using the "sea" program, the -s option is particularly useful
// and allows automatic search for an "ideal" curve. If a curve order is
// exactly divisible by a small prime, that curve is immediately abandoned, 
// and the program moves on to the next, incrementing the B parameter of 
// the curve. This is a fairly arbitrary but simple way of moving on to 
// the "next" curve. 
//
// Note that if a prime q is an Atkin prime, then we know at least that q 
// is not a factor of NP, in other words that NP mod q != 0.
// This can be easily proved.
//
// NOTE: The output file can be used directly with for example the ECDSA
// programs IF AND ONLY IF an ideal curve is found. If you wish to use
// a less-than-ideal curve, you will first have to factor NP completely, and
// find a random point of large prime order.
//
// ****************
//
// Rev. 1 September 1999  -  Faster and more Kangaroos!
// Rev. 2 October   1999  -  Poly class revamped. Faster search for tau mod lp
// Rev. 3 October   1999  -  Eliminated calculation of X^P
// Rev. 4 December  1999  -  Various optimizations
//
// This implementation is free. No terms, no conditions. It requires 
// version 4.24 or greater of the MIRACL library (a Shareware, Commercial 
// product, but free for non-profit making use), 
// available from ftp://ftp.compapp.dcu.ie/pub/crypto/miracl.zip 
//
// However this program may be freely used (unmodified!) to generate curves 
// for commercial use. It may be recompiled for this purpose on any hardware.
//
// 32-bit build only
//
// Note that is a little faster to use an integer-only build of MIRACL.
// See mirdef.hio
//
// Copyright Shamus Software Ltd. 1999
//

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cstring>
#include "ecn.h"         // Elliptic Curve Class
#include "crt.h"         // Chinese Remainder Theorem Class

//
// poly.h implements polynomial arithmetic. FFT methods are used for maximum
// speed, as the polynomials can get very big. 
// But all that gruesome detail is hidden away.
//
// polymod.h implements polynomial arithmetic wrt to a preset poynomial 
// modulus. This looks a bit neater. Function setmod() sets the modulus 
// to be used. Again fast FFT methods are used.
//
// polyxy.h implements a bivariate polynomial class
//

#include "poly.h"
#include "polymod.h"
#include "polyxy.h"

using namespace std;

#ifndef MR_NOFULLWIDTH
Miracl precision=18;            // max. 18x32 bits per big number
#else
Miracl precision(18,MAXBASE); 
#endif

PolyMod MY2,MY4;

ZZn A,B;         // Here ZZn are integers mod the prime p
                 // Montgomery representation is used internally

// Elliptic curve Point duplication formula

void elliptic_dup(PolyMod& X,PolyMod& Y,PolyMod& Z)
{ // (X,Y,Z)=2.(X,Y,Z)
    PolyMod W1,W2,W3,W4;

    W2=Z*Z;           // 1
    W3=A*(W2*W2);     // 2
    W1=X*X;           // 3
    W4=3*W1+W3;
    Z*=(2*Y);          // 4   Z has an implied y
    W2=MY2*(Y*Y);     // 5
    W3=4*X*W2;        // 6
    W1=W4*W4;         // 7
    X=W1-2*W3;
    W2*=W2;
    W2*=8;     // 8
    W3-=X;
    W3*=W4;         // 9  polynomial multiplications
    Y=W3-W2;    
    X*=MY2;          // fix up - move implied y from Z to Y
    Y*=MY2;
    Z*=MY2;
}

//
// This is addition formula for two distinct points on an elliptic curve
// Works with projective coordinates which are automatically reduced wrt a 
// polynomial modulus
// Remember that the expression for the Y coordinate of each point 
// (a function of X) is implicitly multiplied by Y.
// We know Y^2=X^3+AX+B, but we don't have an expression for Y
// So if Y^2 ever crops up - substitute for it 
//

void elliptic_add(PolyMod& XT,PolyMod& YT,PolyMod& ZT,PolyMod& X,PolyMod& Y)
{ // add (X,Y,1) to (XT,YT,ZT)
  // on an elliptic curve 
    PolyMod W1,W2,W4,W5,W6;

    W1=XT;
    W6=ZT*ZT;       // 1
    W4=X*W6;        // 2  * 
    W1-=W4;

    W2=YT;          // W2 has an implied y
    W6*=ZT;         // 3
    W5=Y*W6;        // 4  * W5 has an implied y 
    W2-=W5;
    if (iszero(W1))
    {
        if (iszero(W2)) 
        { // should have doubled
            elliptic_dup(XT,YT,ZT);
            return;
        }
        else
        { // point at infinity
            ZT.clear();
            return;    
        }
    }

    W4=W1+2*W4;     // W4=2*W4+W1 
    W5=W2+2*W5;     // W5=2*W5+W2

    ZT*=W1;       // 5

    W6=W1*W1;       // 6
    W1*=W6;         // 7
    W6*=W4;         // 8
    W4=MY2*(W2*W2);   // 9 Substitute for Y^2

    XT=W4-W6;

    W6=W6-2*XT;
    W2*=W6;       // 10
    W1*=W5;       // 11  polynomial multiplications
    W5=W2-W1;

    YT=W5/(ZZn)2;   

    return;
}

//
//   Program to compute the order of a point on an elliptic curve
//   using Pollard's lambda method for catching kangaroos. 
//
//   As a way of counting points on an elliptic curve, this
//   has complexity O(p^(1/4))
//
//   However Schoof puts a spring in the step of the kangaroos
//   allowing them to make bigger jumps, and lowering overall complexity
//   to O(p^(1/4)/sqrt(L)) where L is the product of the Schoof primes
//
//   See "Monte Carlo Methods for Index Computation"
//   by J.M. Pollard in Math. Comp. Vol. 32 1978 pp 918-924
//
//   This code has been considerably speeded up using ideas from
//   "Parallel Collision Search with Cryptographic Applications", van Oorchot 
//   and Wiener, J. Crypto., Vol. 12, 1-28, 1999
//

#define STORE 80
#define HERD 5

ECn wild[STORE],tame[STORE];
Big wdist[STORE],tdist[STORE];
int wname[STORE],tname[STORE];

Big kangaroo(Big p,Big order,Big ordermod)
{
    ECn ZERO,K[2*HERD],TE[2*HERD],X,P,G,table[128],trap;
    Big start[2*HERD],txc,wxc,mean,leaps,upper,lower,middle,a,b,x,y,n,w,t,nrp;
    int i,jj,j,m,sp,nw,nt,cw,ct,k,distinguished,nbits;
    Big D[2*HERD],s,distance[128],real_order;
    BOOL bad,collision,abort;
    forever
    {
// find a random point on the curve
        do
        {
            x=rand(p);
        } while (!P.set(x,x));

        lower=p+1-2*sqrt(p)-3; // lower limit of search
        upper=p+1+2*sqrt(p)+3; // upper limit of search

        w=1+(upper-lower)/ordermod;
        leaps=sqrt(w);
        mean=HERD*leaps/2;      // ideal mean for set S=1/2*w^(0.5)
        nbits=bits(leaps/16);
        if (nbits>30) nbits=30;
        distinguished=1<<nbits;
        for (s=1,m=1;;m++)
        { /* find table size */
            distance[m-1]=s*ordermod;
            s*=2;
            if ((2*s/m)>mean) break;
        }
        table[0]=ordermod*P;
        for (i=1;i<m;i++)
        { // double last entry
            table[i]=table[i-1];
            table[i]+=table[i-1];