schoof.cpp

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

// Probably not an interesting curve for Cryptographic purposes anyway.....
// But if 20 random points are all "killed" by nrp, its almost
// certain to be the true NP, and not a multiple of a small order.

        bad=FALSE;
        for (i=0;i<20;i++)
        { 
            do
            {
                x=rand(p);
            } while (!P.set(x,x));
            G=P;
            G*=nrp;
            if (G!=ZERO)
            {
                bad=TRUE;
                break;
            }
        }
        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;
}

// Code to parse formula in command line
// This code isn't mine, but its public domain
// Shamefully I forget the source
//
// NOTE: It may be necessary on some platforms to change the operators * and #

#define TIMES '*'
#define RAISE '#'

Big tt;
static char *ss;

void eval_power (Big& oldn,Big& n,char op)
{
        if (op) n=pow(oldn,toint(n));    // power(oldn,size(n),n,n);
}

void eval_product (Big& oldn,Big& n,char op)
{
        switch (op)
        {
        case TIMES:
                n*=oldn; 
                break;
        case '/':
                n=oldn/n;
                break;
        case '%':
                n=oldn%n;
        }
}

void eval_sum (Big& oldn,Big& n,char op)
{
        switch (op)
        {
        case '+':
                n+=oldn;
                break;
        case '-':
                n=oldn-n;
        }
}

void eval (void)
{
        Big oldn[3];
        Big n;
        int i;
        char oldop[3];
        char op;
        char minus;
        for (i=0;i<3;i++)
        {
            oldop[i]=0;
        }
LOOP:
        while (*ss==' ')
        ss++;
        if (*ss=='-')    /* Unary minus */
        {
        ss++;
        minus=1;
        }
        else
        minus=0;
        while (*ss==' ')
        ss++;
        if (*ss=='(' || *ss=='[' || *ss=='{')    /* Number is subexpression */
        {
        ss++;
        eval ();
        n=tt;
        }
        else            /* Number is decimal value */
        {
        for (i=0;ss[i]>='0' && ss[i]<='9';i++)
                ;
        if (!i)         /* No digits found */
        {
                cout <<  "Error - invalid number" << endl;
                exit (20);
        }
        op=ss[i];
        ss[i]=0;
        n=atoi(ss);
        ss+=i;
        *ss=op;
        }
        if (minus) n=-n;
        do
        op=*ss++;
        while (op==' ');
        if (op==0 || op==')' || op==']' || op=='}')
        {
        eval_power (oldn[2],n,oldop[2]);
        eval_product (oldn[1],n,oldop[1]);
        eval_sum (oldn[0],n,oldop[0]);
        tt=n;
        return;
        }
        else
        {
        if (op==RAISE)
        {
                eval_power (oldn[2],n,oldop[2]);
                oldn[2]=n;
                oldop[2]=RAISE;
        }
        else
        {
                if (op==TIMES || op=='/' || op=='%')
                {
                eval_power (oldn[2],n,oldop[2]);
                oldop[2]=0;
                eval_product (oldn[1],n,oldop[1]);
                oldn[1]=n;
                oldop[1]=op;
                }
                else
                {
                if (op=='+' || op=='-')
                {
                        eval_power (oldn[2],n,oldop[2]);
                        oldop[2]=0;
                        eval_product (oldn[1],n,oldop[1]);
                        oldop[1]=0;
                        eval_sum (oldn[0],n,oldop[0]);
                        oldn[0]=n;
                        oldop[0]=op;
                }
                else    /* Error - invalid operator */
                {
                        cout <<  "Error - invalid operator" << endl;
                        exit (20);
                }
                }
        }
        }
        goto LOOP;
}

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

int main(int argc,char **argv)
{
    ofstream ofile;
    int low,lower,ip,pbits,lp,i,j,jj,m,n,nl,L,k,tau,lambda;
    mr_utype t[100];
    Big a,b,p,nrp,x,y,d,s;
    PolyMod XX,XP,YP,XPP,YPP;
    PolyMod Pf[100],P2f[100],P3f[100];
    Poly G,P[100],P2[100],P3[100],Y2,Y4,Fl;
    miracl *mip=&precision;
    BOOL escape,search,fout,dir,gotP,gotA,gotB,eigen,anomalous;
    BOOL permisso[100];
    ZZn delta,j_invariant;
    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(P), P a prime" << endl;
        cout << "The Elliptic Curve has the equation Y^2 = X^3 + AX + B mod P" << endl;
        cout << "schoof <prime number P> <A> <B>" << endl;
        cout << "OR" << endl;
        cout << "schoof -f <formula for P> <A> <B>" << endl;
        cout << "e.g. schoof -f 2#192-2#64-1 -3 35317045537" << endl;
        cout << "To output to a file, use flag -o <filename>" << endl;
        cout << "To search downwards for a prime, use flag -d" << endl;
        cout << "To input P, A and B in Hex, precede with -h" << endl;
        cout << "To search for NP 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 << "Also faster Schoof-Elkies-Atkin implementation" << 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=dir=gotP=gotA=gotB=FALSE;
    p=0; a=0; b=0;

// Interpret command line
    Base=10;
    while (ip<argc)
    {
        if (strcmp(argv[ip],"-f")==0)
        {
            ip++;
            if (!gotP && ip<argc)
            {

                ss=argv[ip++];
                tt=0;
                eval();
                p=tt;
                gotP=TRUE;
                continue;
            }
            else
            {
                cout << "Error in command line" << endl;
                return 0;
            }
        }
        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],"-d")==0)
        {
            ip++;
            dir=TRUE;
            continue;
        }
        if (strcmp(argv[ip],"-s")==0)
        {
            ip++;
            search=TRUE;
            continue;
        }
        if (strcmp(argv[ip],"-h")==0)
        {
            ip++;
            Base=16;
            continue;
        }
     
        if (!gotP) 
        {
            mip->IOBASE=Base;
            p=argv[ip++];
            mip->IOBASE=10;
            gotP=TRUE;
            continue;
        }
        if (!gotA) 
        {
            mip->IOBASE=Base;
            a=argv[ip++];
            mip->IOBASE=10;
            gotA=TRUE;
            continue;
        }
        if (!gotB) 
        {
            mip->IOBASE=Base;
            b=argv[ip++];
            mip->IOBASE=10;
            gotB=TRUE;
            continue;
        }
        cout << "Error in command line" << endl;
        return 0;
    }    

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

    if (!prime(p))
    {
        int incr=0;
        cout << "That number is not prime!" << endl;
        if (dir)
        {
            cout << "Looking for next lower prime" << endl;
            p-=1; incr++;
            while (!prime(p)) { p-=1;  incr++; }
            cout << "Prime P = P-" << incr << endl;
        }
        else
        {
            cout << "Looking for next higher prime" << endl;
            p+=1; incr++;
            while (!prime(p)) { p+=1;  incr++; }
            cout << "Prime P = P+" << incr << endl;
        }
        cout << "Prime P = " << p << endl;
    }
    pbits=bits(p);
    cout << "P mod 24 = " << p%24 << endl;
    cout << "P is " << pbits << " bits long" << endl;

// loop for "-s" search option

    forever {

    fft_reset();   // reset FFT tables

    ecurve(a,b,p,MR_AFFINE);    // initialise Elliptic Curve

    A=a;
    B=b;

// The elliptic curve as a Polynomial

    Y2=0;
    Y2.addterm(B,0);
    Y2.addterm(A,1);
    Y2.addterm((ZZn)1,3);

    Y4=Y2*Y2;
    cout << "Counting the number of points (NP) on the curve" << endl;
    cout << "y^2= " << Y2 << " mod " << p << endl;   

    delta=-16*(4*A*A*A+27*B*B);
    if (delta==0)
    {
        cout << "Not Allowed! 4A^3+27B^2 = 0" << endl;
        if (search) {b+=1; continue; }
        else return 0;
    }
    anomalous=FALSE;
    j_invariant=(-1728*64*A*A*A)/delta;

    cout << "j-invariant= " << j_invariant << endl;

    if (j_invariant==0 || j_invariant==1728)
    {
        anomalous=TRUE;
        cout << "Warning: j-invariant is " << j_invariant << endl;
    }
    if (pbits<14)
    { // do it the simple way
        nrp=1;
        x=0;
        while (x<p)
        {
            nrp+=1+jacobi((x*x*x+a*x+b)%p,p);
            x+=1;
        }
        cout << "NP= " << nrp << endl;
        if (prime(nrp)) cout << "NP is Prime!" << endl;
        else if (search) {b+=1; continue; }
        break;
    } 
    if (pbits<56) 
    { // do it with kangaroos
        nrp=kangaroo(p,(Big)0,(Big)1);
        if (!prime(nrp) && search) {b+=1; continue; }
        break;
    }

    if (pbits<=100) d=pow((Big)2,48);
    if (pbits>100 && pbits<=120) d=pow((Big)2,56);
    if (pbits>120 && pbits<=140) d=pow((Big)2,64);
    if (pbits>140 && pbits<=200) d=pow((Big)2,72);
    if (pbits>200) d=pow((Big)2,80);
    
/*
    if (pbits<200) d=pow((Big)2,72);
    else           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)
            { // 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)
        {
            mr_utype p=qpow(l[i],pp[i]);
            for (j=0;l[j]<p;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 - Schoof 1985 p.485
// Set the first few by hand....
    
    P[1]=1; P[2]=2; P[3]=0; P[4]=0;

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

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

    P[3].addterm(-(A*A),0); P[3].addterm(12*B,1);
    P[3].addterm(6*A,2)   ; P[3].addterm((ZZn)3,4);

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