// Supannika Rongsopa Spring 1996
// ecurve.C
// To compile:	g++ ecurve.C -o ecurve
// To encrypt:  ecurve datain dataout
// To decrypt:  ecurve dataout plain

#include<iostream.h>
#include<stdlib.h>
#include<stdio.h>
#include<Integer.h>
#include<fstream.h>
#include<ctype.h>
#include<time.h>
#include<math.h>

const Integer PRIME = 557;
const int TIMESTESTPRIME = 50;
Integer Tpoints;
struct pair                      
{  Integer x;
   Integer y;
};       
pair AB;                      // for storing a and b as in x^3 + ax + b
pair ALPHA;                   // for storing the generator
pair BETA;                    // for storing the multiples of the generator

char pairToLetter(pair n)
{  
char t;
if (n.x==1) { t = '\n'; return t; } if (n.x==3) { t = ' '; return t; }
if (n.x==5) { t = '!'; return t; }  if (n.x==7) { t = '"'; return t; }
if (n.x==9) { t = '#'; return t; }  if (n.x==11) { t = '$'; return t; }
if (n.x==13) { t = '%'; return t; }  if (n.x==15) { t = '&'; return t; }
if (n.x==17) { t = '\''; return t; } if (n.x==19) { t = '('; return t; }
if (n.x==21) { t = ')'; return t; }  if (n.x==23) { t = '*'; return t; }
if (n.x==25) { t = '+'; return t; }  if (n.x==27) { t = ','; return t; }
if (n.x==29) { t = '-'; return t; }  if (n.x==31) { t = '.'; return t; }
if (n.x==33) { t = '/'; return t; }  if (n.x==35) { t = '0'; return t; }
if (n.x==37) { t = '1'; return t; }  if (n.x==39) { t = '2'; return t; }
if (n.x==41) { t = '3'; return t; }  if (n.x==43) { t = '4'; return t; }
if (n.x==45) { t = '5'; return t; }  if (n.x==47) { t = '6'; return t; }
if (n.x==49) { t = '7'; return t; }  if (n.x==51) { t = '8'; return t; }
if (n.x==53) { t = '9'; return t; }  if (n.x==55) { t = ':'; return t; }
if (n.x==57) { t = ';'; return t; }  if (n.x==59) { t = '<'; return t; }
if (n.x==61) { t = '='; return t; }  if (n.x==63) { t = '>'; return t; }
if (n.x==65) { t = '?'; return t; }  if (n.x==67) { t = '@'; return t; }
if (n.x==69) { t = 'A'; return t; }  if (n.x==71) { t = 'B'; return t; }
if (n.x==73) { t = 'C'; return t; }  if (n.x==75) { t = 'D'; return t; }
if (n.x==77) { t = 'E'; return t; }  if (n.x==79) { t = 'F'; return t; }
if (n.x==81) { t = 'G'; return t; }  if (n.x==83) { t = 'H'; return t; }
if (n.x==85) { t = 'I'; return t; }  if (n.x==87) { t = 'J'; return t; }
if (n.x==89) { t = 'K'; return t; }  if (n.x==91) { t = 'L'; return t; }
if (n.x==93) { t = 'M'; return t; }  if (n.x==95) { t = 'N'; return t; }
if (n.x==97) { t = 'O'; return t; }  if (n.x==99) { t = 'P'; return t; }
if (n.x==101) { t = 'Q'; return t; }  if (n.x==103) { t = 'R'; return t; }
if (n.x==105) { t = 'S'; return t; }  if (n.x==107) { t = 'T'; return t; }
if (n.x==109) { t = 'U'; return t; }  if (n.x==111) { t = 'V'; return t; }
if (n.x==113) { t = 'W'; return t; }  if (n.x==115) { t = 'X'; return t; }
if (n.x==117) { t = 'Y'; return t; }  if (n.x==119) { t = 'Z'; return t; }
if (n.x==121) { t = '['; return t; }  if (n.x==123) { t = '\\'; return t; }
if (n.x==125) { t = ']'; return t; }  if (n.x==127) { t = '^'; return t; }
if (n.x==129) { t = '_'; return t; }  if (n.x==131) { t = '`'; return t; }
if (n.x==133) { t = 'a'; return t; }  if (n.x==135) { t = 'b'; return t; }
if (n.x==137) { t = 'c'; return t; }  if (n.x==139) { t = 'd'; return t; }
if (n.x==141) { t = 'e'; return t; }  if (n.x==143) { t = 'f'; return t; }
if (n.x==145) { t = 'g'; return t; }  if (n.x==147) { t = 'h'; return t; }
if (n.x==149) { t = 'i'; return t; }  if (n.x==151) { t = 'j'; return t; }
if (n.x==153) { t = 'k'; return t; }  if (n.x==155) { t = 'l'; return t; }
if (n.x==157) { t = 'm'; return t; }  if (n.x==159) { t = 'n'; return t; }
if (n.x==161) { t = 'o'; return t; }  if (n.x==163) { t = 'p'; return t; }
if (n.x==165) { t = 'q'; return t; }  if (n.x==167) { t = 'r'; return t; }
if (n.x==169) { t = 's'; return t; }  if (n.x==171) { t = 't'; return t; }
if (n.x==173) { t = 'u'; return t; }  if (n.x==175) { t = 'v'; return t; }
if (n.x==177) { t = 'w'; return t; }  if (n.x==179) { t = 'x'; return t; }
if (n.x==181) { t = 'y'; return t; }  if (n.x==183) { t = 'z'; return t; }
if (n.x==185) { t = '{'; return t; }  if (n.x==187) { t = '|'; return t; }
if (n.x==189) { t = '}'; return t; }  if (n.x==191) { t = '~'; return t; }
else t = '?'; return t; 
}  
pair charToNum(char x)
{  
pair n;
switch(x) { 
case '\n': n.x = 1; n.y = 2; break; case ' ': n.x = 3; n.y = 4; break;   
case '!':  n.x = 5; n.y = 6; break; case '"': n.x = 7; n.y = 8; break;    
case '#':  n.x = 9; n.y = 10; break; case '$': n.x = 11; n.y = 12; break;    
case '%':  n.x = 13; n.y = 14; break; case '&': n.x = 15; n.y = 16; break; 
case '\'': n.x = 17; n.y = 18; break; case '(': n.x = 19; n.y = 20; break; 
case ')':  n.x = 21; n.y = 22; break; case '*': n.x = 23; n.y = 24; break;
case '+':  n.x = 25; n.y = 26; break; case ',': n.x = 27; n.y = 28; break; 
case '-':  n.x = 29; n.y = 30; break; case '.': n.x = 31; n.y = 32; break;
case '/':  n.x = 33; n.y = 34; break; case '0': n.x = 35; n.y = 36; break;
case '1':  n.x = 37; n.y = 38; break; case '2': n.x = 39; n.y = 40; break; 
case '3':  n.x = 41; n.y = 42; break; case '4': n.x = 43; n.y = 44; break; 
case '5':  n.x = 45; n.y = 46; break; case '6': n.x = 47; n.y = 48; break; 
case '7':  n.x = 49; n.y = 50; break; case '8': n.x = 51; n.y = 52; break; 
case '9':  n.x = 53; n.y = 54; break; case ':': n.x = 55; n.y = 56; break; 
case ';':  n.x = 57; n.y = 58; break; case '<': n.x = 59; n.y = 60; break;
case '=':  n.x = 61; n.y = 62; break; case '>': n.x = 63; n.y = 64; break;
case '?':  n.x = 65; n.y = 66; break; case '@': n.x = 67; n.y = 68; break;
case 'A':  n.x = 69; n.y = 70; break; case 'B': n.x = 71; n.y = 72; break;
case 'C':  n.x = 73; n.y = 74; break; case 'D': n.x = 75; n.y = 76; break;
case 'E':  n.x = 77; n.y = 78; break; case 'F': n.x = 79; n.y = 80; break;
case 'G':  n.x = 81; n.y = 82; break; case 'H': n.x = 83; n.y = 84; break;
case 'I':  n.x = 85; n.y = 86; break; case 'J': n.x = 87; n.y = 88; break;
case 'K':  n.x = 89; n.y = 90; break; case 'L': n.x = 91; n.y = 92; break;
case 'M':  n.x = 93; n.y = 94; break; case 'N': n.x = 95; n.y = 96; break;
case 'O':  n.x = 97; n.y = 98; break; case 'P': n.x = 99; n.y = 100; break;
case 'Q':  n.x = 101; n.y = 102; break; case 'R': n.x = 103; n.y = 104; break;
case 'S':  n.x = 105; n.y = 106; break; case 'T': n.x = 107; n.y = 108; break;
case 'U':  n.x = 109; n.y = 110; break; case 'V': n.x = 111; n.y = 112; break;
case 'W':  n.x = 113; n.y = 114; break; case 'X': n.x = 115; n.y = 116; break;
case 'Y':  n.x = 117; n.y = 118; break; case 'Z': n.x = 119; n.y = 120; break;
case '[':  n.x = 121; n.y = 122; break; case '\\': n.x = 123; n.y = 124; break;
case ']':  n.x = 125; n.y = 126; break; case '^': n.x = 127; n.y = 128; break;
case '_':  n.x = 129; n.y = 130; break; case '`': n.x = 131; n.y = 132; break;
case 'a':  n.x = 133; n.y = 134; break; case 'b': n.x = 135; n.y = 136; break;
case 'c':  n.x = 137; n.y = 138; break; case 'd': n.x = 139; n.y = 140; break;
case 'e':  n.x = 141; n.y = 142; break; case 'f': n.x = 143; n.y = 144; break;
case 'g':  n.x = 145; n.y = 146; break; case 'h': n.x = 147; n.y = 148; break;
case 'i':  n.x = 149; n.y = 150; break; case 'j': n.x = 151; n.y = 152; break;
case 'k':  n.x = 153; n.y = 154; break; case 'l': n.x = 155; n.y = 156; break;
case 'm':  n.x = 157; n.y = 158; break; case 'n': n.x = 159; n.y = 160; break;
case 'o':  n.x = 161; n.y = 162; break; case 'p': n.x = 163; n.y = 164; break;
case 'q':  n.x = 165; n.y = 166; break; case 'r': n.x = 167; n.y = 168; break;
case 's':  n.x = 169; n.y = 170; break; case 't': n.x = 171; n.y = 172; break;
case 'u':  n.x = 173; n.y = 174; break; case 'v': n.x = 175; n.y = 176; break;
case 'w':  n.x = 177; n.y = 178; break; case 'x': n.x = 179; n.y = 180; break;
case 'y':  n.x = 181; n.y = 182; break; case 'z': n.x = 183; n.y = 184; break;
case '{':  n.x = 185; n.y = 186; break; case '|': n.x = 187; n.y = 188; break;
case '}':  n.x = 189; n.y = 190; break; case '~': n.x = 191; n.y = 192; break;
default :  n.x = 3; n.y = 4; break;  } 
return n;
} 
// Jacobi function 
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);
                  } 
                } 
        } 
}
// Converts the negative number in a mod system into its equivalent positive
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"; exit(1); }
 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;
}
// Computes the multiple of the given pair (a generator point on curve)
pair Multiples(Integer Px, Integer Py, Integer a)
{   // cout<<"\n "<<a<<"  *  "<<Px<<","<<Py;
    pair multiple;
    Integer Qx, Qy, lambda, QylessPy, QxlessPx, x3, y3, tempx3, tempy3;
    Integer tempPx = Px;
    Integer tempPy = Py;
    Qx = Px; 
    Qy = Py;
    Integer tempinv;
    for (Integer i = 2; i < (a+1); i++)  
      { if ((Px == Qx) && ((Py + Qy) % PRIME == 0))       // P + Q = infinity
	  { if (i == a) { cout<<"\nPoint at infinity! a = "<<a; exit(1); }
	    else { i = i + 1; 
		   Px = tempPx; 
		   Py = tempPy; }
	  }
        else
	{ if ((Px == Qx) && (Py == Qy))                        // if P == Q
	   {  tempinv = Inverse(PRIME,(2*Py)%PRIME);
              lambda = ( ( ((3*Px*Px)%PRIME) + AB.x) * tempinv) % PRIME;
	   }
          else                                                 // if P != Q
	   { if ((Qy-Py) < 0) QylessPy = makepositive((Qy-Py),PRIME);
		else QylessPy = (Qy-Py)%PRIME;
	     if ((Qx-Px) < 0) QxlessPx = makepositive((Qx-Px),PRIME);
                else QxlessPx = (Qx-Px)%PRIME;
             tempinv = Inverse(PRIME,(QxlessPx));
             lambda = ( (QylessPy) * tempinv ) % PRIME;
           }
         x3 = ((lambda*lambda) - Px - Qx) % PRIME; 
         y3 = ((lambda * (Px - x3)) - Py) % PRIME; 
         if (x3 < 0) { tempx3 = makepositive(x3,PRIME); x3 = tempx3; }
         if (y3 < 0) { tempy3 = makepositive(y3,PRIME); y3 = tempy3; }	
         Px = x3; Py = y3;
	}
      }
    multiple.x = x3; multiple.y = y3;
    return multiple;
}

char 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 'n';
  Integer a, jacobi, result;
  int count = 0;
  cout<<"\n\t\t\tPrime testing: "<<n;
  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 'n';
	       }
        }  
  cout<<" YES, "<<n<<" is prime.";
  return 'y';
}
// NEEDS SCHOOF'S ALGORITHM FOR THIS PROCEDURE !!!!!!!!
// Very slowly counts how many points there are on the elliptic curve
char CountPoints()
{ Integer temp, rightside;   // right side of equation y^2 = x^3 + ax + b
  Integer count = 0; 
  Integer C = (PRIME-1)/2;
  for (Integer i = 0; i < PRIME; i++)
   { rightside = ((i*i*i) + (AB.x*i) + (AB.y)) % PRIME;
     if (rightside == 0) count = count+1;
     else
     { temp = FastExp(rightside,C,PRIME); 
       if (temp == 1) 
       { count = count + 2;
         ALPHA.x = i; 
       }
     }
   }
  temp = count+1;      // each count is a pair, plus point o in infinity 
  Tpoints = temp;
  cout<<"Y^2 = x^3 + "<<AB.x<<"x + "<<AB.y;
  cout<<" has *****"<<Tpoints<<"***** points on it.";
  return 'y';
}
// Generates an elliptic curve equation
void getEquation()
{ Integer a, b, equation; 
  char OKpoints = 'n';
  while(OKpoints == 'n')
    { a = (Integer) (random() % PRIME);
      b = (Integer) (random() % PRIME);
      equation = ( (4*a*a*a)+(27*b*b) ) % PRIME;
      if (equation != 0) { cout<<"\n\nTesting the equation: ";
                           cout<<" y^2 = x^3 + "<<a<<"x + "<<b<<endl;
			   AB.x = a; AB.y = b;
			   OKpoints = CountPoints();
                         }
    }
}
// Computes y0 = k(ALPHA), (c1,c2) = k(BETA), and y1 = c1x1, y2 = c2x2
void encipherInput(char X, ofstream& OUTFILE)
{  pair y0, c, xinput;
   Integer  temp, k, y1, y2;
   while(1)
    { temp = (Integer) (random() % (Tpoints-2));           
      k = temp+1;
      y0 = Multiples(ALPHA.x, ALPHA.y, k);
      c = Multiples(BETA.x, BETA.y, k);
      if ((temp % Tpoints != 0) && (c.x * c.y != 0)) break;
    }
//   y0 = Multiples(ALPHA.x, ALPHA.y, k);          y0 = k * ALPHA
//   c = Multiples(BETA.x, BETA.y, k);           (c1,c2) = k * BETA
     cout<<"  "<<X;
   xinput = charToNum(X);
   y1 = (c.x * xinput.x) % PRIME;
   y2 = (c.y * xinput.y) % PRIME;
   OUTFILE<<y0.x<<" ";     OUTFILE<<y0.y<<" ";
   OUTFILE<<y1<<" ";       OUTFILE<<y2<<" "; 
}
void encipher(ifstream& infile, ofstream& outfile)
{   Integer secretA, temp, K;
    char letter, testresult = 'n';
    while(testresult == 'n')     // checks if total number of points is prime
    {          getEquation();
               testresult = PrimeTest(Tpoints);
    } 
    cout<<"\n\n\n----------------------------------------------------------";
    cout<<"\nTHE EQUATION IS FOUND: y^2 = x^3 + "<<AB.x<<"x + "<<AB.y; 
    cout<<"   (mod "<<PRIME<<")"; 
    
    cout<<"\n\nALPHA is ("<<ALPHA.x;
    temp = ((ALPHA.x * ALPHA.x * ALPHA.x) + (AB.x * ALPHA.x) + AB.y) % PRIME;
    for (Integer j = 1; j < PRIME; j++)
         if ((j*j) % PRIME == temp) ALPHA.y = j; 
    cout<<","<<ALPHA.y<<")";
    while(1)
    { temp = (Integer) (random() % (Tpoints-2));     // secret exponent  
      secretA = temp+1;
      if (temp % Tpoints != 0) break;
    }
    cout<<"\n\nSECRET NUMBER IS NEEDED FOR DECIPHERING: ***** ";
    cout<<secretA<<" *****"; 
    cout<<"\n----------------------------------------------------------\n\n";
    BETA = Multiples(ALPHA.x, ALPHA.y, secretA);      // BETA=secret*ALPHA
    outfile<<AB.x<<" ";        // writes 'a' of x^3 + ax + b into output file
    while(infile)
     { infile.get(letter);
       //  cout<<"\nread in: "<<letter;
       encipherInput(letter, outfile); 
     }
    outfile.close();
}
void decipherInput(Integer aa, Integer Y0x, Integer Y0y, 
                               Integer Y1, Integer Y2, ofstream& FILEOUT)
{       Integer c1inverse, c2inverse;
        pair numpair; char letter; 
	pair c = Multiples(Y0x,Y0y,aa); 
	c1inverse = Inverse(PRIME,c.x);  
	c2inverse = Inverse(PRIME,c.y); 
	numpair.x = (Y1 * c1inverse) % PRIME; 
        numpair.y = (Y2 * c2inverse) % PRIME; 
	letter = pairToLetter(numpair); 
	cout<<letter;
        FILEOUT<<letter;
} 
void decipher(ifstream& infile, ofstream& outfile)
{   Integer abx, a, y0x, y0y, y1, y2;
    cout<<"\n\n\t\tPlease enter your secret number a: "; cin>>a;
    infile>>abx; AB.x = abx;
    while(infile)
      { infile>>y0x; 
        infile>>y0y; 
	infile>>y1;  
	infile>>y2; 
	decipherInput(a,y0x,y0y,y1,y2,outfile); 
      }
    outfile.close();
    cout<<"\n\n";    
}  
void main(int argc, char** argv)
{      if (argc != 3) 
	 { cout<<"\nUsage: ecurve filein fileout\n"; exit(1); }
       int choice;
       cout<<"\n              --------------------------------------------";	
       cout<<"\n              MENEZES-VANSTONE ELLIPTIC CURVE CRYPTOSYSTEM";
       cout<<"\n              --------------------------------------------\n";
       cout<<"\n\t\tMENU: 1. encipher";
       cout<<"\n\t\t      2. decipher";
       cout<<"\n\t\t      3. exit";
       cout<<"\n\n\t\tMake your selection now, 1, 2 or 3: "; cin>>choice;
       ifstream filein(argv[1]); 
       ofstream fileout(argv[2]);
       srandom(time(NULL));	
       switch(choice) { case 1: encipher(filein, fileout); break;
		        case 2: decipher(filein, fileout); break;
			case 3: exit(1);
			default: cout<<"\nIllegal selection\n"; exit(1);
		      }
       cout<<endl;
}