ecc.cpp

椭圆曲线密码系统教程和相关源代码,是非常难得的学术资料

// --------------------------------------------------------------------
// ECDLP solver using Pohlig-Hellman algorithm to reduce problem
//       size and Pollard's rho algorithm to solve ECDLP over sub
//       -group (+ a routine to brute-force group with small order
//       to avoid pb with order 2...)
//
// 
//   Pollard's rho part is a (very) modified version of 'Miracl/index.c'(DLP)
//   You will be able to find more info reading HoAC or any others
//   good crypto-papers... :) 
//                                   Amenesia//TKM! 
// --------------------------------------------------------------------
// Info: You have to use Miracl to be able to compile it...
// ATM parameters permit to solve the ECDLP used in pDriLl's keygenme #5
// but it's quite easy to modify it in order to solve your own ECDLP
// (is well commented in a bad english :p)
// --------------------------------------------------------------------
                  
#include <stdlib.h>
#include <miracl.h>

#define NPRIMES 10
static miracl *mip;
static big order, rholim;



// --------------------------- Pollard rho ---------------------------

void iterate(epoint* X,epoint* Q,epoint* R,big a,big b)
{ // Random walk...
//       = a(i)     if X in S1
// a(i+1) = a(i) + 1   if X in S2 
//       = 2*a(i)   if X in S3
//
//       = b(i)     if X in S2
// b(i+1) = b(i) + 1   if X in S1 
//       = 2*b(i)   if X in S3
//
//   X(i)   = a(i)*Q + b(i)*R
//   X(i+1) = f(X(i))
//
// BTW: the criteria used by Dimedrol (ecdlp.solver) is quite
//     good (simple and fast to compute) so i take the same :)

  big p = mirvar(0);   
  epoint_get(X,p,p);

  if (remain(p,7)>4)
  { 
    ecurve_add(Q,X);
    incr(a,1,a);
    if (compare(a,order)==0) zero(a);
    mirkill(p);     
    return;
  }
  if (remain(p,7)>2) 
  {   
    ecurve_add(X,X);
    premult(a,2,a);
    if (compare(a,order)>=0) subtract(a,order,a);
    premult(b,2,b);
    if (compare(b,order)>=0) subtract(b,order,b);
    mirkill(p);           
    return;
  }     
  ecurve_add(R,X);
  incr(b,1,b);
  if (compare(b,order)==0) zero(b);
  mirkill(p);   
}



// To avoid endless loops...
void ecdlpbf(epoint* Q,epoint* R,big k)
{
  epoint* T = epoint_init();
  copy(order,k);
  negify(k,k);
  do
  {
    incr(k,1,k);
    ecurve_mult(k,Q,T);     
  } while (!epoint_comp(T,R));
  epoint_free(T);
}


void rho(epoint* Q,epoint* R,big k)
{ // Find ax,ay,bx and by with: ax*Q + bx*R == ay*Q + by*R
// So : (ax-ay)*R = (by-bx)*Q
// ECDLP => k exists such as k*R = Q
// So: (ax-ay) = (by-bx)*k mod order
//   k = (ax-ay)*(by-bx)^1 mod order
// BTW: check if (by-bx) is inversible
//     (order is prime so (by-bx) is 
//     inversible <=> (by-bx) != 0)

  long rr,i;
  big ax,bx,ay,by,m,n;
  epoint* Y = epoint_init();
  epoint* X = epoint_init();
  m=mirvar(0);
  n=mirvar(0);   
  ax=mirvar(0);
  bx=mirvar(0);
  ay=mirvar(0);
  by=mirvar(2);
  
  if (compare(rholim,order)>=0) 
  {
    ecdlpbf(Q,R,k);
  }
  else 
  { 
    do
        {
          rr=1L;   
          bigrand(order,ay);     
          bigrand(order,by);     
          ecurve_mult2(ay,Q,by,R,Y);
          do
          { // Search a collision...
                    epoint_copy(Y,X);
                    copy(ay,ax);
                    copy(by,bx);
                    rr*=2;
                    for (i=1;i<=rr;i++)
                    {
                      iterate(Y,Q,R,ay,by);
                      if (epoint_comp(X,Y)) break;
                    }
          } while (!epoint_comp(X,Y));
        } while (compare(bx,by)==0);

        subtract(ay,ax,m);
        if(size(m)<0){add(m,order,m);}
        subtract(bx,by,n); 
        if(size(m)<0){add(m,order,m);}           
        xgcd(n,order,n,n,n); 
        mad(m,n,n,order,order,k); 

// I don't know why but it seems that
//   -k*A != (-k+ord(A))*A 
// If you are able to explain me why
// feel free to contact me... :)
        
ecurve_mult(k,Q,X);     
if (!epoint_comp(X,R)){ subtract(k,order,k); }                                                                                                                                                                                                                                                                                                                                                     
ecurve_mult(k,Q,X); 
if (!epoint_comp(X,R)){ printf("Error !!!"); } 
  }   
epoint_free(Y);
  epoint_free(X);
  mirkill(by);
  mirkill(ay);
  mirkill(bx);
  mirkill(ax);
}

// --------------------------- End Pollard rho ---------------------------





int main()
{
  mip =mirsys(100,0);
  unsigned char mod_str[10];
  int i,j;
  int pw[NPRIMES];   
  big pf[NPRIMES],logmod[NPRIMES];
  for (i=0;i<NPRIMES;i++) 
  {
    pf[i]=mirvar(0);
    logmod[i]=mirvar(0);
  }   
  
  big_chinese bc;   
  
  big k=mirvar(0);
  big p=mirvar(0);
  big a=mirvar(0);
  big b=mirvar(0);
  big p1=mirvar(0);     
  big xA=mirvar(0);
  big xB=mirvar(0);

  epoint* A = epoint_init();   
  epoint* At = epoint_init(); 
  epoint* B = epoint_init(); 
  epoint* Q= epoint_init(); 
  epoint* R= epoint_init(); 
  epoint* Rj= epoint_init(); 

  big c=mirvar(0);   
  big N=mirvar(0);
  order=mirvar(0);
  rholim=mirvar(2);
  irand(41563436);   
  
  mip->IOBASE=16; 
  // --------------------- Elliptic Curve Definition --------------------- 
  printf("In progress...\n");
                  
  cinstr(a,"1");   
  cinstr(b,"0");   
  cinstr(p,"ACC00CF0775153B19E037CE879D332BB");   
    
  cinstr(xA,"18650A0922FC3EC0B778347B20EE6619"); 
  cinstr(xB,"0D85FA1C5BC8982485ACD9FA9B1BEBEC");     
  

  ecurve_init(a,b,p,MR_AFFINE);
  epoint_set(xA,xA,0,A);
  epoint_set(xB,xB,1,B);


  // --------------------- Order factors ---------------------   
  mip->IOBASE=16;
  cinstr(p1,"566006783BA8A9D8CF01BE743CE9995E");                 
  cinstr(pf[0],"2"); pw[0]=1;
  cinstr(pf[1],"3"); pw[1]=1;
  cinstr(pf[2],"7"); pw[2]=1;
  cinstr(pf[3],"D"); pw[3]=2;     
  cinstr(pf[4],"7F"); pw[4]=1;
  cinstr(pf[5],"D3"); pw[5]=1;
  cinstr(pf[6],"1DF75"); pw[6]=1;
  cinstr(pf[7],"5978F"); pw[7]=1;
  cinstr(pf[8],"1F374C47"); pw[8]=1;   
  cinstr(pf[9],"5F73FD8D3"); pw[9]=1;   
            
  
  // --------------------- Pohlig Hellman --------------------- 
  // If ord(A) = p1^e1 * p2^e2 * ... * pn^en there is a group 
  // isomorphism such:   f : <A> -> Cp1^e1 * ... * Cpn^en 
  // where Cpi^ei are subgroup of <A> of order pi^ei.
  //
  // The projection of f to Cpi^ei is given by:
  //       fi : <A> -> Fpi^ei
  //           R -> (N/pi^ei)*R
  //   Moreover fi is a group homomorphism so because of linearity
  //   ("lin閍rit? en francais mais j'ai quelques doutes sur son
  //   equivalent en anglais): B = k*A so fi(B) = fi(k*A) = k*fi(A)
  //   but we are in Cpi^ei so k is determined modulo (pi^ei). 
  //
  //   Using CRT we will find p mod (p1^e1* ... * pn^en)...
  //   If you are really interested in PH algo you sould read more :)

  for (i=0;i<NPRIMES;i++)
  {       
    // ui = 0
    zero(logmod[i]); 
    // Q = (n/pi)*A               
    copy(p1,N);     
    divide(N,pf[i],N);
    ecurve_mult(N,A,Q);
    // R(0) = B
    epoint_copy(B,Rj);
                
    // For j=0 to (e-1)
    for (j=0;j<pw[i];j++)
    {           
          // R = (n/pi^(j+1))*Rj         
          ecurve_mult(N,Rj,R);
          // Solve R = kj*Q   
          copy(pf[i],order);                                                 
          rho(Q,R,k); 
          // c = kj*pi^j     
          powltr(j,pf[i],p1,c);
          power(pf[i],j,p1,c);                     
          multiply(k,c,c);
          // Rj+1 = Rj - kj*pi^j*A
          ecurve_mult(c,A,At);
          ecurve_sub(At,Rj);           
          // ui = ui + kj*pi^j       
          add(logmod[i],c,logmod[i]);
          // N = (n/pi^(j+2))           
          divide(N,pf[i],N);                   
    }
    power(pf[i],pw[i],p1,pf[i]); 
    // Interface
    cotstr(pf[i],mod_str);           
    printf("# Solved over C(%s)\n#> ",mod_str); 
    cotnum(logmod[i],stdout);     
  }           
  
  // --------- Chinese remainder thereom ---------     
  crt_init(&bc,NPRIMES,pf);
  crt(&bc,logmod,k);   
  crt_end(&bc);

  // -------------------- User-friendly Interface :) -------------------- 
      
  ecurve_mult(k,A,Q);         
  if (epoint_comp(Q,B))
  {
    printf("\nDiscrete logarithm: "); 
    cotnum(k,stdout); 
  }
  else
  { 
    printf("Wrong solution !? :x");
  }           
  getch();


  // -----------------------------------------------------------------
  epoint_free(A);
  epoint_free(At);   
  epoint_free(B); 
  epoint_free(Q);
  epoint_free(R); 
  epoint_free(Rj);
  
  for (i=0;i<NPRIMES;i++) 
  {
    mirkill(pf[i]);
    mirkill(logmod[i]);
  }     
  mirkill(k);
  mirkill(p);
  mirkill(a);
  mirkill(b);
  mirkill(p1);     
  mirkill(xA);
  mirkill(xB);
  mirkill(c);
  mirkill(N);
  mirexit(); 
  return(0); 
}