//Supannika Rongsopa 
// myfactor.C 
#include<stdio.h>
#include<Integer.h>
#include<fstream.h>
#include<ctype.h>
#include<time.h>
#include<math.h>
#include<stdlib.h>
#include<new.h>
#include<iostream.h>  

const int MAXANSWERLENGTH = 1000;
const int TIMESTESTPRIME = 50;
const int B = 8;
const int TIMES = 700000;
struct pair
{ Integer x;
  Integer y;
};
pair Point;
Integer A, N, K = 1;

//------------------------- *** ARRAY CLASS *** --------------------------    
class Array                                       
{ private: Integer * dataptr;		              //pointer to element 
           unsigned long size;                    //size of Array instance
  public:  
  Array(unsigned long);                                        	       
 ~Array() { delete [] dataptr; }                                
  Array & operator = (const Array & Value);       
  Integer & operator [] (unsigned long) const; 
  void assignNum(unsigned long, Integer);    
  unsigned long getsize() const;		  
  friend ostream & operator << (ostream & Out, const Array & Value );    
};                               
   //Constructor for size n Array.          
   Array::Array(unsigned long n)        
   { if (n<=0) { cerr<<"1.illegal array size: "<<n<<endl;
		 exit(1); }
     size = n;
     dataptr = new Integer[size];
     for (unsigned long i = 0; i < size; i++)
        dataptr[i] = i;   
   }   
   //Overloaded subscript operator. It returns the value stored in cell i.
   Integer & Array::operator [] (unsigned long i) const       
   { if ( i < 0 || i > size+1) { cerr<<"3.illegal array size: "<< i <<endl;
			         exit(1); }
     return (dataptr[i]);            
   }                
   //Overloaded assignment operator. It first checks if its only argument 
   //is itself. If not, it proceeds to delete itself (lhs) and creates a new 
   //Array instance with same size and copies every element. It returns
   //itself (*this - lhs). 
   Array & Array::operator = (const Array & Value)    
   { if (this != &Value)          //if not, deletes lhs and creates new Array  
       { delete [] dataptr;        //with elements copied one by one from rhs
         size = Value.size;
         dataptr = new Integer[size];
         for (unsigned long i = 0; i < size; i++)
	    dataptr[i] = Value.dataptr[i]; 
         return *this;                    
       }
   }
   //Member function for assigning value to a cell [n].
   void Array::assignNum(unsigned long n, Integer value)
   { if ( n < 0 || n > size+1 ) {cerr<<"4.illegal array subscript:"<<n<<endl;
                                 exit(1); }
     dataptr[n] = value; 
   }
   //The operator << is overloaded to print the values stored in the Array.
   //This function is a friend function. Both arguments are sent as reference
   //(we don't want to make a copy of ostream). The second argument is sent 
   //as constant so it can't be changed. It returns ostream &.
   ostream & operator << (ostream & Out, const Array & Value)
   { for (unsigned long i=0; i < Value.getsize(); i++)
       Out <<Value[i]<<"  ";
     return Out;
   }
   //Accessor for Array's size.                 
   unsigned long Array::getsize() const 
   { return size; } 
// ---------------------------------------------------------------------------

Array kValues(B);

Integer makepositive(Integer n, const Integer &P)
{ Integer temp = 0;
  while ( (temp + n) < 0 ) 
       temp = temp + P;
  n = temp + n; 
  return n;
}
// Computes the inverse of 'num' in the modulo 'm' system
Integer Inverse (const Integer &m, const Integer &num)
{ if (num==0) 
    { cout<<"\n0 has no inverse in a mod "<<m<<" system.\n\n"; exit(1); }
  Integer n0, inv, b0, t, t0, q, r, temp, tempnum;
  if (num < 0) tempnum = makepositive(num,m);
    else tempnum = num; 
  n0 = m;               b0 = tempnum;                  t0 = 0;
  t = 1;                q = (Integer) (n0/b0);
  r = n0 - q * b0;
  while (r > 0)
       { temp = t0 - q * t;
         if (temp >= 0)       temp = temp % m;
         else if (temp < 0)   temp = m - ((-temp) % m);
         t0 = t;              t = temp; n0 = b0;
         b0 = r;              q = (Integer) (n0/b0);
         r = n0 - q * b0;
       }
 if (b0 != 1) 
  { cout<<"\n\n"<<"\t"<<num<<" has no inverse modulo "<<m<<"\n\n"; return 0; }
 else inv = t % m;  
 return inv;
}
// Computes num^exp mod n
Integer FastExp (const Integer &num, const Integer &exp, const Integer &n)
{       Integer c1, exp1, x;
 	c1 = num; 	exp1 = exp;        x = 1;
 	while (exp1 != 0)
   	{    while ((exp1 % 2) == 0)
        	{ exp1 = exp1/2;
          	  c1 = (c1 * c1) % n;
        	}
             exp1 = exp1 - 1;
             x = (x * c1) % n;
   	}
 	return x;
}
Integer Jacobi (const Integer & rand_num, const Integer & num)
{	Integer d, p1, temp;
 	if (rand_num == 1)      return 1;
 	d = gcd(rand_num, num);
 	if (d > 1)              return 0;
 	else
   	{	if (rand_num % 2== 0)
      		{ temp = (num*num-1)/8;   
                  if (temp % 2 == 0)
	 	  { p1 = rand_num/2;
          	    return Jacobi(p1,num);         	  }
                  else
	          { p1 = rand_num/2;
                    return -1*Jacobi(p1,num);             }
                }
                else
                { temp = (rand_num - 1)*(num - 1)/4;
       	          if (temp % 2 == 0)
	          { p1 = num % rand_num ;
                    return Jacobi(p1,rand_num);           }
                  else
                  { p1 = num % rand_num ;  
                    return -1*Jacobi(p1,rand_num);        } 
                } 
        } 
}
void foundfactor(Integer denominator)
{   Integer factor = gcd(denominator, N);
    cout<<"\n*** One of the factors of "<<N<<" is: "<<factor<<" ***\n\n";
}
// Computes the multiple given a pair - pair^mult mod PRIME.
pair Multiples (const pair &p, const Integer &mult)
{       cout<<"\nComputing "<<mult<<" * "<<p.x<<","<<p.y;
        pair answer, P; 
        P.x = p.x; P.y = p.y; answer.x = p.x; answer.y = p.y;
        Integer tempinv, lambda, tempx3, tempy3, M = mult-1;
	Integer QylessPy, QxlessPx, x3, y3;
	while (M != 0)
   	{   
	     while ((M % 2) == 0)
             { M = M/2;         // c1 = (c1 * c1) % n;   
	       tempinv = Inverse(N, (2*P.y) % N); 
	       if (tempinv == 0) { foundfactor(2*P.y); exit(1); }
               lambda = ((((3*P.x*P.x) % N) + A) * tempinv) % N;

	       x3 = ((lambda*lambda) - P.x - P.x) % N; 
               y3 = ((lambda * (P.x - x3)) - P.y) % N; 
               if (x3 < 0) { tempx3 = makepositive(x3, N); x3 = tempx3; }
               if (y3 < 0) { tempy3 = makepositive(y3, N); y3 = tempy3; }
	       P.x = x3; P.y = y3; 
             }
  
             M = M - 1;            // x = (x * c1) % n; 
             if (P.x == answer.x && P.y == answer.y)
	     { tempinv = Inverse(N, (2*P.y) % N); 
  	       if (tempinv == 0) { foundfactor(2*P.y); exit(1); }
               lambda = ((((3*P.x*P.x) % N) + A) * tempinv) % N;
	     }
	     else
	     { if ((answer.y-P.y) < 0) 
	        QylessPy = makepositive((answer.y-P.y), N);
               else QylessPy = (answer.y-P.y) % N;
	       if ((answer.x-P.x) < 0) 
	        QxlessPx = makepositive((answer.x-P.x), N);
               else QxlessPx = (answer.x-P.x) % N;
               tempinv = Inverse(N, (QxlessPx));
               if (tempinv == 0) { foundfactor(QxlessPx); exit(1); }
	       lambda = ( (QylessPy) * tempinv ) % N; 
	     }
	     x3 = ((lambda*lambda) - P.x - answer.x) % N; 
             y3 = ((lambda * (P.x - x3)) - P.y) % N; 
             if (x3 < 0) { tempx3 = makepositive(x3, N); x3 = tempx3; }
             if (y3 < 0) { tempy3 = makepositive(y3, N); y3 = tempy3; }
	     answer.x = x3; answer.y = y3;
              cout << x3 << "  " << y3 << endl;
	}
	answer.x = x3; answer.y = y3;
        
 	return answer;
}
int PrimeTest(Integer n)
{ /* if (n % 2 == 0 || n % 3 == 0 || n % 5 == 0 || n % 7 == 0 || n % 11 == 0 ||
      n % 13 == 0 || n % 17 == 0 || n % 19 == 0 || n % 23 == 0 || n % 29 == 0)
     return 1;  */ 
  Integer a, jacobi, result;
  int count = 0;
  while(count < TIMESTESTPRIME)
        { a = (Integer) (random() % (n-1));             
	  if (gcd(a, n) == 1)
               { count++;
		 jacobi = Jacobi(a, n);
                 result = FastExp(a, (n - 1)/2, n);
		 if (result == n-1) result = -1;
                 if (jacobi != result)  
		   return 0;         // not prime
	       }
        }  
  return 1;
}
void getEquation()
{   srandom(time(NULL));
    Integer a, b, x, y, equation, temp;
    while(1)
    { a = (Integer) (random() % N);
      x = (Integer) (random() % N);
      y = (Integer) (random() % N);
      temp = ((y*y) - (x*x*x) - (a*x)) % N; 
      b = makepositive(temp,N);
      equation = ((4*a*a*a)+(27*b*b)) % N;
      if (equation != 0) 
        { Point.x = x; Point.y = y; A = a;
          cout<<"\n\nThe equation is: y^2 = x^3 + "<<a<<"x + "<<b;
          break;
        }
     }
}
void getK()
{   // Array kValues(B);    
    cout<<"\nFactoring number N = "<<N;
    Integer C = (sqrt(N)+1) + ( 2*(sqrt(sqrt(N)+1)+1));
    cout<<"\nC is "<<C;
    for (unsigned long i = 2; i < B; i++)
       if (i == 2 || i == 3 || i == 5 || i == 7)
	  {   Integer temp = kValues[i] * i;
	      while(temp <= C)
	      { Integer value = kValues[i] * i;
	        kValues.assignNum(i, value);
	        temp = kValues[i] * i;
	      }
          }
    cout<<"\nThe array contains: "<<kValues;
    for (i = 2; i < B; i++)
       if (kValues[i] != i)   K = K*kValues[i];    
    cout<<"\nK is "<<K;
  //pair temp;
  //temp = Multiples(Point, K);    cout<<"\nK*Point is "<<temp.x<<","<<temp.y;
}
void factor()
{   char answer[MAXANSWERLENGTH]; int tries = 1;
    cout<<"\nEnter your number to be factored: "; cin>>answer;
    N = atoI(answer);
    while (N <= 0)
      { cout<<"\n\tNegative integer is not allowed.";
        cout<<"\n\tEnter a POSITIVE, NON-ZERO integer: "; cin>>answer;
	N = atoI(answer);
      }
    if (N == 1 || N == 2 || N == 3)
      { cout<<"\n\t"<<N<<" is prime."; return; }
    else if (PrimeTest(N) == 1)
      { cout<<"\n\t"<<N<<" is prime."; return; }
    else 
      { getK();
        while(tries < TIMES)
        { getEquation();
          pair temp;
	  temp = Multiples(Point, K); 
	  cout<<"\nK*Point is "<<temp.x<<","<<temp.y;
	  tries++;
        }
        cout<<"\n*** Cannot find a factor. "<<N<<" could be prime. ***";
      }
}    
void main()
{      int selection;     
       while(1)
	 { cout<<"\n\n\t--------------------------------";
           cout<<"\n\tELLIPTIC CURVE FACTORING PROGRAM";
           cout<<"\n\t--------------------------------\n";
           cout<<"\n\t1. Enter a number to be factored.";
           cout<<"\n\t2. Exit.";
           cout<<"\n\n\tMake your selection now: "; cin>>selection;
	   switch(selection) { case 1: factor(); break;
			       case 2: exit(1);
			       default: cout<<"\nIllegal entry\n"; exit(1);
		             }
	 }       
       cout<<endl;
}