main.cpp

00001 #include <cstdlib>
00002 #include <iostream>
00003 #include <vector>
00004 using namespace std;
00005 
00006 #include <math.h>
00007 #include "FiniteFieldElement.hpp"
00008 
00009 namespace Cryptography
00010 {    
00011         /*
00012             Elliptic Curve over a finite field of order P:
00013             y^2 mod P = x^3 + ax + b mod P
00014             
00015             NOTE: this implementation is simple and uses normal machine integers for its calculations.
00016                   No special "big integer" manipulation is supported so anything _big_ won't work. 
00017                   However, it is a complete finite field EC implementation in that it can be used 
00018                   to learn and understand the behaviour of these curves and in particular to experiment with them 
00019                   for their use in cryptography
00020                   
00021             Template parameter P is the order of the finite field Fp over which this curve is defined
00022         */
00023         template<int P>
00024         class   EllipticCurve
00025         {
00026             public:
00027                 // this curve is defined over the finite field (Galois field) Fp, this is the 
00028                 // typedef of elements in it
00029                 typedef FiniteFieldElement<P> ffe_t;
00030                 
00031                 /*
00032                     A point, or group element, on the EC, consisting of two elements of the field FP
00033                     Points can only created by the EC instance itself as they have to be 
00034                     elements of the group generated by the EC
00035                 */
00036                 class   Point
00037                 {
00038                     friend  class   EllipticCurve<P>;
00039                     typedef FiniteFieldElement<P> ffe_t;
00040                     ffe_t  x_;
00041                     ffe_t  y_;
00042                     EllipticCurve    *ec_;
00043                     
00044                     // core of the doubling multiplier algorithm (see below)
00045                     // multiplies acc by m as a series of "2*acc's"
00046                     void   addDouble(int m, Point& acc)
00047                     {        
00048                         if ( m > 0 )
00049                         {       
00050                             Point r = acc; 
00051                             for ( int n=0; n < m; ++n )
00052                             {
00053                                 r += r;     // doubling step                          
00054                             }
00055                             acc = r;
00056                         }        
00057                     }
00058                     // doubling multiplier algorithm
00059                     // multiplies a by k by expanding in multiplies by 2
00060                     // a is also an accumulator that stores the intermediate results
00061                     // between the "1s" of the binary form of the input scalar k
00062                     Point scalarMultiply(int k, const Point& a)
00063                     {
00064                         Point acc = a;
00065                         Point res = Point(0,0,*ec_);
00066                         int i = 0, j = 0;
00067                         int b = k;
00068                         
00069                         while( b )
00070                         {
00071                             if ( b & 1 )
00072                             {
00073                                 // bit is set; acc = 2^(i-j)*acc
00074                                 addDouble(i-j,acc);
00075                                 res += acc;           
00076                                 j = i;  // last bit set
00077                             }
00078                             b >>= 1;
00079                             ++i;
00080                         }
00081                         return res;
00082                     }
00083                     // adding two points on the curve
00084                     void    add(ffe_t x1, ffe_t y1, ffe_t x2, ffe_t y2, ffe_t & xR, ffe_t & yR) const
00085                     {
00086                          // special cases involving the additive identity                     
00087                         if ( x1 == 0 && y1 == 0 ) 
00088                         {
00089                             xR = x2;
00090                             yR = y2;
00091                             return;
00092                         }
00093                         if ( x2 == 0 && y2 == 0 )
00094                         {
00095                             xR = x1;
00096                             yR = y1;
00097                             return;
00098                         }
00099                         if ( y1 == -y2 ) 
00100                         {
00101                             xR = yR = 0;
00102                             return;
00103                         }
00104                                                 
00105                         // the additions
00106                         ffe_t s;                                                
00107                         if ( x1 == x2 && y1 == y2 )
00108                         {
00109                             //2P                          
00110                             s = (3*(x1.i()*x1.i()) + ec_->a()) / (2*y1);
00111                             xR = ((s*s) - 2*x1);                            
00112                         }
00113                         else
00114                         {   
00115                             //P+Q                            
00116                             s = (y1 - y2) / (x1 - x2);
00117                             xR = ((s*s) - x1 - x2);
00118                         }
00119                         
00120                         if ( s != 0 )
00121                         {
00122                             yR = (-y1 + s*(x1 - xR));
00123                         }
00124                         else
00125                         {
00126                             xR = yR = 0;
00127                         }           
00128                     }
00129                     
00130                     Point(int x, int y)
00131                     : x_(x),
00132                       y_(y),
00133                       ec_(0)
00134                     {}
00135                                   
00136                     Point(int x, int y, EllipticCurve<P> & EllipticCurve)
00137                      : x_(x),
00138                        y_(y),
00139                        ec_(&EllipticCurve)
00140                     {}
00141                     
00142                     Point(const ffe_t& x, const ffe_t& y, EllipticCurve<P> & EllipticCurve)
00143                      : x_(x),
00144                        y_(y),
00145                        ec_(&EllipticCurve)
00146                     {}
00147                     
00148                 public:                    
00149                     static  Point   ONE;
00150                     
00151                     // copy ctor
00152                     Point(const Point& rhs)
00153                     {
00154                         x_ = rhs.x_;
00155                         y_ = rhs.y_;
00156                         ec_ = rhs.ec_;
00157                     }
00158                     // assignment
00159                     Point& operator=(const Point& rhs)
00160                     {
00161                         x_ = rhs.x_;
00162                         y_ = rhs.y_;
00163                         ec_ = rhs.ec_;
00164                         return *this;
00165                     }
00166                     // access x component as element of Fp
00167                     ffe_t x() const { return x_; }
00168                     // access y component as element of Fp
00169                     ffe_t y() const { return y_; }
00170                     // calculate the order of this point by brute-force additions
00171                     // WARNING: this can be VERY slow if the period is long and might not even converge 
00172                     // so maxPeriod should probably be set to something sensible...
00173                     unsigned int     Order(unsigned int maxPeriod = ~0) const
00174                     {
00175                         Point r = *this;
00176                         unsigned int n = 0;
00177                         while( r.x_ != 0 && r.y_ != 0 )
00178                         {
00179                             ++n;
00180                             r += *this;
00181                             if ( n > maxPeriod ) break;
00182                         }
00183                         return n;
00184                     }
00185                     // negate
00186                     Point   operator-()
00187                     {
00188                         return Point(x_,-y_);
00189                     }                                        
00190                     // ==
00191                     friend bool    operator==(const Point& lhs, const Point& rhs)
00192                     {
00193                         return (lhs.ec_ == rhs.ec_) && (lhs.x_ == rhs.x_) && (lhs.y_ == rhs.y_);
00194                     }
00195                     // !=
00196                     friend bool    operator!=(const Point& lhs, const Point& rhs)
00197                     {
00198                         return (lhs.ec_ != rhs.ec_) || (lhs.x_ != rhs.x_) || (lhs.y_ != rhs.y_);
00199                     }                    
00200                     // a + b         
00201                     friend Point operator+(const Point& lhs, const Point& rhs)
00202                     {       
00203                         ffe_t xR, yR;
00204                         lhs.add(lhs.x_,lhs.y_,rhs.x_,rhs.y_,xR,yR);
00205                         return Point(xR,yR,*lhs.ec_);    
00206                     }
00207                     // a * int
00208                     friend  Point operator*(int k, const Point& rhs)
00209                     {
00210                         return Point(rhs).operator*=(k);
00211                     }
00212                     // +=
00213                     Point& operator+=(const Point& rhs)
00214                     {   
00215                         add(x_,y_,rhs.x_,rhs.y_,x_,y_);
00216                         return *this;  
00217                     }
00218                     // a *= int
00219                     Point& operator*=(int k)
00220                     {
00221                         return (*this = scalarMultiply(k,*this));
00222                     }                    
00223                     // ostream handler: print this point
00224                     friend ostream& operator <<(ostream& os, const Point& p)
00225                     {
00226                         return (os << "(" << p.x_ << ", " << p.y_ << ")");
00227                     }
00228                 };
00229                 
00230                 // ==================================================== EllipticCurve impl
00231                 
00232                 typedef EllipticCurve<P> this_t;
00233                 typedef class EllipticCurve<P>::Point point_t;
00234                 
00235                 // ctor
00236                 // Initialize EC as y^2 = x^3 + ax + b
00237                 EllipticCurve(int a, int b)
00238                 : a_(a),
00239                   b_(b),
00240                   m_table_(),
00241                   table_filled_(false)
00242                 {                    
00243                 }
00244                 // Calculate *all* the points (group elements) for this EC
00245                 //NOTE: if the order of this curve is large this could take some time...
00246                 void    CalculatePoints()
00247                 {
00248                     int x_val[P];
00249                     int y_val[P];
00250                     for ( int n = 0; n < P; ++n )
00251                     {
00252                         int nsq = n*n;
00253                         x_val[n] = ((n*nsq) + a_.i() * n + b_.i()) % P;
00254                         y_val[n] = nsq % P;                        
00255                     }
00256                     
00257                     for ( int n = 0; n < P; ++n )
00258                     {
00259                         for ( int m = 0; m < P; ++m )
00260                         {
00261                             if ( x_val[n] == y_val[m] )
00262                             {
00263                                 m_table_.push_back(Point(n,m,*this));
00264                             }
00265                         }
00266                     }
00267                     
00268                     table_filled_ = true;
00269                 }
00270                 // get a point (group element) on the curve 
00271                 Point   operator[](int n)
00272                 {
00273                     if ( !table_filled_ )
00274                     {
00275                         CalculatePoints();
00276                     }
00277                     
00278                     return m_table_[n];
00279                 }
00280                 // number of elements in this group
00281                 size_t  Size() const { return m_table_.size(); }
00282                 // the degree P of this EC
00283                 int     Degree() const { return P; }
00284                 // the parameter a (as an element of Fp)
00285                 FiniteFieldElement<P>  a() const { return a_; }
00286                 // the paramter b (as an element of Fp)
00287                 FiniteFieldElement<P>  b() const { return b_; }
00288                 
00289                 // ostream handler: print this curve in human readable form
00290                 template<int T>
00291                 friend ostream& operator <<(ostream& os, const EllipticCurve<T>& EllipticCurve);                       
00292                 // print all the elements of the EC group
00293                 ostream&    PrintTable(ostream &os, int columns=4);
00294                 
00295                 private:
00296                     typedef std::vector<Point>  m_table_t;
00297                     
00298                     m_table_t                   m_table_;   // table of points
00299                     FiniteFieldElement<P>       a_;         // paramter a of the EC equation
00300                     FiniteFieldElement<P>       b_;         // parameter b of the EC equation
00301                     bool    table_filled_;                  // true if the table has been calculated
00302         };
00303         
00304         template<int T>
00305             typename EllipticCurve<T>::Point EllipticCurve<T>::Point::ONE(0,0);
00306                                
00307         template<int T>
00308         ostream& operator <<(ostream& os, const EllipticCurve<T>& EllipticCurve)
00309         {
00310             os << "y^2 mod " << T << " = (x^3" << showpos;
00311             if ( EllipticCurve.a_ != 0 )
00312             {
00313                 os << EllipticCurve.a_ << "x";                
00314             }
00315             
00316             if ( EllipticCurve.b_.i() != 0 )
00317             {
00318                 os << EllipticCurve.b_; 
00319             }
00320             
00321             os << noshowpos << ") mod " << T;
00322             return os;
00323         }
00324         
00325         template<int P>
00326         ostream&    EllipticCurve<P>::PrintTable(ostream &os, int columns)
00327         {
00328             if ( table_filled_ )
00329             {
00330                 int col = 0;
00331                 typename EllipticCurve<P>::m_table_t::iterator iter = m_table_.begin();
00332                 for ( ; iter!=m_table_.end(); ++iter )
00333                 {
00334                     os << "(" << (*iter).x_.i() << ", " << (*iter).y_.i() << ") ";
00335                     if ( ++col > columns )
00336                     {
00337                         os << "\n";
00338                         col = 0;
00339                     }
00340                 }
00341             }
00342             else
00343             {
00344                 os << "EllipticCurve, F_" << P;
00345             }
00346             return os;
00347         }                        
00348 }
00349 
00350 namespace   utils
00351 {    
00352     float   frand()
00353     {
00354         static float norm = 1.0f / (float)RAND_MAX;
00355         return (float)rand()*norm;
00356     }
00357     
00358     int irand(int min, int max)
00359     {
00360         return min+(int)(frand()*(float)(max-min));
00361     }
00362 }
00363 
00364 using namespace Cryptography;
00365 using namespace utils;
00366 
00367 int main(int argc, char *argv[])
00368 {
00369     typedef EllipticCurve<263> ec_t;
00370     ec_t   myEllipticCurve(1,1);
00371     
00372    /* cout << "A little Elliptic Curve cryptography example\nby Jarl Ostensen, 2007\n\n";
00373     
00374     // print out a little info and test some properties
00375     cout << "The elliptic curve: " << myEllipticCurve << "\n";
00376     
00377     // calulate all the points for this curve. NOTE: in the real world this would not 
00378     // be a very sensible thing to do. If the period is very large this is big and slow
00379     myEllipticCurve.CalculatePoints();
00380     
00381     cout << "\nPoints on the curve (i.e. the group elements):\n";
00382     myEllipticCurve.PrintTable(cout,5);
00383     cout << "\n\n";
00384     
00385     ec_t::Point P = myEllipticCurve[2];
00386     cout << "some point P  = " << P << ", 2P = " << (P+P) << "\n";    
00387     ec_t::Point Q = myEllipticCurve[3];
00388     cout << "some point Q = " << Q << ", P+Q = " << (P+Q) << "\n"; 
00389     ec_t::Point R = P;
00390     R += Q;
00391     cout << "P += Q = " << R << "\n";
00392     R = P;
00393     R += R;
00394     cout << "P += P = 2P = " << R << "\n";
00395     
00396     cout << "\nEC message encryption example\n===============================================\n\n";
00397     */
00398     // Menes-Vanstone EC message encryption scheme
00399 
00400     // Public: the base point on the curve (i.e. base group element) used to generate keys
00401     // this is a REALLY slow way of picking a base point...
00402     ec_t::Point G = myEllipticCurve[0];
00403     while( (G.y() == 0 || G.x() == 0) || (G.Order()<2) )
00404     {
00405         int n = (int)(frand()*myEllipticCurve.Size());
00406         G = myEllipticCurve[n];
00407     }
00408     
00409     cout << "G = " << G << ", order(G) is " << G.Order() << "\n";
00410     
00411     // Alice
00412     int a = irand(1,myEllipticCurve.Degree()-1);
00413     ec_t::Point Pa = a*G;  // public key
00414     cout << "Alice' public key Pa = " << a << "*" << G << " = " << Pa << endl;    
00415         
00416     // Bob
00417     int b = irand(1,myEllipticCurve.Degree()-1);;
00418     ec_t::Point Pb = b*G;  // public key       
00419     cout << "Bob's public key Pb = " << b << "*" << G << " = " << Pb << endl;    
00420     
00421     // Jane, the eavesdropper
00422     int j = irand(1,myEllipticCurve.Degree()-1);;
00423     ec_t::Point Pj = j*G;
00424     cout << "Jane's public key Pj = " << j << "*" << G << " = " << Pj << endl;    
00425 
00426     cout << "\n\n";
00427     
00428     // Alice encrypts her message to send to Bob    
00429     // NOTE: the message first has to be split up into chunks that are in the Galois field F_p that is 
00430     // the domain of the EC
00431     // With P quite small (like in these examples) this is a serious limitation, but in the real world 
00432     // P could be very sizeable indeed, thus providing enough bits for good chunks
00433     int m1 = 19;
00434     int m2 = 72;
00435     
00436     cout << "Plain text message from Alice to Bob: (" << m1 << ", " << m2 << ")\n";
00437     
00438     // encrypt using Bob`s key
00439     ec_t::Point Pk = a*Pb;
00440     ec_t::ffe_t c1( m1*Pk.x() );
00441     ec_t::ffe_t c2( m2*Pk.y() );
00442     
00443     // encrypted message is: Pa,c1,c2
00444     cout << "Encrypted message from Alice to Bob = {Pa,c1,c2} = {" << Pa << ", " << c1 << ", " << c2 << "}\n\n";
00445     
00446     // Bob now decrypts Alice`s message, using her public key and his session integer "b" which was also used to generate his public key
00447     Pk = b*Pa;
00448     ec_t::ffe_t m1d = c1/Pk.x();
00449     ec_t::ffe_t m2d = c2/Pk.y();
00450     
00451     cout << "\tBob's decrypted message from Alice = (" << m1d << ", " << m2d << ")" << endl;
00452     
00453     // Jane intercepts the message and tries to decrypt it using her key
00454     Pk = j*Pa;
00455     m1d = c1/Pk.x();
00456     m2d = c2/Pk.y();
00457 
00458     cout << "\nJane's decrypted message from Alice = (" << m1d << ", " << m2d << ")" << endl;
00459     
00460     cout << endl;
00461     
00462     system("PAUSE");
00463     return EXIT_SUCCESS;
00464 }