//
// read in all pre-processed bivariate modular polynomials
//
    modulo(p);   // set modulus
    l[0]=2;
    max=0;
    for (i=1;;i++)
    {
        mueller >> lp;
        if (mueller.eof()) break;
        max=i;
        l[i]=lp;
        cout << setw(3) << lp << flush;
        posXY=NULL;
        Gl[i].clear();
        forever 
        {
            mueller >> c >> nx >> ny;
            posXY=Gl[i].addterm((ZZn)c,nx,ny,posXY);
            if (nx==0 && ny==0) break;
        }
        cout << "\b\b\b" << flush;
    }

    mueller.close();

// 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;
    }
    
 // Calculate j-invariant

    j=(-1728*64*A*A*A)/delta;

    if (j==0 || j==1728)
    {
        cout << "Not Allowed! j-invariant = 0 or 1728" << endl;
        if (search) {b+=1; continue; }
        else return 0;
    }


// 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);
    YY=0;
    YY.addterm((ZZn)-1,0);      // Point (X,-Y)
    setmod(Y2);

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

    nl=0; 
    accum=1;           // accumulated product of good primes

    Poly zero;
    PolyMod one,XT,YT,ZT,XL,YL,ZL,ZL2,ZT2;
    one=1;         // polynomial = 1
    zero=0;        // polynomial = 0

    PolyXY dGx,dGy,dGxx,dGxy,dGyy;
    ZZn E0b,E0bd,E2bs,E4b,E6b,Dg,Dj,Dgd,Djd,Dgs,Djs,jl,jld,p1,jd;
    ZZn E4bl,E6bl,deltal,atilde,btilde,gd,f,fd,s,Eg,Ej,Exy;
    int r,v,ld,ld1,discrim;
    ZZn M,cf[500],cft[500],ad,K,RF;
    Poly Fl,T,WP[500],H,X,Y,R;
    term *pos;

    E4b=-(A/3);
    E6b=-(B/2);
    delta=(E4b*E4b*E4b-E6b*E6b)/1728;

//
// find out how many bits we are going to need
// before Kangaroos can take over
//
    first=5;
    sl[0]=3; sl[1]=5; sl[2]=7; sl[3]=8; sl[4]=9; sl[5]=0; // do low prime powers
    SCHP=9;

    if (pbits<=256) d=pow((Big)2,64);
    else            d=pow((Big)2,72); // kangaroos do more work

    d=sqrt(p/d);

    escape=FALSE;


// Calculate Divisor Polynomials for small primes - 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];
    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];

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

    for (jj=5;jj<=SCHP+1;jj++)
    { // different for even and odd 
        if (jj%2==1)
        { 
            n=(jj-1)/2;
            if (n%2==0)
                P[jj]=P[n+2]*P3[n]*Y4-P3[n+1]*P[n-1];
            else
                P[jj]=P[n+2]*P3[n]-Y4*P3[n+1]*P[n-1];
        }
        else
        {
            n=jj/2;
            P[jj]=P[n]*(P[n+2]*P2[n-1]-P[n-2]*P2[n+1])/(ZZn)2;
      
        }
        if (jj <= 1+(SCHP+1)/2)
        { // precalculate for later
            P2[jj]=P[jj]*P[jj];
            P3[jj]=P2[jj]*P[jj];
        }
    }

//
// Schoof's original method for small primes
//

    for (i=0;;i++)
    {
        lp=sl[i];
        if (lp==0) break;

        if (lp>=first)
        {
            good[nl]=lp;
            accum*=lp;
        }
        k=p%lp;

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

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

        cout << "\b\b\b\bX^PP" << flush;
        XPP=compose(XP,XP);   // this is faster
        cout << "\b\b\b\bY^PP" << flush;
        YPP=YP*compose(YP,XP);
        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^(P*P),Y^(P*P)) + 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);

// 
// 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)
            if (lp>=first)  t[nl++]=0;
            cout << setw(3) << (p+1)%lp;
            if ((p+1)%lp==0)
            {      
                cout << " ***" << endl;
                if (search) {escape=TRUE; break;}
            }
            else cout << endl;
            continue;         
        }
//
// Now keep finding tau.(X^P,Y^P), until equality is detected
// This is very simple!
//
        XL=XP;
        YL=YP;
        ZL=1;
        ZT2=ZT*ZT;
        for (tau=1;tau<=(lp/2);tau++)
        { 
            cout << setw(3) << (p+1-tau)%lp << flush;     

            ZL2=ZL*ZL;
            if (iszero(XT*ZL2-ZT2*XL))   // LHS == RHS??
            { // LHS = RHS
                if (!iszero(YT*ZL2*ZL-YL*ZT*ZT2))
                { // point doubled - change sign
                    tau=lp-tau;
                    cout << "\b\b\b";
                    cout << setw(3) << (p+1-tau)%lp << flush;     
                }
                if (lp>=first) t[nl++]=tau;
                if ((p+1-tau)%lp==0)
                {
                    cout << " ***" << endl;
                    if (search) escape=TRUE;
                }
                else cout << endl;
                break;
            }
            elliptic_add(XL,YL,ZL,XP,YP);
            cout << "\b\b\b";            
        }
        if (escape) break;       
    }

    if (!escape) for (i=1;accum<=d;i++)      
    {
        if (i>max)
        {
            cout << "WARNING: Ran out of Modular Polynomials!" << endl;
            break;
        }
        lp=l[i];

        if (lp<=SCHP) continue;

        k=p%lp;
        for (is=1;;is++)
            if (is*(lp-1)%12==0) break;

        el=lp;
        s=is;

// Get next Modular Polynomial

        MP=Gl[i];

// Evaluate bivariate polynomial at Y=j-invariant
// and use as polynomial modulus

        F=MP.F(j);
        setmod(F);
        cout << setw(3) << lp << flush;
        XP=pow(XX,p);
//
// Determine "Splitting type" 
//
        cout << "\b\b\bGCD" << flush;
        G=gcd(XP-XX);

        if (degree(G)==lp+1) 
        {
            cout << "\b\b\b" << flush;
            continue;      // pathological case
        }
        if (degree(G)==0)  // Atkin Prime
        {
            if (!atkin && lp>100)    // Don't process large Atkin Primes
            {
                cout << "\b\b\b" << flush;
                continue;
            }
            BOOL useful=FALSE;
            cout << "\b\b\b" << flush;
            cout << "ATK" << flush;

            PolyMod u[20];
            int max_r,lim=1;
            u[0]=XP;  
            u[1]=compose(u[0],u[0]);

//
// The code for processing Atkin Primes is in here, but currently largely 
// unused. However the simplest case is used, as it suggests only one 
// value for NP mod lp, and so can be used just like an Elkies prime
// 
//
            if (atkin) max_r=lp+1;
            else       max_r=2;
            for (r=2;r<=max_r;r++)
            {
                PolyMod C;
                int kk,m;
                BOOL first;
                if ((lp+1)%r!=0) continue;    // only keep good ones!
                v=jac(k,lp);     // check Schoof Prop. 6.3
                jj=(lp+1)/r;
                if (jj%2==0 && v==(-1)) continue;
                if (jj%2==1 && v==1) continue;
       //         if (phi(r)>8) continue;       // > 8 candidates
       //         if (lp<30 && phi(r)>2) continue;
       //         if (lp<60 && phi(r)>4) continue;
                kk=r; m=0;
                first=TRUE;
//
// Right-to-Left Power Composition - find X^(P^r)
//
                forever
                {
                    if (kk%2!=0)
                    {
                        if (first) C=u[m]; 
                        else       C=compose(u[m],C);
                        first=FALSE;
                    }
                    kk/=2;
                    if (kk==0) break;
                    m++;
                    if (m>lim) 
                    { // remember them for next time