/*
 * cm.cpp - Finding an elliptic curve and point of nearly prime order
 * See IEEE 1363 Annex A for documentation!
 *
 * !!! New much improved version - March 2002
 * !!! Now tested for D up to 10^7
 *
 * !!! New faster version - uses Floats instead of Flashs - November 2003
 * !!! Now 100 times faster!
 * !!! Now tested for D up to and beyond 10^9 
 *
 * !!! New invariants - gamma2, w3, w5, w7 and w13. Thanks to Marcel Martin.
 *
 * Sometimes its better to use the -IEEE flag, as this can often be faster 
 * than the Gamma2 invariant, as although it requires a bigger "class number", 
 * it needs less precision.
 *
 * Uses functions from the MIRACL multiprecision library, specifically
 * classes:-
 * Float   - Big floating point 
 * Complex - Big Complex float
 * Big     - Big integer
 * ZZn     - Big integers mod an integer
 * FPoly   - Big Float polynomial
 * Poly    - Big ZZn polynomial
 *
 * Written by Mike Scott, Dublin, Ireland. March 1998 - March 2004
 *
 * Full MIRACL source is available from 
 * ftp.computing.dcu.ie/pub/crypto/miracl.zip
 *
 * MIRACL is a shareware source code product.
 * However it is free for educational and non-profit making use
 * Queries to mike@computing.dcu.ie
 * Web page http://indigo.ie/~mscott
 */

#include <iostream>
#include <fstream>
#include <cmath>
#include <cstring>
#include "ecn.h"
#include "complex.h"
#include "flpoly.h"
#include "poly.h"

using namespace std;

miracl *mip;

FPoly T[25];  // Reduced class Polynomial. 
static char *s;
BOOL fout,suppress,three;

// F(z) Function A.13.3

Complex Fz(Complex z)
{
    Complex t;
    int sign=1;
    Complex osum,sum=(Float)1;

    if (z.iszero()) return sum;

    Complex zi=z;
    Complex zj=z*z;
    Complex r=z;
    Complex z3=zj*z;

    forever
    { // do 2 terms at a time....
        t=zi+zj;
        osum=sum;
        if (sign) sum-=t;
        else      sum+=t;

        if (sum==osum) break;

        r*=z3;        
        zi*=r; zj*=r; zj*=z;

        sign=1-sign;
    }

    return sum;
}

// Fj(A,B,C) function A.13.3

Complex F(int j,Big A,Big B,Big C,Big D,int N)
{
    Complex t,theta24,theta,theta2,theta3,theta6;
    Float sd;

    if (j>=3)
    { // Gamma 2 and Morain's invariants
        sd=-sqrt((Float)D);
        t=Complex(sd*fpi(),(fpi()*(Float)B));    
        t/=((Float)A);
        theta=exp(t);
        if (j==3)
        {
            theta3=theta;
            theta3*=theta3;
            theta3*=theta;
            theta6=theta3;
            theta6*=theta6;
            t=theta*pow((Fz(theta6)/Fz(theta3)),8);
            return (256*t*t*t+1)/t;
        }
        if (j==4)
        {
              return (pow(Fz(theta)/Fz(pow(theta,N)),24/(N-1))/theta);   
        }
        return 0;
    }

    sd=-sqrt((Float)D);
    t=Complex(sd*fpi(),(fpi()*(Float)B));    

    t/=(Float)(24*A);
    theta24=exp(t);                                // theta^(1/24)
    theta=pow(theta24,24);                         // theta

    if (j!=2) 
        t=recip(theta24);      // -24th root
    else  
        t=theta24*theta24;     // 12th root

    theta2=theta;
    theta2*=theta2;

    if (j==0) return (t*Fz(-theta)/Fz(theta2));
    if (j==1) return (t*Fz(theta)/Fz(theta2));
    if (j==2) return (sqrt((Float)2)*t*Fz(theta2*theta2)/Fz(theta2));
    return 0;
}

int geti(Big D)
{
    Big d=D%8;
    if (d==1 || d==2 || d==6 || d==7) return 3;
    if (d==3)
    {
       if (D%3==0) return 2;
       else        return 0;
    } 
    if (d==5) return 6;
    return 0;
}

int getk(Big D)
{
    Big d=D%8;
    if (d==1 || d==2 || d==6) return 2;
    if (d==3 || d==7) return 1;
    if (d==5) return 4;
    return 0; 
}

int getN(Big D)
{
    if (D%13==0) return 13;
    if (D%7==0)  return 7;
    if (D%5==0)  return 5;
    if (D%3==0)  return 3;
    return 2;
}

// A.13.3

void class_poly(Complex& lam,Float *Fi2,Big A,Big B,Big C,Big D,BOOL conj,BOOL P1363)
{
    Big ac,l,t;
    int i,j,k,g,e,m,n,dm8;
    Complex cinv;

    if (P1363)
    {
        g=1;
        if (D%3==0) g=3;

        ac=A*C;
        if (ac%2==1) 
        {
            j=0;
            l=A-C+A*A*C;
        }
        if (C%2==0)  j=1;
        if (A%2==0)  j=2;

        if (A%2==0)
        {
            t=(C*C-1)/8;
            if (t%2==0) m=1;
            else        m=-1;
        }
        else
        {
            t=(A*A-1)/8;
            if (t%2==0) m=1;
            else        m=-1;
        }
    
        dm8=D%8;
        i=geti(D);
        k=getk(D);
        switch (dm8)
        {
        case 1: 
        case 2:    n=m;
                   if (C%2==0) l=A+2*C-A*C*C;
                   if (A%2==0) l=A-C-A*C*C;
                   break;
        case 3:    if (ac%2==1) n=1;
                   else         n=-m;
                   if (C%2==0) l=A+2*C-A*C*C;
                   if (A%2==0) l=A-C+5*A*C*C;
                   break;
        case 5:    n=1;
                   if (C%2==0) l=A-C+A*A*C;
                   if (A%2==0) l=A-C-A*C*C;
                   break;
        case 6:    n=m;
                   if (C%2==0) l=A+2*C-A*C*C;
                   if (A%2==0) l=A-C-A*C*C;
                   break;
        case 7:    if (ac%2==0) n=1;
                   else         n=m;
                   if (C%2==0) l=A+2*C-A*C*C;
                   if (A%2==0) l=A-C-A*C*C;
                   break;
               
        default: break;
        }
        e=(k*B*l)%48;
        if (e<0) e+=48;
        cinv=pow(lam,e);
        cinv*=(n*Fi2[i]);
        cinv=pow(cinv*pow(F(j,A,B,C,D,0),k),g);
    }
    else
    {
        int N=getN(D);   
// adjust A and B

        if (N==2)
        {
            j=3;
            if (A%3!=0)
            {
                if (B%3!=0)
                {
                    if ((B+A+A)%3!=0) B+=(4*A);
                    else              B+=(A+A);
                }
            }
            else
            {
                if (B%3!=0)
                {
                    if (C%3!=0)
                    {
                        if ((C+B)%3!=0) B+=(4*A);
                        else            B+=(A+A);
                    }
                }
            }
            A*=3;
        }
        else
        {
            j=4;
            if ((A%N)==0) 
            {
                A=C;
                B=-B;
            }
            while (B%N!=0) B+=(A+A);
            A*=N;
        }
        cinv=F(j,A,B,C,D,N);
    }


 // multiply polynomial by new term(s)

    FPoly F;
    if (conj)
    { // conjugate pair
      // t^2-2a+(a^2+b^2) , where cinv=a+ib
        F.addterm((Float)1,2);
        F.addterm(-2*real(cinv),1);
        F.addterm(real(cinv)*real(cinv)+imaginary(cinv)*imaginary(cinv),0);
    }
    else 
    { // t-cinv
        F.addterm((Float)1,1);
        F.addterm(-real(cinv),0);

// store as a linear polynomial, or combine 2 to make a quadratic
        if (T[0].iszero())
        {
            T[0]=F;
            return;
        }
        else 
        {
            F=T[0]*F;      // got a quadratic
            T[0].clear();
        }
    }

// accumulate Polynomial as 2^m degree components
// This allows the use of karatsuba via the "special" function
// This is the time critical bit....

    for (i=1;;i++)
    {
        if (T[i].iszero())
        {
            T[i]=F;             // store this 2^i degree polynomial
            break;
        }
        else
        {
            F=special(T[i],F);  // e.g. if i=1 two quadratics make a quartic..
            T[i].clear();
        }
    }
}

// A.13.2
// Set P1363 to False to calculate class number a la Cohen

int groups(Complex& lam,Float *Fi2,unsigned long D,BOOL doit,BOOL P1363)
{
    unsigned long s,t,A,C,B,TB,lim;
    int cn=0;
    s=lsqrt(D/3,1);
    for (B=0;B<=s;B+=1)
    {
//  cout << "B= " << B << " s= " << s << endl;
        t=D+B*B;
        if (!P1363)
        {
            if (t%4!=0) continue;
            else t/=4;
        }
// cout << "t= " << t << endl;
        lim=lsqrt(t,1);
// cout << "lim= " << lim << endl;
        if (P1363) A=2*B;
        else       A=B;
        if (A==0) A+=1;
        for(;;)
        {
            while (t%A!=0) 
            {
                A+=1;
                if (A>lim) break;
            }

            if (A>lim) break;
            C=t/A;
           
            TB=B;
            if (P1363) TB*=2;
            if (lgcd(lgcd(A,TB),C)==1)
            { // output more class group members
                BOOL conj;
                if (TB>0 && C>A && A>TB) 
                { 
                    conj=TRUE;
                    cn+=2;
                    if (doit)
                    { 
                        if (!suppress) cout << ".." << flush; 
                    }  
                }
                else
                {
                    conj=FALSE;
                    cn+=1;
                    if (doit)
                    {   
                        if (!suppress) cout << "." << flush; 
                    }
                } 
                if (doit) class_poly(lam,Fi2,(Big)A,(Big)B,(Big)C,(Big)D,conj,P1363);
            }
            A+=1;
        }        
    }
    return cn;         // class number
}

// check congruence conditions A14.2.1

BOOL isaD(unsigned long d,int pm8,Big k)
{
    unsigned int dm8;
    unsigned long i;
    BOOL sqr;
    dm8=d%8;
    if (k==1 && dm8!=3) return FALSE;
    if ((k==2 || k==3) && dm8==7) return FALSE;
    if (pm8==3 && (dm8==1 || dm8==4 || dm8==5 || dm8==6)) return FALSE;
    if (pm8==5 && dm8%2==0) return FALSE;
    if (pm8==7 && (dm8==1 || dm8==2 || dm8==4 || dm8==5)) return FALSE;
    sqr=FALSE;
    for (i=2;;i++)
    {
        if (d%(i*i)==0)
        {
            sqr=TRUE;
            break;
        }
        if (i*i>d) break;
    }
    if (sqr) return FALSE;
    return TRUE;
}

// Testing for CM discriminants A.14.2.2

Big floor(Big N,Big D)
{
   Big R;
   if (N==0) return 0;
   if (N>0 && D>0) return N/D;
   if (N<0 && D<0) return (-N)/(-D);
   R=N/D;
   if (N%D!=0) R-=1;
   return R;
}

BOOL isacm(Big p,unsigned long D,Big &W,Big &V)
{
    Big B2,A,B,C,t,X,Y,ld,delta;
    B=sqrt(p-(Big)D,p);
    A=p;
    C=(B*B+(Big)D)/p;
    X=1;
    Y=0;
    ld=0;

    while (1)
    {
        if (C>=A)
        {
            B2=2*B;
            if (B2<0) B2=-B2;
            if (A>=B2) break;
        } 
        delta=floor(2*B+C,2*C);