}

        middle=(upper+lower)/2;
        if (ordermod>1)
            middle+=(ordermod+order-middle%ordermod);  // middle must be
                                                       // order mod ordermod

        for (i=0;i<HERD;i++) start[i]=middle+13*ordermod*i; // tame ones
        for (i=0;i<HERD;i++) start[HERD+i]=13*ordermod*i;   // wild ones

        for (i=0;i<2*HERD;i++) 
        {
            K[i]=start[i]*P;         // on your marks ....
            D[i]=0;                  // distance counter
        }
        cout << "Releasing " << HERD << " Tame and " << HERD << " Wild Kangaroos\n";

        nt=0; nw=0; cw=0; ct=0;
        collision=FALSE;  abort=FALSE;
        forever
        { 
            for (jj=0;jj<HERD;jj++)
            {
                K[jj].get(txc);
                i=txc%m;  /* random function */

                if (txc%distinguished==0)
                {
                    if (nt>=STORE)
                    {
                        abort=TRUE;
                        break;
                    }
                    cout << "." << flush;
                    tame[nt]=K[jj];
                    tdist[nt]=D[jj];
                    tname[nt]=jj;
                    for (k=0;k<nw;k++)
                    {
                        if (wild[k]==tame[nt])
                        {
                            ct=nt; cw=k;
                            collision=TRUE;
                            break;
                        }
                    }
                    if (collision) break;
                    nt++;
                }
                D[jj]+=distance[i];
                TE[jj]=table[i];
            }
            if (collision || abort) break;
            for (jj=HERD;jj<2*HERD;jj++)
            {
                K[jj].get(wxc);
                j=wxc%m;

                if (wxc%distinguished==0)
                {
                    if (nw>=STORE)
                    {
                        abort=TRUE;
                        break;
                    }
                    cout << "." << flush;
                    wild[nw]=K[jj];
                    wdist[nw]=D[jj];
                    wname[nw]=jj;
                    for (k=0;k<nt;k++)
                    {
                        if (tame[k]==wild[nw]) 
                        {
                            ct=k; cw=nw;
                            collision=TRUE;
                            break;
                        }
                    }
                    if (collision) break;
                    nw++;
                }
                D[jj]+=distance[j];
                TE[jj]=table[j];
            }
            if (collision || abort) break;

            multi_add(2*HERD,TE,K);   // jump together - its faster
        }
        cout << endl;
        if (abort)
        {
            cout << "Failed - this should be rare! - trying again" << endl;
            continue;
        }
        nrp=start[tname[ct]]-start[wname[cw]]+tdist[ct]-wdist[cw]; 
                                                   // = order mod ordermod
        G=P;
        G*=nrp;
        if (G!=ZERO)
        {
            cout << "Sanity Check Failed. Please report to mike@compapp.dcu.ie" << endl;
            exit(0);
        }
        if (prime(nrp))
        {
            cout << "NP= " << nrp << endl;
            cout << "NP is Prime!" << endl;
            break;
        }

// final checks....
        real_order=nrp; i=0; 
        forever
        {
            sp=get_mip()->PRIMES[i];
            if (sp==0) break;
            if (real_order%sp==0)
            {
                G=P;
                G*=(real_order/sp);
                if (G==ZERO) 
                {
                    real_order/=sp;
                    continue;
                }
            }
            i++;
        }   
        if (real_order <= 4*sqrt(p))
        { 
            cout << "Low Order point used - trying again" << endl;
            continue;
        }
        real_order=nrp;
        for (i=0;(sp=get_mip()->PRIMES[i])!=0;i++)
            while (real_order%sp==0) real_order/=sp;
        if (real_order==1)
        { // all factors of nrp were considered
            cout << "NP= " << nrp << endl;
            break;
        }
        if (prime(real_order))
        { // all factors of NP except for one last one....
            G=P;
            G*=(nrp/real_order);
            if (G==ZERO) 
            {
                cout << "Failed - trying again" << endl;
                continue;
            }
            else 
            {
                cout << "NP= " << nrp << endl;
                break;
            }
        }
// Couldn't be bothered factoring nrp completely
// Probably not an interesting curve for Cryptographic purposes anyway.....
// But if 10 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<10;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;
}

//
// Coefficients Ck from Mueller's lemma 6.2
// 3. page 90
// 4. page 127-128
//

void get_ck(int terms,ZZn a,ZZn b,ZZn *c)
{
    int k,h;
    if (terms==0) return;
    c[1]=-a/5;
    if (terms==1) return;
    c[2]=-b/7;
    for (k=3;k<=terms;k++)
    {
        c[k]=0;
        for (h=1;h<=k-2;h++) c[k]+=c[h]*c[k-1-h];
        c[k]*=((ZZn)3/(ZZn)((k-2)*(2*k+3)));
    }
}

//
// quadratic field arithmetic
// needed for Atkin Primes
//

void mulquad(int p,int qnr,int x,int y,int& a,int& b)
{
    int olda=a;
    a=(a*x+b*y*qnr)%p;
    b=(olda*y+b*x)%p;     
}

void powquad(int p,int qnr,int x,int y,int e,int& a,int& b)
{
    int k=e;
    a=1;
    b=0;
    if (k==0) return;

    for(;;)
    {
        if (k%2!=0)
            mulquad(p,qnr,x,y,a,b);        
        k/=2;
        if (k==0) return;
        mulquad(p,qnr,x,y,x,y);
    }
}

//
// Euler Totient function
//

int phi(int n)
{
    int i,r=1;
    for (i=2;i<n;i++) if (igcd(i,n)==1) r++;
    return r;
}

int main(int argc,char **argv)
{
    ofstream ofile;
    ifstream mueller;
    int SCHP,first,max,lower,ip,parity,pbits,lp,i,jj,is,n,nx,ny,nl;
    int sl[8],k,tau,lambda;
    mr_utype good[100],l[100],t[100];
    Big a,b,c,p,nrp,x,y,d,accum;
    PolyMod XX,YY,XP,XPP,YP,YPP;
    PolyXY MP,Gl[200];
    termXY *posXY;
    Poly F,G,P[500],P2[500],P3[500],Y2,Y4;
    miracl *mip=&precision;
    BOOL escape,search,fout,gotI,gotA,gotB,atkin;
    ZZn j,g,qb,qc,delta,el;
    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 << "sea <A> <B> <-i file>" << endl;
        cout << "where the input file contains pre-processed modular polynomials" << endl;
        cout << "for a particular prime modulus P" << endl;
        cout << "e.g. sea -3 35317045537 -i mueller.pol" << endl;
        cout << "To input A and B in Hex, precede with -h" << endl;
        cout << "To output to a file, use flag -o <filename>" << endl;
        cout << "To observe Atkin prime processing, use flag -a" << endl;
        cout << "NOTE: Atkin prime information is not currently used" << 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 << "http://indigo.ie/~mscott for details" << endl;
        cout << "or email mscott@indigo.ie" << endl;
        return 0;
    }

    ip=0;
    gprime(10000);                 // generate small primes < 1000
    search=fout=gotI=gotA=gotB=atkin=FALSE;
    a=0; b=0;

// Interpret command line
    Base=10;
    while (ip<argc)
    {
        if (!fout && 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 (!gotI && strcmp(argv[ip],"-i")==0)
        {
            ip++;
            if (ip<argc)
            {
                gotI=TRUE;

                mueller.open(argv[ip],ios::in);
                if (!mueller)
                {
                    cout << "input file " << argv[ip] << " could not be opened" << endl;
                    return 0;
                }
                ip++;
                continue;
            }
            else
            {
                cout << "Error in command line" << endl;
                return 0;
            }
        }

        if (strcmp(argv[ip],"-s")==0)
        {
            ip++;
            search=TRUE;
            continue;
        }

        if (strcmp(argv[ip],"-a")==0)
        {
            ip++;
            atkin=TRUE;
            continue;
        }

        if (strcmp(argv[ip],"-h")==0)
        {
            ip++;
            Base=16;
            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 (!gotI || !gotA || !gotB)
    {
        cout << "Error in command line" << endl;
        return 0;
    }

// get prime modulus from .pol file

    mueller >> p;
    pbits=bits(p);
    if (pbits<64)
    {
        cout << "Prime modulus too small for this program!" << endl;
        cout << "Use SCHOOF program (ftp://ftp.compapp.dcu.ie/pub/crypto/schoof.exe)" << endl;
        exit(0);
    }
    cout << "P= " << p << endl;
    cout << "P mod 24 = " << p%24 << endl;
    cout << "P is " << pbits << " bits long" << endl;
    cout << "Reading in pre-processed Modular Polynomials... " << endl;