schoof2.cpp

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

            }
        }
        if (bad)
        {
            cout << "Failed - trying again" << endl;
            continue;                       
        }
        cout << "NP is composite and not ideal for Cryptographic use" << endl;
        cout << "NP= " << nrp << " (probably)" << endl;
        break;
    }
    return nrp;
}

int qpow(int x,int y)
{ // quick and dirty power function
    int r=x;
    for (int i=1;i<y;i++) r*=x;
    return r;
}

int main(int argc,char **argv)
{
    ofstream ofile;
    int TR,M,a,b,c,low,lower,ip,lp,i,j,jj,m,n,nl,L,k,tau,lambda,cf;
    mr_utype t[100];
    Big aa,bb,p,nrp,x,y,d,s;
    Poly2Mod X2,XP,YP,YPy,XP2,XPP,YPP,YPPy,TT,SX,XK[600];
    Poly2Mod Pf[600],P2f[600],P3f[600];
    Poly2 Fl,X,G,P[100],P2[100],P3[100],FX,Y4;
    EC2 GG;
    miracl *mip=&precision;
    BOOL eigen,found,escape,search,fout,gotM,gotA,gotB,gota,gotb,gotc;
    static BOOL permisso[600];
    int Base;

    argv++; argc--;
    if (argc<1)
    {
        cout << "Incorrect Usage" << endl;
        cout << "Program finds the number of points (NP) on an Elliptic curve" << endl;
        cout << "which is defined over the Galois field GF(2^M)" << endl;
        cout << "The Elliptic Curve has the equation Y^2 +XY = X^3 + AX^2 + B " << endl;
        cout << "A suitable trinomial or Pentanomial basis must" << endl;
        cout << "also be specified t^m+t^a+1 or t^m+t^a+t^b+t^c+1" << endl;
        cout << "These can be found in e.g. the IEEE P1363 standard" << endl;
        cout << "or can be generated using the MIRACL findbase example program" << endl;
        cout << "schoof2 <A> <B> <M> <a> <b> <c>" << endl << endl;
        cout << "e.g. schoof2 1 52 191 9" << endl << endl;
        cout << "To input A and B in Hex, precede with -h" << endl;
        cout << "To output to a file, use flag -o <filename>" << endl;
        cout << "To search for NP a near-prime, incrementing B, use flag -s" << endl;
        cout << "\nFreeware from Shamus Software, Dublin, Ireland" << endl;
        cout << "Full C++ source code and MIRACL multiprecision library available" << endl;
        cout << "http://indigo.ie/~mscott for details" << endl;
        cout << "or email [e-mail]" << endl;
        return 0;
    }

    ip=0;
    gprime(10000);                 // generate small primes < 1000
    search=fout=gotM=gotA=gotB=gota=gotb=gotc=FALSE;
    M=a=b=c=0;

// Interpret command line
    Base=10;
    while (ip<argc)
    {
        if (strcmp(argv[ip],"-o")==0)
        {
            ip++;
            if (ip<argc)
            {
                fout=TRUE;
                ofile.open(argv[ip++]);
                continue;
            }
            else
            {
                cout << "Error in command line" << endl;
                return 0;
            }
        }
        if (strcmp(argv[ip],"-s")==0)
        {
            ip++;
            search=TRUE;
            continue;
        }
        if (strcmp(argv[ip],"-h")==0)
        {
            ip++;
            Base=16;
            continue;
        }
        if (!gotA) 
        {
            mip->IOBASE=Base;
            aa=argv[ip++]; 
            mip->IOBASE=10;
            gotA=TRUE;
            continue;
        }
        if (!gotB) 
        {
            mip->IOBASE=Base;
            bb=argv[ip++];
            mip->IOBASE=10;
            gotB=TRUE;
            continue;
        }
        if (!gotM) 
        {
            M=atoi(argv[ip++]);
            gotM=TRUE;
            continue;
        }
        if (!gota)
        {
            a=atoi(argv[ip++]);
            gota=TRUE;
            continue;
        } 
        if (!gotb)
        {
            b=atoi(argv[ip++]);
            gotb=TRUE;
            continue;
        } 
        if (!gotc)
        {
            c=atoi(argv[ip++]);
            gotc=TRUE;
            continue;
        } 
        cout << "Error in command line" << endl;
        return 0;
    }    

    if ((!gotM || !gotA || !gotB) || a==0)
    {
        cout << "Error in command line" << endl;
        return 0;
    }

// loop for "-s" search option
    p=pow((Big)2,M);
    forever {
 
    A=aa;
    if (bb>=p) bb%=p;
  
    B=bb;
    if (!ecurve2(M,a,b,c,A,B,TRUE,MR_AFFINE))    // initialise Elliptic Curve
    {
        cout << "Illegal Curve Parameters" << endl;
        return 0;
    }

    D=mip->C;
    GF2m At=A;
    TR=trace(At);

// The elliptic curve as a Polynomial

    FX=0;
    FX.addterm(B,0);
    FX.addterm((GF2m)A,2);

    FX.addterm((GF2m)1,3);

    X=0; X.addterm(1,1);

    cout << "Counting the number of points (NP) on the curve" << endl;
    cout << "y^2 + xy = " << FX << " mod 2^" << M << endl;   

    if (B==0)
    {
        cout << "Not Allowed! B = 0" << endl;
        if (search) {bb+=1; continue; }
        else return 0;
    }

    if (M<12)
    { // do it the simple way
        nrp=2;
        x=1;
        while (x< (1<<M))
        { /* check if two points exist for each x */
            if (GG.set(x,0)) nrp+=2;
            x+=1;
        }
 
        do 
        {
            x=rand(p);
        } while (!GG.set(x,x));

        GG*=nrp;
        if (!GG.iszero())
        {
            cout << "Sanity Check Failed. Please report to [e-mail]" << endl;
            exit(0);
        }

        if (prime(nrp/2))
        {
            cout << "NP is 2*Prime!" << endl;
            cout << "NP= 2*" << nrp/2 << endl;
            break;
        }
        if (nrp%4==0)
        {
            if (A==0 && prime(nrp/4))
            {
                cout << "NP is 4*Prime!" << endl;
                cout << "NP= 4*" << nrp/4 << endl;
                break;
            }
        }
        cout << "NP= " << nrp << endl;
        if (search) {bb+=1; continue; }

        break;
        
    } 
    else if (M<56) 
    {
        nrp=kangaroo(p,(Big)0,(Big)2,TR,found);
        if (found) break;
        if (search) {bb+=1; continue; }
        break;
    }

    if (M<=100) d=pow((Big)2,48);  
    if (M>100 && M<=120) d=pow((Big)2,56); 
    if (M>120 && M<=160) d=pow((Big)2,64); 
    if (M>160 && M<=200) d=pow((Big)2,72);
    if (M>200) d=pow((Big)2,80);
    d=sqrt(p/d);
    if (d<256) d=256;

    mr_utype l[100]; 
    int pp[100];   // primes and powers 

// see how many primes will be needed
// l[.] is the prime, pp[.] is the power

    for (s=1,nl=0;s<=d;nl++)
    {
        int tp=mip->PRIMES[nl];
        pp[nl]=1;     // every prime included once...
        s*=tp;
        for (i=0;i<nl;i++)
        { // if a new prime power is now in range, include its contribution
            int cp=mip->PRIMES[i];
            int p=qpow(cp,pp[i]+1);
            if (p<tp  || (cp==2 && p<16*tp) )
            { // new largest prime power
                s*=cp;
                pp[i]++;
            } 
        }
    }  
    L=mip->PRIMES[nl-1];

    cout << nl << " primes used (plus largest prime powers), largest is " << L << endl;

    for (i=0;i<nl;i++)
        l[i]=mip->PRIMES[i];

    int start_prime;   // start of primes & largest prime powers
    for (i=0;;i++)
    {
        if (pp[i]!=1)
        {
            int p=qpow(l[i],pp[i]);
            for (j=0;l[j]<p && j<nl;j++) ;
            nl++;
            for (m=nl-1;m>j;m--) l[m]=l[m-1];
            l[j]=p;   // insert largest prime power in table
        }
        else 
        {
            start_prime=i;
            break;
        } 
    }

// table of primes and prime powers now looks like:-
// 2 3 5 7 9 11 13 16 17 19 ....
//     S   p        p

// where S is start_prime, and p marks the largest prime powers in the range
// CRT uses primes starting from S, but small primes are kept in anyway, 
// as they allow quick abort if searching for prime NP.

// Calculate Divisor Polynomials 
// Set the first few by hand....

    P[0]=0; P[1]=1; P[2]=0; P[3]=0; P[4]=0;

    P2[1]=1; P3[1]=1;

    P[2].addterm(1,1);

    P2[2]=P[2]*P[2];
    P3[2]=P2[2]*P[2];

    P[3].addterm(B,0);   
    P[3].addterm(1,3);
    P[3].addterm(1,4);  

    P2[3]=P[3]*P[3];
    P3[3]=P2[3]*P[3];

    P[4].addterm(B,2);
    P[4].addterm(1,6);

    P2[4]=P[4]*P[4];
    P3[4]=P2[4]*P[4];

    lower=5;     // next one to be calculated

    t[0]=0;

    Poly2Mod zero,one,XT,XTy,YT,YTy,ZT,XL,YL,ZL,ZL2,ZT2;
    one=1;                  // polynomial = 1
    zero=0;                 // polynomial = 0
    
    Crt CRT(nl-start_prime,&l[start_prime]);  // initialise for application of 
                                              // chinese remainder thereom

// now look for trace%prime for prime=3,5,7,11 etc
// actual trace is found by combining these via CRT

    l[0]=4;
    if (TR==0) l[0]=8;

    escape=FALSE;
    for (i=0;i<nl;i++)      
    {
        lp=l[i];       // next prime
        k=p%lp;

// generation of Divisor polynomials as needed
// See Schoof p. 485

        if (lp<=8 || lp%2!=0)
        {
            for (j=lower;j<=lp+1;j++)
            { // different for even and odd 
                if (j%2==1)
                { 
                    n=(j-1)/2;
                    P[j]=P[n+2]*P3[n]+P3[n+1]*P[n-1];
                }
                else
                {
                    n=j/2;
                    P[j]=P[n]*(P[n+2]*P2[n-1]+P[n-2]*P2[n+1]);
                    P[j]=divxn(P[j],1);
                }
                if (j <= 1+(L+1)/2)
                { // precalculate for later
                    P2[j]=P[j]*P[j];
                    P3[j]=P2[j]*P[j];
                }
            }

            if (lp+2>lower) lower=lp+2;
        }    
        for (tau=0;tau<=lp/2;tau++) permisso[tau]=TRUE;

        eigen=FALSE;
        if (lp%2==0)
        { // its 2^c
            Poly2 G,G2,PI;
            GF2m S;
            int t=lp,c=0;
            while (t!=1) {t/=2; c++;}

//
// When lp is a power of 2, we can construct a root of the Division Polynomial
// directly. See Menezes P.107. This is much quicker. Then the eigenvalue
// heuristic can be used.
//

            S=sqrt(sqrt((GF2m)B));
            G=X+S;
            G2=G*G;
            PI=1;

            for (jj=2;jj<c;jj++)
            {
                S=sqrt(S);
                G=G2+S*X*PI;
                PI=PI*G2;
                G2=G*G;
            }  
            eigen=TRUE;
            Fl=G;          // G is a root of degree lp/4
        }
        else
        {
            Poly2Mod Xcoord,batch;
            setmod(P[lp]);
            MFX=FX;
            XX=X;
            XP2=1;
  
//
// Calculate X^P. Save intermediates, as they can be used to
// calculate Y^P (if we need to.. )
//

            cout << "X^P " << flush;
            for (jj=0;jj<M-1;jj++)
            {
                XP2*=XX;
                XP2*=XP2;
                XK[jj]=XP2;
            }
            XP=XP2*(XX*XX);

// Eigenvalue search - see Menezes
// Batch the GCDs as they are slow. 
// This gives us product of both eigenvalues - a polynomial of degree (lp-1)
// But thats a lot better than (lp^2-1)/2

            if (prime((Big)lp))
            {
                batch=1;
                cout << "\b\b\b\bGCD " << flush;  
                for (tau=1;tau<=(lp-1)/2;tau++)
                {
                    Xcoord=(XP+XX)*P2[tau]+(Poly2Mod)P[tau-1]*P[tau+1];
                    batch*=Xcoord;
                }
                Fl=gcd(batch);
                if (degree(Fl)==(lp-1)) eigen=TRUE;
            }
            cout << "\b\b\b\b";
        }   

        if (eigen)
        { // eigenvalue found!
            Poly2Mod one;
            setmod(Fl);
            one=1;
            MFX=FX;