schoof.cpp

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

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

    P[4].addterm((ZZn)(-4)*(8*B*B+A*A*A),0);
    P[4].addterm((ZZn)(-16)*(A*B),1);
    P[4].addterm((ZZn)(-20)*(A*A),2);
    P[4].addterm((ZZn)80*B,3);
    P[4].addterm((ZZn)20*A,4);
    P[4].addterm((ZZn)4,6);

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

    lower=5;     // next one to be calculated

  // Finding the order modulo 2
  // If GCD(X^P-X,X^3+AX+B) == 1 , trace=1 mod 2, else trace=0 mod 2

    XX=0;
    XX.addterm((ZZn)1,1);

    setmod(Y2);
    XP=pow(XX,p);
    G=gcd(XP-XX);         
    t[0]=0;
    if (isone(G)) t[0]=1;
    cout << "NP mod 2 =   " << (p+1-(int)t[0])%2;
    if ((p+1-(int)t[0])%2==0)
    {                                  
        cout << " ***" << endl;
        if (search) {b+=1; continue; }
    }
    else cout << endl;

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

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

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

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

        for (j=lower;j<=lp+1;j++)
        { // different for even and odd 
            if (j%2==1)
            { 
                n=(j-1)/2;
                if (n%2==0)
                    P[j]=P[n+2]*P3[n]*Y4-P3[n+1]*P[n-1];
                else
                    P[j]=P[n+2]*P3[n]-Y4*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])/(ZZn)2;
            }
            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;

        setmod(P[lp]);
        MY2=Y2;
        MY4=Y4;
// These next are time-consuming calculations of X^P, Y^P, X^(P*P) and Y^(P*P)

        cout << "X^P " << flush;
        XP=pow(XX,p);

// 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
        eigen=FALSE;
        if (!anomalous && prime((Big)lp))
        {
            PolyMod Xcoord,batch;
            batch=1;
            cout << "\b\b\b\bGCD " << flush;
            for (tau=1;tau<=(lp-1)/2;tau++)
            {
                if (tau%2==0)
                    Xcoord=(XP-XX)*P2[tau]*MY2+(PolyMod)P[tau-1]*P[tau+1];
                else
                    Xcoord=(XP-XX)*P2[tau]+(PolyMod)P[tau-1]*P[tau+1]*MY2;
                batch*=Xcoord;     
            }
            Fl=gcd(batch);  // just one GCD!
            if (degree(Fl)==(lp-1)) eigen=TRUE;
        }

        if (eigen)
        {
            setmod(Fl);
            MY2=Y2;
            MY4=Y4;

//
// Only the Y-coordinate is calculated. No need for X^P !
//
            cout << "\b\b\b\bY^P" << flush;
            YP=pow(MY2,(p-1)/2);
            cout << "\b\b\b";

//
// Now looking for value of lambda which satisfies
// (X^P,Y^P) = lambda.(XX,YY). 
//
// Note that it appears to be sufficient to only compare the Y coordinates (!?)
//
            cout << "NP mod " << lp << " = " << flush;
            Pf[0]=0; P2f[0]=0; P3f[0]=0;
            Pf[1]=1; P2f[1]=1; P3f[1]=1;    
            low=2;
            for (lambda=1;lambda<=(lp-1)/2;lambda++)
            {
                int res=0;
                PolyMod Ry,Ty;
                tau=(lambda+invers(lambda,lp)*p)%lp;

                cout << setw(3) << (p+1-tau)%lp << flush; 

// Get Divisor Polynomials as needed - this time mod the new (small) modulus Fl

                for (jj=low;jj<=lambda+2;jj++)
                    Pf[jj]=(PolyMod)P[jj]; 
                if (lambda+3>low) low=lambda+3;

// compare Y-coordinates - 5 polynomial mod-muls required

                P2f[lambda+1]=Pf[lambda+1]*Pf[lambda+1];
                P3f[lambda]=P2f[lambda]*Pf[lambda];
                if (lambda%2==0)
                {
                    Ry=(Pf[lambda+2]*P2f[lambda-1]-Pf[lambda-2]*P2f[lambda+1])/4;
                    Ty=MY4*YP*P3f[lambda];
                }
                else
                {
                    if (lambda==1) Ry=(Pf[lambda+2]*P2f[lambda-1]+P2f[lambda+1])/4;
                    else           Ry=(Pf[lambda+2]*P2f[lambda-1]-Pf[lambda-2]*P2f[lambda+1])/4;
                    Ty=YP*P3f[lambda];
                }

                if (degree(gcd(Ty-Ry))!=0) res=1;
                if (degree(gcd(Ty+Ry))!=0) res=2;
                if (res!=0) 
                {  // has it doubled, or become point at infinity?
                    if (res==2)
                    { // it doubled - wrong sign
                        tau=(lp-tau)%lp;
                        cout << "\b\b\b";
                        cout << setw(3) << (p+1-tau)%lp << flush; 
                    }
                    t[i]=tau;
                    if ((p+1-tau)%lp==0)
                    {
                        cout << " ***" << endl;
                        if (search) escape=TRUE;
                    }
                    else cout << endl;
                    break;
                }
                cout << "\b\b\b";
            }
            for (jj=0;jj<low;jj++)
            {  
                Pf[jj].clear();
                P2f[jj].clear();
                P3f[jj].clear();
            }
            if (escape) break;
            continue;
        }

 // no eigenvalue found, but some tau values can be eliminated...

        if (!anomalous && prime((Big)lp))
        {
            if (degree(Fl)==0) 
            {
                for (tau=0;tau<=lp/2;tau++)
                {
                    jj=(lp+tau*tau-(4*p)%lp)%lp;
                    if (jac(jj,lp)!=(-1)) permisso[tau]=FALSE;
                }
            }
            else
            { // Fl==P[lp] so tau=+/- sqrt(p) mod lp
                jj=(int)(2*sqrmp((p%lp),lp))%lp;
                for (tau=0;tau<=lp/2;tau++) permisso[tau]=FALSE;
                if (jj<=lp/2) permisso[jj]=TRUE;
                else          permisso[lp-jj]=TRUE;
            }
        }
        if (!prime((Big)lp))
        { // prime power
            for (jj=0;jj<start_prime;jj++)
                if (lp%(int)l[jj]==0)
                {
                    for (tau=0;tau<=lp/2;tau++)
                    {
                        permisso[tau]=FALSE;
                        if (tau%(int)l[jj]==(int)t[jj])      permisso[tau]=TRUE;
                        if ((lp-tau)%(int)l[jj]==(int)t[jj]) permisso[tau]=TRUE;
                    }
                    break;
                }  
        }

        cout << "\b\b\b\bY^P " << flush;
        YP=pow(MY2,(p-1)/2);
        cout << "\b\b\b\bX^PP" << flush;
       

        if (lp<40)
            XPP=compose(XP,XP);               // This is faster!
        else XPP=pow(XP,p);

        cout << "\b\b\b\bY^PP" << flush;
        if (lp<40)
            YPP=YP*compose(YP,XP);            // This is faster!
        else YPP=pow(YP,p+1);
        cout << "\b\b\b\b";

        PolyMod Pk,P2k,PkP1,PkM1,PkP2;
        Pk=P[k]; PkP1=P[k+1]; PkM1=P[k-1]; PkP2=P[k+2];

        P2k=(Pk*Pk);
//
// This is Schoof's algorithm, stripped to its bare essentials
//
// Now looking for the value of tau which satisfies 
// (X^PP,Y^PP) + k.(X,Y) =  tau.(X^P,Y^P)
// 
// Note that (X,Y) are rational polynomial expressions for points on
// an elliptic curve, so "+" means elliptic curve point addition
// 
// k.(X,Y) can be found directly from Divisor polynomials
// Schoof Prop (2.2)
//
// Points are converted to projective (X,Y,Z) form
// This is faster (x2). Observe that (X/Z^2,Y/Z^3,1) is the same
// point in projective co-ordinates as (X,Y,Z)
//

        if (k%2==0)
        {
            XT=XX*MY2*P2k-PkM1*PkP1;
            YT=(PkP2*PkM1*PkM1-P[k-2]*PkP1*PkP1)/4;
            XT*=MY2;      // fix up, so that Y has implicit y multiplier
            YT*=MY2;      // rather than Z
            ZT=MY2*Pk;
        }
        else
        {
            XT=(XX*P2k-MY2*PkM1*PkP1);
            if (k==1) YT=(PkP2*PkM1*PkM1+PkP1*PkP1)/4;
            else      YT=(PkP2*PkM1*PkM1-P[k-2]*PkP1*PkP1)/4;
            ZT=Pk;
        }

        elliptic_add(XT,YT,ZT,XPP,YPP,one);
// 
// Test for Schoof's case 1 - LHS (XT,YT,ZT) is point at infinity
//

        cout << "NP mod " << lp << " = " << flush;
        if (iszero(ZT))
        { // Is it zero point? (XPP,YPP) = - K(X,Y)
            t[i]=0;
            cout << setw(3) << (p+1)%lp;
            if ((p+1)%lp==0)
            {      
                cout << " ***" << endl;
                if (search) {escape=TRUE; break;}
            }
            else cout << endl;
            continue;         
        }

 // try all candidates one after the other

        PolyMod XP2,XP3,XP4,XP6,YP2,YP4;
        PolyMod ZT2XP,ZT2YP2,XPYP2,XTYP2,ZT3YP;

        ZT2=ZT*ZT;
        ZT3=ZT2*ZT;
        XP2=XP*XP;
        XP3=XP*XP2;
        XP4=XP2*XP2;
        XP6=XP3*XP3;
        YP2=MY2*(YP*YP);
        YP4=YP2*YP2;

        ZT2XP=ZT2*XP; ZT2YP2=ZT2*YP2; XPYP2=XP*YP2; XTYP2=XT*YP2; ZT3YP=ZT3*YP;

        Pf[0]=0; Pf[1]=1; Pf[2]=2;
        P2f[1]=1; P3f[1]=1;
        P2f[2]=Pf[2]*Pf[2];
        P3f[2]=P2f[2]*Pf[2];
    
        Pf[3]=3*XP4+6*A*XP2+12*B*XP-A*A;
        P2f[3]=Pf[3]*Pf[3];
        P3f[3]=P2f[3]*Pf[3];

        Pf[4]=(4*XP6+20*A*XP4+80*B*XP3-20*A*A*XP2-16*A*B*XP-32*B*B-4*A*A*A);
        P2f[4]=Pf[4]*Pf[4];
        P3f[4]=P2f[4]*Pf[4];
        low=5;

        for (tau=1;tau<=lp/2;tau++)
        {
            int res=0;
            PolyMod Rx,Tx,Ry,Ty;
            if (!permisso[tau]) continue;

            cout << setw(3) << (p+1-tau)%lp << flush;     

            for (jj=low;jj<=tau+2;jj++)
            { // different for odd and even
                if (jj%2==1)
                { /* 3 mod-muls */
                    n=(jj-1)/2;
                    if (n%2==0)
                       Pf[jj]=Pf[n+2]*P3f[n]*YP4-P3f[n+1]*Pf[n-1];
                    else
                       Pf[jj]=Pf[n+2]*P3f[n]-YP4*P3f[n+1]*Pf[n-1];
                }
                else
                { /* 3 mod-muls */
                    n=jj/2;
                    Pf[jj]=Pf[n]*(Pf[n+2]*P2f[n-1]-Pf[n-2]*P2f[n+1])/(ZZn)2;
                }
                P2f[jj]=Pf[jj]*Pf[jj];                            // square
                if (jj<=1+(1+(lp/2))/2) P3f[jj]=P2f[jj]*Pf[jj];   // cube
            } 
            if (tau+3>low) low=tau+3;

            if (tau%2==0)
            { // 4 mod-muls
                Rx=ZT2*(XPYP2*P2f[tau]-Pf[tau-1]*Pf[tau+1]);
                Tx=XTYP2*P2f[tau];
            }
            else
            { // 4 mod-muls
                Rx=(ZT2XP*P2f[tau]-ZT2YP2*Pf[tau-1]*Pf[tau+1]);
                Tx=XT*P2f[tau];
            }
            if (iszero(Rx-Tx))
            { // we have a result. Now compare Y's
                if (tau%2==0)
                {
                    Ry=ZT3YP*(Pf[tau+2]*P2f[tau-1]-Pf[tau-2]*P2f[tau+1]);
                    Ty=4*YT*YP4*P2f[tau]*Pf[tau];
                }
                else
                {
                    if (tau==1) Ry=ZT3YP*(Pf[tau+2]*P2f[tau-1]+P2f[tau+1]);
                    else        Ry=ZT3YP*(Pf[tau+2]*P2f[tau-1]-Pf[tau-2]*P2f[tau+1]);
                    Ty=4*YT*P2f[tau]*Pf[tau];
                }

                if (iszero(Ry-Ty)) res=1;
                else               res=2;
            }

            if (res!=0) 
            {  // has it doubled, or become point at infinity?
                if (res==2)
                { // it doubled - wrong sign
                    tau=lp-tau;
                    cout << "\b\b\b";
                    cout << setw(3) << (p+1-tau)%lp << flush; 
                }
                t[i]=tau;
                if ((p+1-tau)%lp==0)
                {
                    cout << " ***" << endl;
                    if (search) escape=TRUE;
                }
                else cout << endl;
                break;
            }
            cout << "\b\b\b";
        }
        for (jj=0;jj<low;jj++)
        {
            Pf[jj].clear();
            P2f[jj].clear();
            P3f[jj].clear();
        }
        if (escape) break;
    }
    Modulus.clear();

    for (i=0;i<=L+1;i++) 
    {
        P[i].clear(); // reclaim space
        P2[i].clear();
        P3[i].clear();
    }

    if (escape) {b+=1; continue;}
    Big order,ordermod;
    ordermod=1; for (i=0;i<nl-start_prime;i++) ordermod*=(int)l[start_prime+i];
    order=(p+1-CRT.eval(&t[start_prime]))%ordermod;    // get order mod product of primes

    nrp=kangaroo(p,order,ordermod);

    if (!prime(nrp) && search) {b+=1; continue; }
    else break;
    }

    if (fout) 
    {
        ECn P;
        ofile << bits(p) << endl;
        mip->IOBASE=16;
        ofile << p << endl;

        ofile << a << endl;
        ofile << b << endl;
    // generate a random point on the curve 
    // point will be of prime order for "ideal" curve, otherwise any point
        do {
            x=rand(p);
        } while (!P.set(x,x));
        P.get(x,y);
        ofile << nrp << endl;
        ofile << x << endl;
        ofile << y << endl;
        mip->IOBASE=10;
    }
    if (p==nrp) 
    {
        cout << "WARNING: Curve is anomalous" << endl;
        return 0;
    }

    if (p+1==nrp)
    {
        cout << "WARNING: Curve is supersingular" << endl;
    }

// check MOV condition for curves of Cryptographic interest
// if (pbits<128) return 0;

    d=1;
    for (i=1;i<50;i++)
    {
        d=modmult(d,p,nrp);
        if (d==1) 
        {
           if (i==1 || prime(nrp)) cout << "WARNING: Curve fails MOV condition - K = " << i << endl;
           else                    cout << "WARNING: Curve fails MOV condition - K <= " << i << endl;

           return 0;
        }
    }
    
    return 0;
}