/*
 * Factorisation de grands entiers à l'aide de courbes elliptiques.
 *
 * Ce code est sous licence GPL.
 *
 * Samuel Mimram
 */

/*
 * ECFacto is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * ECFacto is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Ocaml-mad; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * @author Samuel Mimram
 */


#include "stdafx.h"

#include "largeint.h"

#define B_incr 1
#define a_incr 17
#define nb_essais	1000
#define x_init 1
#define y_init 1

spLargeint a, b, p;
bool error;
spLargeint error_factor=new largeint;

#define O_neutre ec_point(true)

void PrintTime()
{
    char cTmp[100];

    _strtime(cTmp);
    cout << cTmp << endl;
}

class ec_point;

#define spEcpoint RefCountPtr<ec_point>

class ec_point
{
    public:
        spLargeint x;
        spLargeint y;
        bool Infty;

        ec_point(){Init(); Infty=false;};
        ec_point(spLargeint X, spLargeint Y){Init(); x=X; y=Y; Infty=false;};
        ec_point(bool isO){Init(); Infty=isO;};
        ~ec_point(){};
        void Init();
        void Print();

        friend bool operator ==(spEcpoint P1, spEcpoint P2);
        friend spEcpoint operator +(spEcpoint P1, spEcpoint P2);
        friend spEcpoint operator * (spLargeint m, spEcpoint P);
        friend spEcpoint operator * (int m, spEcpoint P);
};

inline void CopyEcpoint(spEcpoint dst, spEcpoint src){CopyLargeint(dst->x, src->x); CopyLargeint(dst->y, src->y); dst->Infty=src->Infty;};

void ec_point::Init()
{
    x=new largeint;
    y=new largeint;
}

void ec_point::Print()
{
    if(Infty)
        cout << "point : O" << endl;
    else
    {
        cout << "point : " << endl;
        x->Print();
        y->Print();
    }
}

bool operator ==(spEcpoint P1, spEcpoint P2)
{
    if(P1->Infty && P2->Infty)
        return true;
    else if(P1->Infty || P2->Infty)
        return false;
    else
        return (P1->x==P2->x && P1->y==P2->y);
}

spEcpoint operator + (spEcpoint P1, spEcpoint P2)
{
    spEcpoint P3=new ec_point;
    spLargeint lambda;

    if(P1->Infty)
        CopyEcpoint(P3, P2);
    else if(P2->Infty)
        CopyEcpoint(P3, P1);
    else if(P1->x==P2->x && P1->y==mod(-P2->y, p))
        P3->Infty=true;
    else
    {
        if(P1==P2)
        {
            if(gcd(P1->y, p)!=1)
            {
                CopyEcpoint(error_factor, P1->y);
                error=true;
            }
            lambda=mod((3*P1->x*P1->x+a)*invModN(2*P1->y, p), p);
        }
        else
        {
            if(gcd(P2->x-P1->x, p)!=1)
            {
                error_factor=P2->x-P1->x;
                error=true;
            }
            lambda=mod((P2->y-P1->y)*invModN(P2->x-P1->x, p), p);
        }
        P3->x=mod((lambda*lambda-P1->x-P2->x), p);
        P3->y=mod((lambda*(P1->x-P3->x)-P1->y), p);
    }
    return P3;
}

spEcpoint fast_scalarmult (spLargeint m, spEcpoint P)
{
    spEcpoint tmpP=new ec_point, retP=new ec_point(true);

    CopyEcpoint(tmpP, P);
    while(m!=0)
    {
        if(error)
        {
            tmpP=new ec_point(true);
            return tmpP;
        }
        if(getbit(m, 0))
            retP=retP+tmpP;
        m->shr(1);
        tmpP=2*tmpP;
    }

    return retP;
}

spEcpoint operator * (spLargeint m, spEcpoint P)
{
    spLargeint n=new largeint;
    CopyLargeint(n, m);
    return fast_scalarmult(n, P);
}

spEcpoint operator * (int m, spEcpoint P)
{
    if(m==2)
        return P+P;
    else
    {
        spLargeint tmp=new largeint(m);
        return (tmp*P);
    }
}

spLargeint calc_k(spLargeint B)
{
    spLargeint TODO;

    return lcm(B, B-1);	//max(lcm(B, B-1), 2);
}

spLargeint TrouveFacteur(spLargeint n, spLargeint B)
{
    spEcpoint P, Q;
    spLargeint k=new largeint, tmp=new largeint, ret;
    int a_essais;
    bool nouvelle_courbe=false;

    error=false;
    if((n%2)==0)
    {
        PrintTime();
        cout << "facteur : 2";
        cout.flush();
        cout << "\n\n";
        cout.flush();
        ret=new largeint(2);
        return ret;
    }
    if((n%3)==0)
    {
        PrintTime();
        cout << "facteur : 3";
        cout.flush();
        cout << "\n\n";
        cout.flush();
        ret=new largeint(3);
        return ret;
    }

    p=n;
    k=calc_k(B);
    b=new largeint(1);
    a=new largeint(1); a_essais=0;

    P=new ec_point(new largeint(x_init), new largeint(y_init));
    while(1)
    {
        b=(P->y*P->y)-(P->x*P->x*P->x)-(a*P->x);

        tmp=gcd(4*a*a*a+27*b*b, n);
        while(tmp!=1)
        {
            if(tmp==n)
            {
                a=a+new largeint(a_incr);
                b=(P->y*P->y)-(P->x*P->x*P->x)-(a*P->x);
                tmp=gcd(4*a*a*a+27*b*b, n);
            }
            else
                return tmp;
        }
        if(a_essais>nb_essais)
            nouvelle_courbe=true;
        Q=k*P;
        if(error)
        {
            spLargeint g=gcd(error_factor, n);

            if(1<g && g<n)
            {
                PrintTime();
                cout << "facteur : ";
                cout.flush();
                g->Print();
                cout.flush();
                cout << "a : ";
                cout.flush();
                a->Print();
                cout.flush();
                cout << "b : ";
                cout.flush();
                b->Print();
                cout.flush();
                cout << "p : ";
                cout.flush();
                p->Print();
                cout.flush();
                cout << "B : ";
                cout.flush();
                B->Print();
                cout.flush();
                cout << "k : ";
                cout.flush();
                k->Print();
                cout.flush();
                cout << "\n";
                cout.flush();
                return g;
            }
            else
                nouvelle_courbe=true;
        }
        if(nouvelle_courbe)
        {
            a=new largeint(1);
            a_essais=0;
            B=B+B_incr;
            while(k==calc_k(B)){B=B+B_incr;};
            k=calc_k(B);
        }
        else
        {
            a=a+a_incr;
            a_essais++;
        }
    }

    ret=new largeint(-1);
    return ret;
}

int main(int argc, char* argv[])
{
    /*a=-1;b=0;p=1283;

      int i;
      ec_point P(121, 30);
      for(i=1;i<1285; i++)
      {
      cout << i << " : ";
      (i*P).Print();
      cout << endl;
      }
      getch();
      return 0;*/

    spLargeint n=new largeint, B=new largeint, f;

    cout << "Nombre a factoriser : ";
    cout.flush();
    n->Input();
    cout << "B-lissite : ";
    cout.flush();
    B->Input();

    cout << endl;
    PrintTime();
    cout << endl;

    f=new largeint(2);
    while(f>1)
    {
        f=TrouveFacteur(n, B);
        if(f>1)
        {
            n=n/f;
        }
    }

    cout << "\n\nFactorisation terminee." << endl;
    getch();
    return 0;
}