schoof.cpp

miracl-大数运算库,大家使用有什么问题请多多提意见

// Schoof's Original Algorithm!
// Mike Scott June 1999   [e-mail]
// Counts points on GF(p) Elliptic Curve, y^2=x^3+Ax+B a prerequisite 
// for implemention of  Elliptic Curve Cryptography
//
// !!!!!!!!!!!!!!!!!!!!!!!!
// NOTE! September 1999. This program has been rendered effectively obsolete 
// by the implementation of the superior Schoof-Elkies-Atkin Algorithm
// This has O(log(p)^5) complexity compared with O(log(p)^6), and so will
// always be faster.
//
// Read the comments at the head of the file 
// ftp://ftp.computing.dcu.ie/pub/crypto/sea.cpp
// for more details
//
// However this program is self-contained and hence easier to use.
// It is quite adequate for counting points on up to 256 bit curves (that is 
// for primes p < 256 bits in length) assuming a reasonably powerful PC
// It also works for the smallest curves, and so is ideal for educational
// use
//
// It is now straightforward to create an executable that runs under Linux.
// See ftp://ftp.computing.dcu.ie/pub/crypto/README for details
//
// For elliptic curves defined over GF(2^m), see the utility schoof2.exe
//
// !!!!!!!!!!!!!!!!!!!!!!!!
//
// The self-contained Windows Command Prompt executable for this program may 
// obtained from ftp://ftp.computing.dcu.ie/pub/crypto/schoof.exe 
//
// Basic algorithm is due to Schoof
// "Elliptic Curves Over Finite Fields and the Computation of Square Roots 
// mod p", Rene Schoof, Math. Comp., Vol. 44 pp 483-494
// Expression for Mod 2 Cardinality from "Counting Points on Elliptic
// Curves over Finite Fields", Rene Schoof, Jl. de Theorie des Nombres
// de Bordeaux 7 (1995) pp 219-254
// Another useful reference is
// "Elliptic Curve Public Key Cryptosystems", Menezes, 
// Kluwer Academic Publishers, Chapter 7
// Thanks are due to Richard Crandall for the tip about using prime powers
//
// **
// This program implements Schoof's original algorithm, augmented by
// the use of prime powers. By finding the Number of Points mod the product
// of many small primes and large prime powers, the final search for NP is 
// greatly speeded up.
// Final phase search uses Pollard Lambda ("kangaroo") algorithm
// This final phase effectively stretches the range of Schoof's
// algorithm by about 70-80 bits.
// This approach is only feasible due to the use of fast FFT methods
// for large degree polynomial multiplication
// **
//
// Ref "Monte Carlo Methods for Index Computation"
// by J.M. Pollard in Math. Comp. Vol. 32 1978 pp 918-924
// Ref "A New Polynomial Factorisation Algorithm and its implementation",
// Victor Shoup, Jl. Symbolic Computation, 1996 
//
//
// An "ideal" curve is defined as one with with prime number of points.
//
// When using the "schoof" 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: 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.
//
// 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.computing.dcu.ie/pub/crypto/miracl.zip 
//
// However this program may be used (unmodified) to generate curves for 
// commercial use.
//
// 32-bit build only
//
// Revision history
//
// Rev. 1 June 1999       - included prime powers
// Rev. 2 June 1999       - tweaked some inner loops
// Rev. 3 July 1999       - changed from rational to projective coordinates
//                          power windowing helps X^PP, Y^PP calculations
// Rev. 4 August 1999     - suppressed unnecessary creation of ZZn temporaries 
//                          in poly.cpp and polymod.cpp
// Rev. 5 September 1999  - Use fast modular composition to calculate X^PP 
//                          Half required size of ZZn's
//                          More & Faster Kangaroos!   
// Rev. 6 October 1999    - Revamped Poly class
//                          25% Faster technique for finding tau mod lp
// Rev. 7 November 1999   - Calculation of Y^PP eliminated!
//
// Rev. 8 November 1999   - Find Y^PP using composition - faster tau search
// Rev. 9 December 1999   - Various optimisations
// Rev. 10 March 2000     - Eigenvalue Heuristic introduced
//
// Timings for test curve Y^2=X^3-3X+49 mod 65112*2^144-1 (160 bits)
// Pentium Pro 180MHz
//
// Rev. 0  - 100 minutes
// Rev. 1  - 67 minutes
// Rev. 2  - 57 minutes
// Rev. 3  - 51 minutes
// Rev. 4  - 46 minutes
// Rev. 5  - 31 minutes
// Rev. 6  - 25 minutes
// Rev. 7  - 22 minutes
// Rev. 8  - 21 minutes
// Rev. 9  - 18 minutes
// Rev. 10 - 13 minutes
//
// 160 bit curve -  13 minutes
// 192 bit curve -  60 minutes
// 224 bit curve - 176 minutes
// 256 bit curve - 355 minutes
//
// This execution time can be related directly to 
// the O(log(p)^6) complexity of the algorithm as implemented here.
//
// Note that a small speed-up can be obtained by using 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 get very big. For example when searching for the
// curve cardinality mod the prime 31, the polynomials are of degree (31*31-1)/2
// = 480. 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.

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

using namespace std;

#ifndef MR_NOFULLWIDTH
Miracl precision=10;            // max. 10x32 bits per big number
#else
Miracl precision(10,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
    if (A==(ZZn)(-3))
    {
        W4=3*(X-W2)*(X+W2); // 2  (and save 1)  
    }
    else
    {
        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;           // 8
    W2*=8;
    W3-=X;
    W3*=W4;          // 9
    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,PolyMod& Z)
{ // add (X,Y,Z) to (XT,YT,ZT) on an elliptic curve 

    PolyMod W1,W2,W3,W4,W5,W6;

    if (!isone(Z))
    { // slower if Z!=1
        W3=Z*Z;    // 1x
        W1=XT*W3;  // 2x
        W3*=Z;     // 3x
        W2=YT*W3;  // 4x
    }
    else
    {
        W1=XT;
        W2=YT;      // W2 has an implied y
    }
    W6=ZT*ZT;       // 1
    W4=X*W6;        // 2
    W6*=ZT;         // 3
    W5=Y*W6;        // 4 W5 has an implied y

    W1-=W4;
    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;     
    W5=W2+2*W5;     

    ZT*=W1;       // 5

    if (!isone(Z)) ZT*=Z; // 5x

    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
    W5=W2-W1;

    YT=W5/(ZZn)2;   
}

//
//   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", 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[50],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;
    Big D[2*HERD],s,distance[50],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)
        distinguished=1<<(bits(leaps/16));
        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];
        }

        middle=(upper+lower)/2;
        if (ordermod>1)
            middle+=(ordermod+order-middle%ordermod);

        for (i=0;i<HERD;i++) start[i]=middle+13*ordermod*i;
        for (i=0;i<HERD;i++) start[HERD+i]=13*ordermod*i;

        for (i=0;i<2*HERD;i++) 
        {
            K[i]=start[i]*P;         // on your marks
            D[i]=0;                  // distance counter
        }
        cout << "Releasing " << HERD << " Tame and " << HERD << " Wild Kangaroos\n";

        nt=0; nw=0; cw=0; ct=0;
        collision=FALSE;  abort=FALSE;
        forever
        { 
            for (jj=0;jj<HERD;jj++)
            {
                K[jj].get(txc);
                i=txc%m;  /* random function */

                if (txc%distinguished==0)
                {
                    if (nt>=STORE)
                    {
                        abort=TRUE;
                        break;
                    }
                    cout << "." << flush;
                    tame[nt]=K[jj];
                    tdist[nt]=D[jj];
                    tname[nt]=jj;
                    for (k=0;k<nw;k++)
                    {
                        if (wild[k]==tame[nt])
                        {
                            ct=nt; cw=k;
                            collision=TRUE;
                            break;
                        }
                    }
                    if (collision) break;
                    nt++;
                }
                D[jj]+=distance[i];
                TE[jj]=table[i];
            }
            if (collision || abort) break;
            for (jj=HERD;jj<2*HERD;jj++)
            {
                K[jj].get(wxc);
                j=wxc%m;

                if (wxc%distinguished==0)
                {
                    if (nw>=STORE)
                    {
                        abort=TRUE;
                        break;
                    }
                    cout << "." << flush;
                    wild[nw]=K[jj];
                    wdist[nw]=D[jj];
                    wname[nw]=jj;
                    for (k=0;k<nt;k++)
                    {
                        if (tame[k]==wild[nw]) 
                        {
                            ct=k; cw=nw;
                            collision=TRUE;
                            break;
                        }
                    }
                    if (collision) break;
                    nw++;
                }
                D[jj]+=distance[j];
                TE[jj]=table[j];
            }
            if (collision || abort) break;

            multi_add(2*HERD,TE,K);   // jump together - its faster
        }
        cout << endl;
        if (abort)
        {
            cout << "Failed - this should not happen! - trying again" << endl;
            continue;
        }
        nrp=start[tname[ct]]-start[wname[cw]]+tdist[ct]-wdist[cw];
        G=P;
        G*=nrp;
        if (G!=ZERO)
        {
            cout << "Sanity Check Failed. Please report to [e-mail]" << endl;
            exit(0);
        }

        if (prime(nrp))
        {
            cout << "NP= " << nrp << endl;
            cout << "NP is Prime!" << endl;
            break;
        }

// final checks....
        real_order=nrp; i=0; 
        forever
        {
            sp=get_mip()->PRIMES[i];
            if (sp==0) break;
            if (real_order%sp==0)
            {
                G=P;
                G*=(real_order/sp);
                if (G==ZERO) 
                {
                    real_order/=sp;
                    continue;
                }
            }
            i++;
        }   
        if (real_order <= 4*sqrt(p))
        { 
            cout << "Low Order point used - trying again" << endl;
            continue;
        }
        real_order=nrp;
        for (i=0;(sp=get_mip()->PRIMES[i])!=0;i++)
            while (real_order%sp==0) real_order/=sp;
        if (real_order==1)
        { // all factors of nrp were considered
            cout << "NP= " << nrp << endl;
            break;
        }
        if (prime(real_order))
        { // all factors of NP except for one last one....
            G=P;
            G*=(nrp/real_order);
            if (G==ZERO) 
            {
                cout << "Failed - trying again" << endl;
                continue;
            }
            else 
            {
                cout << "NP= " << nrp << endl;
                break;
            }
        }
// Couldn't be bothered factoring nrp completely