// are we stuck in hopeless loop? This can't happen now - thanks Marcel

        if (delta==0 && ld==0) return FALSE;
        ld=delta;
        t=X;
        X=(delta*X+Y);
        Y=-t;

        t=C;

        C=(A+C*delta*delta-2*B*delta);
        B=(t*delta-B);
        A=t;        
    }
    if (D==11 && A==3)
    {
        t=X;
        X=Y;
        Y=-t;
        B=-B;
        t=C;
        C=A;
        A=t;        
    }      
    if (D==1 || D==3) 
    {
       W=2*X;
       V=2*Y;
       return TRUE;
    }
    else V=0;

    if (A==1) 
    {
        W=2*X;
        return TRUE;
    }

    if (A==4) 
    {
        W=4*X+B*Y;
        return TRUE;
    }

    return FALSE;
}


// Constructing a Curve (and Point if order is prime)
// A.14.4

BOOL get_curve(Complex& lam,Float *Fi2,ofstream& ofile,Big r,Big other_r,Big ord,unsigned long D,Big p,Big W,int m,BOOL always)
{
    Poly polly;
    Big A0,B0,k;
    BOOL unstable;
    char c;
    int i,j,terms,class_number;
    BOOL P1363,pord=prime(ord);

    P1363=TRUE;
    if (!always && D>3 && D%8==3) P1363=FALSE;

    k=r;
    while (k%ord==0)  k/=ord;

    if (!suppress) 
    {
        if (D>1000 && D%3==0) cout  << "D= " << D << " (Divisible by 3 - might need extra precision)" << endl;
        else                  cout  << "D= " << D << endl;
        if (P1363) cout << "IEEE-1363 invariant" << endl;
        else
        {
            switch (getN(D))
            {
            case 2:  cout << "Gamma2 invariant" << endl;
                     break;
            case 3:  cout << "w3 invariant" << endl;
                     break;
            case 5:  cout << "w5 invariant" << endl;
                     break;
            case 7:  cout << "w7 invariant" << endl;
                     break;
            case 13:  cout << "w13 invariant" << endl;
                     break;
            }
        }
        cout << "D mod 24 = " << D%24 << endl;
        if (prime(ord)) cout << "K= " << k << endl;
    }

    class_number=groups(lam,Fi2,D,FALSE,P1363);
    
    cout << "class number= " << class_number << endl;

    cout << "Group order= " << ord << endl;
    cout << "do you want to proceed with this one (Y/N)? " ;
    cin >> c;
    if (c=='N' || c=='n') 
    {
        if (!suppress) cout << "finding a curve..." << endl;
        return FALSE;
    }

    modulo(p);
                                   
// A.14.4.1
    
    if (D==1)
    {
        A0=1; B0=0;
    }
    if (D==3)
    {
        A0=0; B0=1;
    }

    if (D!=1 && D!=3)
    {
        FPoly Acc;

        for (i=0;i<25;i++) T[i].clear();
        if (!suppress) cout << "constructing polynomial";
        groups(lam,Fi2,D,TRUE,P1363);

// Construct Polynomial from its 2^m degree components..

        for (i=24;i>=0;i--)
        {
            if (!T[i].iszero())
            {
                Acc=T[i];              // find the first component..
                T[i].clear();
                break;
            }  
        }
        if (i>0)
        {
            for (j=i-1;j>0;j--)
            {
                if (!T[j].iszero())
                {
                    Acc=special(Acc,T[j]); // special karatsuba function
                    T[j].clear();          // multiply into accumulator
                }   
            }
            if (!T[0].iszero())
            {                              // check for a left-over linear poly
                Acc=Acc*T[0];
                T[0].clear();
            }
        }
        for (i=0;i<25;i++) T[i].clear();

        terms=degree(Acc);
        Float f,rem;
        Big whole;
        int nbits,maxbits=0;

        unstable=FALSE;
        for (i=terms;i>=0;i--)
        {
            f=Acc.coeff(i);
            if (f>0)
                 f+=makefloat(1,10000);
            else f-=makefloat(1,10000);
            whole=f.trunc(&rem);
            nbits=bits(whole);
            if (nbits>maxbits) maxbits=nbits;
            polly.addterm((ZZn)whole,i);
            if (fabs(rem)>makefloat(1,100)) 
            {
                unstable=TRUE; 
                break; 
            }
        }
        Acc.clear();
        if (!suppress) cout << endl;
        if (unstable) 
        {
            if (!suppress) cout << "Curve abandoned - numerical instability!" << endl;
            if (!suppress) cout << "Curve abandoned - double MIRACL precision and try again!" << endl;
            if (!suppress) cout << "finding a curve..." << endl;
            return FALSE;
        }
        if (!suppress) 
        {
            cout << polly << endl;
            cout << "Maximum precision required in bits= " << maxbits << endl;
        }
    }

// save space with smaller miracl

    mirexit();
    mip=mirsys(128,0);
    modulo(p);

    ECn pt,G;
    Big a,b,x,y;
    Big w,eps;
    int at;
    Poly g,spolly=polly;     // smaller polly
    polly.clear();
    forever
    {
        if (D!=1 && D!=3)
        {
            if (!suppress) cout << "Factoring polynomial of degree " << degree(spolly) << " ....." << endl;

            if (P1363)
            {
                if (W%2==0)
                {
                    ZZn V;
                    g=factor(spolly,1);
                    V=-g.coeff(0);
                    V=pow(V,24/(lgcd(D,3)*getk(D)));
                    V*=pow((ZZn)2,(4*geti(D))/getk(D));
                    if (D%2!=0) V=-V;
                    A0=(Big)((ZZn)(-3)*(V+64)*(V+16));   
                    B0=(Big)((ZZn)2*(V+64)*(V+64)*(V-8));
                }
                else
                {
                    Poly V,w,w1,w2,w3,a,b;

                    g=factor(spolly,3);
                    if (D%3!=0)
                        w.addterm(-1,24);
                    else
                        w.addterm(-256,8);
                    V=w%g;
                    w.clear();
                    w1=V; w2=V; w3=V;
                    w1.addterm(64,0);
                    w2.addterm(256,0);
                    w3.addterm(-512,0);
                    a=w1*w2;
                    a.multerm(-3,0);
                    a=a%g;
                    b=w1*w1*w3;
                    b.multerm(2,0);
                    b=b%g;

                    a=((a*a)*a)%g;
                    b=(b*b)%g;
                    for (int d=degree(g)-1;d>=0;d--)
                    {
                        ZZn t,c=a.coeff(d);
                        if (c!=(ZZn)0)
                        {
                            t=b.coeff(d);
                            A0=(Big)(c*t);
                            B0=(Big)(c*t*t);
                            if (d==1) break;
                        }     
                    }
                }
            }
            else
            {
                ZZn X,J,X2,K;

                g=factor(spolly,1);
                X=-g.coeff(0);
                X2=X*X;
                switch (getN(D))
                {
                case 2:  J=X2*X; 
                         break;
                case 3:  J=(pow((X+3),3)*(X+27))/X;
                         break;
                case 5:  J=pow((X2+10*X+5),3)/X;
                         break;
                case 7:  J=(pow((X2+5*X+1),3)*(X2+13*X+49))/X;
                         break;
                case 13: J=(pow((X2*X2+7*X2*X+20*X2+19*X+1),3)*(X2+5*X+13))/X;
                default: break;
                }
                K=J/(J-1728);
                A0=-3*K;
                B0=2*K;
            }

        }

// A.14.4.2
// but try -3 first, followed by small positive values for A

        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;

        if (D!=1 && D!=3 && three && a!=-3) continue;

// 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 )
        {
                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;
        mirexit();
        int precision=(1<<m);
        mip=mirsys(precision+2,0);     // restart with new precision
        setprecision(m);
        gprime(1000000);
        mip->RPOINT=ON;


        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 #

#if defined(unix)
#define TIMES '.'
#define RAISE '^'
#else
#define TIMES '*'
#define RAISE '#'
#endif

void eval_power (Big& oldn,Big& n,char op)
{
        if (op) n=pow(oldn,toint(n));    // power(oldn,size(n),n,n);
}

void eval_product (Big& oldn,Big& n,char op)
{
        switch (op)
        {