cm.cpp

大数运算库

        a=A0;
        b=B0;
        at=-3;
        if (D==3) at=1;
        forever
        {
            if (D==1)
            {
                if (at<100)
                    eps=(ZZn)at/(ZZn)A0;
                else eps=rand(p);
                a=modmult(A0,eps,p);
            }
            if (D==3) 
            {
                if (at<100)
                    eps=(ZZn)at/(ZZn)B0;
                else eps=rand(p);
                b=modmult(B0,eps,p);
            }
            if (D!=1 && D!=3)
            {
                if (at<100)
                { // transform a to be simple if possible
                    w=(ZZn)at/ZZn(A0);
                    if (jacobi(w,p)!=1) 
                    {   
                        if (at==-3) at=1;
                        else at++; 
                        continue;
                    }
                    eps=sqrt(w,p);
                } else eps=rand(p);
                a=modmult(A0,pow(eps,2,p),p);
                b=modmult(B0,pow(eps,3,p),p);   
            } 
            ecurve(a,b,p,MR_PROJECTIVE);
            for (int xc=1;;xc++)
            {
                if (!pt.set((Big)xc)) continue;
                pt*=k;
                if (pt.iszero()) continue;
                break;
            }
            G=pt;                  // check its not the other one...

            if (r!=ord || !(other_r*G).iszero())  
            { 
                pt*=ord;
                if (pt.iszero()) break;
            }

            if (at==-3) at=1;
            else at++;
        }
        if (a==(p-3)) a=-3;
// check MOV condition
// A.12.1

        BOOL MOV=TRUE;
        
        w=1;
        for (i=1;i<100;i++)
        {
            w=modmult(w,p,ord);
            if (w==1) 
            {
               MOV=FALSE;
               if (!suppress) 
               {
                   if (i==1 || pord) cout << "\n*** Failed MOV condition - K = " << i << endl;
                   else                    cout << "\n*** Failed MOV condition - K <= " << i << endl;
               }
               break;
            }
        }
    
        if (!suppress)
        {
            if (MOV)  cout << "MOV condition OK" << endl;
            if (pord) cout << "\nCurve and Point Found" << endl;   
            else      cout << "\nCurve Found" << endl; 
        } 

        cout << "A= " << a << endl;
        cout << "B= " << b << endl;
        G.get(x,y);
        cout << "P= " << p << endl;
        cout << "R= " << ord;
        if (pord) 
        {
            cout << " a " << bits(ord) << " bit prime"  << endl;
            cout << "X= " << x << endl;
            cout << "Y= " << y << endl;
        }
        else            cout << " NOT prime" << endl;
        cout << endl;

        if (D!=1 && D!=3 )
        {
//            polly=polly/g;
//            if (degree(polly)>0)
//            {
                cout << "Try for different random factorisation (Y/N)? ";
                cin >> c;
                if (c=='Y' || c=='y') continue;
//            }
        }
        break;
    }
    if (pord) cout << "\nCurve and Point OK (Y/N)? " ;  
    else      cout << "\nCurve OK (Y/N)? " ;
    cin >> c;
    if (c=='N' || c=='n') 
    {
        if (!suppress) cout << "finding a curve..." << endl;
        return FALSE;
    }
    if (fout)
    {
        ofile << bits(p) << endl;
        mip->IOBASE=16;
        ofile << p << endl;
        ofile << a << endl;
        ofile << b << endl;
        ofile << ord << endl;
        if (pord)
        {
            ofile << x << endl;
            ofile << y << endl;
        }
        mip->IOBASE=10;
    }
    exit(0);
    return TRUE;
}

// Code to parse formula
// 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 '#'

void eval_power (Big& oldn,Big& n,char op)
{
        int i;
        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 (Big& t)
{
        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 (*s==' ')
        s++;
        if (*s=='-')    /* Unary minus */
        {
        s++;
        minus=1;
        }
        else
        minus=0;
        while (*s==' ')
        s++;
        if (*s=='(' || *s=='[' || *s=='{')    /* Number is subexpression */
        {
        s++;
        eval(t);
        n=t;
        }
        else            /* Number is decimal value */
        {
        for (i=0;s[i]>='0' && s[i]<='9';i++)
                ;
        if (!i)         /* No digits found */
        {
                cout <<  "Error - invalid number" << endl;
                exit (20);
        }
        op=s[i];
        s[i]=0;
        n=atol(s);
        s+=i;
        *s=op;
        }
        if (minus) n=-n;
        do
        op=*s++;
        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]);
        t=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;
}

int main(int argc,char **argv)
{
    BOOL good,dir;
    int i,ip,k,D,SD,precision;
    ofstream ofile;

    mip=mirsys(64,0);
    Big p,t;

    argv++; argc--;   
    irand(0L);
    if (argc<1)
    {
      cout << "Incorrect Usage" << endl;
      cout << "Program tries to find Elliptic Curve mod a prime P and point of prime order" << endl;
      cout << "that is a point (X,Y) on the curve y^2=x^3+Ax+B of order R" << endl;
      cout << "where R is large and prime > |P/K|. (K defaults to 256)" << endl;
      cout << "cm <prime number P>" << endl;
      cout << "OR" << endl;
      cout << "cm -f <formula for P>" << endl;
      cout << "e.g. cm -f 2#192-2#64-1" << endl;
      cout << "To suppress the commentary, use flag -s" << endl;
      cout << "(the commentary will make some sense to readers of IEEE 1363 Annex)" << endl;
      cout << "To search downwards for a prime, use flag -d" << endl;
      cout << "To output to a file, use flag -o <filename>" << endl;
      cout << "To set maximum  co-factor size K, use e.g. flag -k8" << endl;
      cout << "To set infinite co-factor size K, use flag -k0" << endl;
      cout << "(In this case the reported R is the number of points on the curve)" << endl;
      cout << "To start searching from a particular D value, use e.g. flag -D10000" << endl;
      cout << "To set MIRACL precision, use e.g. flag -P64 (default -P32)" << endl;
      cout << "If program fails, try again with double precision" << endl;  
      cout << "e.g. cm -f 2#224-2#96+1 -k12 -D400000 -P64 -o common.ecs" << endl;
      cout << "\nFreeware from Shamus Software, Dublin, Ireland" << endl;
      cout << "Full C++ source code and MIRACL multiprecision library available" << endl;
      cout << "Email to [e-mail] for details" << endl;
      return 0;
    }

    gprime(1000);
    precision=32;
    ip=0;
    k=256;
    SD=1;
    suppress=FALSE;
    fout=FALSE;
    dir=FALSE;
    p=0;
    while (ip<argc)
    {
        if (strcmp(argv[ip],"-f")==0)
        {
            if (p==0)
            {
                ip++;
                s=argv[ip++];
                t=0;
                eval(t);
                p=t;
                continue;
            }
            else
            {
                cout << "Error in command line" << endl;
                return 0;
            }
        }    
        if (strcmp(argv[ip],"-s")==0)
        {
            ip++;
            suppress=TRUE;
            continue;
        }
        if (strcmp(argv[ip],"-d")==0)
        {
            ip++;
            dir=TRUE;
            continue;
        }
        if (strncmp(argv[ip],"-k",2)==0)
        {
             k=atoi(argv[ip]+2);
             ip++;
             continue;
        }
        if (strncmp(argv[ip],"-D",2)==0)
        {
             SD=atoi(argv[ip]+2);
             ip++;
             continue;
        }
        if (strncmp(argv[ip],"-P",2)==0)
        {
             precision=atoi(argv[ip]+2);
             ip++;
             continue;
        }
        if (strcmp(argv[ip],"-o")==0)
        {
            ip++;
            fout=TRUE;
            ofile.open(argv[ip++]);
            continue;
        }
        if (p==0) p=argv[ip++];
        else
        {
            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;
    }
    if (p<10)
    {
        cout << "Prime is too small - use another method!" << endl;
        return 0;
    }
    if (bits(p)>1024)
    {
        cout << "Prime is too big - sorry" << endl;
        return 0;
    }

    mirexit();

    if (precision<16) precision=16;

    mip=mirsys(precision,0);     // restart with new precision
    gprime(100000);

    Big W,V,K,r1,r2,ord,rmin;
    Complex lam;

    Flash Fi2[7];
    Flash pi24;
    
    if (k==0) rmin=1;
    else      rmin=(p-2*sqrt(p))/k;
    if (rmin==0) 
    {
        cout << "Bad k co-factor value" << endl;
        return 0;
    }

    W=sqrt(p)+1;
    K=(W*W)/rmin;

    if (!suppress) cout << "P mod 8 = " << p%8 << endl;
    if (!suppress) cout << "P is " << bits(p) << " bits long" << endl;
    if (!suppress) cout << "precomputations..." << endl;

    pi24=pi()/(Flash)24; 
    lam=exp(Complex((Flash)0,pi24));

    Fi2[0]=1;
    Fi2[2]=pow((Flash)2,Flash(-1,3));
    Fi2[3]=pow((Flash)2,Flash(-1,2));
    Fi2[6]=Flash(1,2);

    if (!suppress) cout << "finding a curve..." << endl;

    for (D=SD;;D++)
    {
        if (!isaD(D,p%8,K)) continue;
        if (jacobi(p-D,p)==-1) continue;
        good=TRUE;
// A.14.2.3
        for (i=1;;i++)
        {
            int sp=mip->PRIMES[i];
            if (D==sp || sp*sp>D) break;
            if (D%sp==0 && jacobi(p,(Big)sp)==-1)
            {
                good=FALSE;
                break;
            }
        }
        if (!good) continue;
        if (!isacm(p,D,W,V)) continue;

        r1=p+1+W;
        r2=p+1-W;

        if (k==0) ord=r1;
        else      ord=trial_divide(r1);
        if (ord==1) ord=r1;
        if (ord>rmin && (k==0 || prime(ord))) get_curve(lam,Fi2,ofile,r1,r2,ord,D,p,W);

        if (k==0) ord=r2;
        else      ord=trial_divide(r2);
        if (ord==1) ord=r2;
        if (ord>rmin && (k==0 || prime(ord))) get_curve(lam,Fi2,ofile,r2,r1,ord,D,p,W);

        if (D==1)
        {
            r1=p+1+V;
            r2=p+1-V;
            if (k==0) ord=r1;
            else      ord=trial_divide(r1);
            if (ord==1) ord=r1;
            if (ord>rmin && (k==0 || prime(ord))) get_curve(lam,Fi2,ofile,r1,r2,ord,D,p,W);

            if (k==0) ord=2;
            else      ord=trial_divide(r2);
            if (ord==1) ord=r2;
            if (ord>rmin && (k==0 || prime(ord))) get_curve(lam,Fi2,ofile,r2,r1,ord,D,p,W);
            
        }
        if (D==3)
        {
            r1=p+1+(W+3*V)/2;
            r2=p+1-(W+3*V)/2;
            if (k==0) ord=r1;
            else ord=trial_divide(r1);
            if (ord==1) ord=r1;
            if (ord>rmin && (k==0 || prime(ord))) get_curve(lam,Fi2,ofile,r1,r2,ord,D,p,W);

            if (k==0) ord=r2;
            else ord=trial_divide(r2);
            if (ord==1) ord=r2;
            if (ord>rmin && (k==0 || prime(ord))) get_curve(lam,Fi2,ofile,r2,r1,ord,D,p,W);

            r1=p+1+(W-3*V)/2;
            r2=p+1-(W-3*V)/2;
            if (k==0) ord=r1;
            else ord=trial_divide(r1);
            if (ord==1) ord=r1;
            if (ord>rmin && (k==0 || prime(ord))) get_curve(lam,Fi2,ofile,r1,r2,ord,D,p,W);

            if (k==0) ord=r2;
            else ord=trial_divide(r2);
            if (ord==1) ord=r2;
            if (ord>rmin && (k==0 || prime(ord))) get_curve(lam,Fi2,ofile,r2,r1,ord,D,p,W);
        }
    }
    cout << "No satisfactory curve found" << endl;
    return 0;
}