u[m]=compose(u[m-1],u[m-1]);
                        lim=m;
                    }
                }

                if (iszero(C-XX)) 
                { // found splitting type
                    useful=TRUE;
                    break;
                }
            }
            cout << "\b\b\b" << flush;
            if (!useful) continue;

            cout << "NP mod " << lp << " = " << flush;

            int a,b,candidates,gx,gy,ord,qnr=2;
            BOOL gen;
            while (jac(qnr,lp)!=(-1)) qnr++;

//
// [4] Algorithm VII.4 - find a generator of F(lp^2)
//
            ord=lp*lp-1;
            gy=1;
            for (gx=1;gx<lp;gx++)
            {
                gen=TRUE;
                for (jj=2;jj<=ord/2;jj++)
                {
                    if (ord%jj!=0) continue;
                    powquad(lp,qnr,gx,gy,ord/jj,a,b);
                    if (a==1 && b==0) {gen=FALSE; break;}
                }
                if (gen) break;
            }
//
// (gx,gy) is a generator
//
            candidates=0;
            cout << setw(3);
            for (jj=1;jj<r;jj++)
            {
                if (jj>1 && igcd(jj,r)!=1) continue;
                powquad(lp,qnr,gx,gy,jj*ord/r,a,b);

                tau=((a+1)*k*(int)invers(2,lp))%lp;
                if (tau==0)
                {           // r must be 2 - I can make use of this!
                            // Its an Atkin prime, but only one possibility
                    candidates++;
                    cout << (p+1)%lp << flush;
                    if ((p+1)%lp==0)
                    {
                        cout << " ***" << endl;
                        if (search) escape=TRUE;
                    }
                    else cout << endl; 
                    good[nl]=lp;
                    t[nl]=tau;
                    nl++;
                    accum*=lp;  
                    break;
                }
                else if (jac(tau,lp)==1)
                {
                    candidates+=2;
                    tau=sqrmp(tau,lp);
                    tau=(2*tau)%lp;
                    if (candidates==phi(r))
                    { 
                         cout << (p+1-tau)%lp << " or " << (p+1+tau)%lp << endl;   
                         break;
                    }
                    else cout << (p+1-tau)%lp << "," << (p+1+tau)%lp << "," << flush;
                }  
            }
            if (escape) break;
            continue;  
        }

//
// Good Elkies prime - so use it!
//
// First solve quadratic for a root
//
        if (degree(G)==1)
        {
            discrim=0; 
            g=-G.coeff(0); // Elkies Prime, one root, (2 possibilites)   
        }
        else               // degree(G)==2
        {                  // Elkies Prime, two roots
            discrim=1;
            qb=G.coeff(1);
            qc=G.coeff(0); 
            g=sqrt(qb*qb-4*qc);
            g=(-qb-g)/2;   // pick either one
        }
        cout << "\b\b\bELK" << flush;
//
// Mueller's procedure for finding the atilde, btilde and p1
// parameters of the isogenous curve
// 3. page 111
// 4. page 131-133
// First we need partial differentials of bivariate Modular Polynomial
//

        dGx=diff_dx(MP);
        dGy=diff_dy(MP);
        dGxx=diff_dx(dGx);
        dGxy=diff_dx(dGy);
        dGyy=diff_dy(dGy);

        Eg=dGx.F(g,j);   // Evaluated at (g,j)
        Ej=dGy.F(g,j);
        Exy=dGxy.F(g,j);

        Dg=g*Eg;    
        Dj=j*Ej;

        deltal=delta*pow(g,12/is)/pow(el,12);

        if (Dj==0)
        {
            E4bl=E4b/(el*el);
            atilde=-3*pow(el,4)*E4bl;
            jl=pow(E4bl,3)/deltal;
            btilde=2*pow(el,6)*sqrt((jl-1728)*deltal);
            p1=0;
        }
        else
        {
            E2bs=(-12*E6b*Dj)/(s*E4b*Dg);

            gd=-(s/12)*E2bs*g;
            jd=-E4b*E4b*E6b/delta;
            E0b=E6b/(E4b*E2bs); 

            Dgd=gd*Eg+g*(gd*dGxx.F(g,j)+jd*Exy);
            Djd=jd*Ej+j*(jd*dGyy.F(g,j)+gd*Exy);  
   
            E0bd=((-s*Dgd)/12-E0b*Djd)/Dj;

            E4bl=(E4b-E2bs*(12*E0bd/E0b+6*E4b*E4b/E6b-4*E6b/E4b)+E2bs*E2bs)/(el*el);

            jl=pow(E4bl,3)/deltal;
            f=pow(el,is)/g; fd=s*E2bs*f/12;

            Dgs=dGx.F(f,jl);
            Djs=dGy.F(f,jl);

            jld=-fd*Dgs/(el*Djs);
            E6bl=-E4bl*jld/jl;

            atilde=-3*pow(el,4)*E4bl;
            btilde=-2*pow(el,6)*E6bl;
            p1=-el*E2bs/2;            
        }

//
// Find factor of Division Polynomial from atilde, btilde and p1 
// Here we follow 3. p 116
// Polynomials have been modified s.t x=z^2
//
// Note that all Polynomials can be reduced mod x^(d+1),
// where d=(lp-1)/2, using modxn() function
//

        cout << "\b\b\bFAC" << flush;
        ld=(lp-1)/2;
        ld1=(lp-3)/2;

        get_ck(ld1,A,B,cf);

        WP[1]=1;
        pos=NULL;
        for (k=ld1;k>0;k--)
           pos=WP[1].addterm(cf[k],k+1,pos);
        for (v=2;v<=ld;v++)
            WP[v]=modxn(WP[v-1]*WP[1],ld+1);
//
// WPv have understood multiplier x^-v
//            
        get_ck(ld1,atilde,btilde,cft);

        Y=0;
        pos=NULL;
        for (k=ld1;k>0;k--)
            pos=Y.addterm((lp*cf[k]-cft[k])/(ZZn)((2*k+1)*(2*k+2)),k+1,pos);
        Y.addterm(-p1,1,pos);

        RF=1;
        H=1;
        X=1;
        for (r=1;r<=ld;r++)
        {
            X=modxn(X*Y,ld+1);
            RF*=r;
            H+=(X/RF);
        }
//
//  H has understood multiplier x^-d
//
        ad=1;
        Fl=0;
        pos=Fl.addterm(ad,ld);
        for (v=ld-1;v>=0;v--)
        {
            H-=ad*WP[v+1];
            H=divxn(H,1);
            ad=H.min();
            pos=Fl.addterm(ad,v,pos);
        }

        setmod(Fl);
        MY2=Y2;
        MY4=Y4;

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

// Calculate Divisor Polynomials for small primes - Schoof 1985 p.485
// This time mod the new (small) modulus Fl
// Set the first few by hand....
    
        PolyMod Pf[300],P2f[300],P3f[300];
        Pf[0]=0; Pf[1]=1; Pf[2]=2; Pf[3]=0; Pf[4]=0;

        P2f[1]=1; P3f[1]=1;

        P2f[2]=Pf[2]*Pf[2];
        P3f[2]=P2f[2]*Pf[2];

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

        P2f[3]=Pf[3]*Pf[3];
        P3f[3]=P2f[3]*Pf[3];

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

        P2f[4]=Pf[4]*Pf[4];
        P3f[4]=P2f[4]*Pf[4];
        lower=5;
        
//
// Now looking for value of lambda which satisfies
// (X^P,Y^P) = lambda.(XX,YY). 
// 3. Page 118, Algorithm 7.9
//
// Note that it appears to be sufficient to only compare the Y coordinates (!?)
// For a justification see page 120 of 3. 
// Thank you SYSTRAN translation service! (www.altavista.com)
//
        good[nl]=lp;
        cout << "NP mod " << lp << " = " << flush;
        for (lambda=1;lambda<=(lp-1)/2;lambda++)
        {
            int res=0;
            PolyMod Ry,Ty;
            tau=(lambda+invers(lambda,lp)*p)%lp;

            k=(lp+tau*tau-(4*p)%lp)%lp;
            if (jac(k,lp)!=discrim) continue; 
//
//  Possible values of tau could be eliminated here by an application of 
//  Atkin's algorithm....  
//  
            cout << setw(3) << (p+1-tau)%lp << flush; 

      // This "loop" is usually executed just once
            for (jj=lower;jj<=lambda+2;jj++)
            { // different for even and odd 
                if (jj%2==1)     // 2 mod-muls
                { 
                    n=(jj-1)/2;
                    if (n%2==0)
                        Pf[jj]=Pf[n+2]*P3f[n]*MY4-P3f[n+1]*Pf[n-1];
                    else
                        Pf[jj]=Pf[n+2]*P3f[n]-MY4*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
                P3f[jj]=P2f[jj]*Pf[jj];    // cube
            }
            if (lambda+3>lower) lower=lambda+3;

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

            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 (iszero(Ty-Ry)) res=1;
            if (iszero(Ty+Ry)) 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[nl]=tau;
                if ((p+1-tau)%lp==0)
                {
                    cout << " ***" << endl;
                    if (search) escape=TRUE;
                }
                else cout << endl;
                break;
            }
            cout << "\b\b\b";
        }
        nl++;        
        accum*=lp;
        if (escape) break;
    }
    Modulus.clear();

    if (escape) {b+=1; continue;}

    Crt CRT(nl,good);
    Big order,ordermod;
    ordermod=accum;
    order=(p+1-CRT.eval(t))%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;
        ofile << nrp << 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 << x << endl;
        ofile << y << endl;
        mip->IOBASE=10;
    }
    if (p==nrp) 
    {
        cout << "WARNING: Curve is anomalous" << endl;
        return 0;
    }
// check MOV condition for curves of Cryptographic interest
    d=1;
    for (i=0;i<50;i++)
    {
        d=modmult(d,p,nrp);
        if (d==1) 
        {
           cout << "WARNING: Curve fails MOV condition" << endl;
           return 0;
        }
    }

    return 0;
}