ncurve.cpp

ecppki-cpp curve

// elliptic curve implementation

#include "vlong.hpp"
#include "ncurve.hpp"

#define assert_on        0      // Consistency checking (for debugging)
#define use_field_op     0      // Use operators

#if assert_on
  #include "prime.hpp"
  static void assert0(int x){ if (x==0) while (1/x); }
  #define assert(x) assert0(x)
#else
  #define assert(x)
#endif

small_field::small_field( unsigned L0, unsigned root )
: L(L0),
  BASE( 1u<<L0 ),
  BASE_M1( (1u<<L0)-1 ),
  log( new lunit[ 1u<<L0 ] ),
  alog( new lunit[ 1u<<L0 ] )
{
  assert( L0 <= sizeof(lunit)*8 );
  root |= BASE;
  unsigned j = 1;
  for (unsigned i = 0; i < BASE_M1; i+=1)
  {
    log[j] = i;
    alog[i] = j;
    j *= 2; // 2 is primitive element
    if (j & BASE) j ^= root;
  }
  log[0] = BASE_M1; // used by field::mul
}

small_field::~small_field()
{
  delete [] log;
  delete [] alog;
}

void field::add( const poly a, const poly b, poly c )
{
  if ( a[0] < b[0] )
    add(b,a,c);
  else
  {
    unsigned i;
    for (i=0;i<=b[0];i+=1)
      c[i] = a[i] ^ b[i];
    for (;i<=a[0];i+=1) c[i] = a[i];
    c[0] = a[0];
  }
  while ( c[0] && c[c[0]]==0 ) c[0] -= 1;
}

void field::copy ( const poly a, poly b )
{
  for (unsigned i=0;i<=a[0];i+=1)
    b[i] = a[i];
}

int field::equal( const poly a, const poly b )
{
  unsigned i = 0;
  while (1)
  {
    if (a[i] != b[i]) return 0;
    if ( i == a[0] ) return 1;
    i += 1;
  }
}

void field::set_random( poly a )
{
  a[0] = K;
  for (unsigned i=1;i<=K;i+=1)
    a[i] = rand(BASE);
  while ( a[0] && a[a[0]]==0 )
    a[0] -= 1;
}

void field::reduce( poly a )
{
  for (unsigned i=a[0];i>K;i-=1)
  {
    a[i-K] = a[i-K] ^ a[i];
    a[i+T-K] = a[i+T-K] ^ a[i];
    a[i] = 0;
  }
  if ( a[0] > K )
  {
    a[0] = K;
    while (a[0] && a[a[0]]==0) a[0] -= 1;
  }
}

field_element& field_element:: operator = ( const field_element & x )
{
  f = x.f;
  field::copy( x.v, v );
  return *this;
}

field_element::field_element()
{
  f=0;
  v[0]=0;
}

field_element::field_element( field * F )
{
  f=F;
  v[0]=0;
}

field_element::field_element( const field_element & x )
{
  f = x.f;
  field::copy( x.v, v );
}

int field_element::operator == ( const field_element & x ) const
{
  return field::equal( v, x.v );
}

int field_element::operator == ( unsigned x ) const
{
  return ( v[0] == 0 && x == 0 ) || ( v[1] == 1 && x == v[1] );
}

field_element field_element::operator + ( const field_element & x ) const
{
  field_element result(f);
  field::add( v, x.v, result.v );
  return result;
}

field_element field_element::operator * ( const field_element & x ) const
{
  field_element result(f);
  f->mul( v, x.v, result.v );
  f->reduce( result.v );
  return result;
}

field_element curve::sq( const field_element & x )
{
  field_element result(x.f);
  x.f->square( x.v, result.v );
  x.f->reduce( result.v );
  return result;
}

field_element field_element::operator / ( const field_element & x ) const
{
  field_element result(f);
  field::poly tmp;
  f->inverse( x.v, tmp );
  f->mul( v, tmp, result.v );
  f->reduce( result.v );
  return result;
}

field_element sqrt( const field_element x )
{
  field_element result(x.f);
  x.f->sqrt( x.v, result.v );
  return result;
}

unsigned curve::ybit( const field_element & x )
{
  if ( x.v[0]==0 ) return 0;
  return x.v[1] & 1;
}

int field::slow_trace( const poly a )
{
  field::poly t1;
  copy( a, t1 );
  for (unsigned i=1;i<M;i+=1)
  {
    field::poly t2;
    square( t1, t2 );
    reduce( t2 );
    add( t2, a, t1 );
  }
  return t1[0] != 0;
}

int field::trace( const poly a )
{
  unsigned min = a[0]; if ( min > tm[0] ) min = tm[0];
  unsigned w = 0;
  for (unsigned i=1;i<=min;i+=1)
    w ^= a[i] & tm[i];
  // compute parity of w
  unsigned result = 0;
  while (w)
  {
    result ^= w & 1;
    w >>= 1;
  }
  assert( result == slow_trace(a) );
  return result;
}

void field::quad_solve( const poly a, poly b )
{
  if (M % 2 == 1) // M is odd - compute half-trace
  {
    copy( a, b );
    for (unsigned i=0;i<M/2;i+=1)
    {
      field::poly t1,t2;
      square( b, t1 ); reduce( t1 );
      square( t1, t2 ); reduce( t2 );
      add( t2, a, b );
    }
  }
  else
  {
    field::poly d;
    b[0] = 0;
    copy( a, d );
    for (unsigned i=1;i<M;i+=1)
    {
      field::poly t1,t2;
      square( d, t1 ); reduce( t1 );
      mul( nzt, t1, t2 ); // N.B. make nzt the first parameter, since it is mostly zero
      square( b, t1 );
      add( t1, t2, b ); reduce( b );
      square( d, t1 ); reduce( t1 );
      add( t1, a, d );
    }
  }
}

unsigned field::rand( unsigned base )
{
  // A very bad prng - but should be good enough for initialisation
  // You can use whatever prng you like here
  // but it needs to be standard so that
  // the same point C.P is always generated
  // for a given root,CB
  prng += 1;
  return prng % base;
}

field::field( full_curve_parameter & a )
: small_field( a.L, a.root ), M(a.L*a.K), K(a.K), T(a.T)
{
  assert( K <= MAXK );
  prng = 0;
  if ( a.tm == 0 )
  {
    // Calculate trace mask
    poly x;
    for (unsigned i=1;i<=K;i+=1)
    {
      x[0] = i;
      unsigned m = 1;
      unsigned w = 0;
      for (unsigned j=0;j<L;j+=1)
      {
        x[i] = m;
        if ( slow_trace(x) )
          w ^= m;
        m *= 2;
      }
      x[i] = 0;
      tm[i] = w;
      if (w) tm[0] = i;
    }
    a.tm = pack( tm );
  }
  else
    unpack( a.tm, tm );

  if (M%2==0)
  {
    nzt[0] = 1;
    nzt[1] = 1;
    while ( (nzt[1] & tm[1]) == 0 )
      nzt[1] *= 2;
    assert( trace( nzt ) != 0 );
  }

#if assert_on // test basic field operations
  for (unsigned i=0;i<10;i+=1)
  {
    field_element a(this),b(this),c(this);
    set_random(a.v); set_random(b.v); set_random(c.v);
    assert( a+b == b+a );
    assert( a*b == b*a );
    assert( a*(b+c) == a*b + a*c );
    if ( b.v[0] )
      assert( a == b*(a/b) );
    assert( ::sqrt(a)*::sqrt(a) == a );
  }
#endif
}

void field::mul( const poly a, const poly b, poly c )
{
  assert( a[0] <= K );
  assert( b[0] <= K );

  if ( a[0]==0 || b[0]==0 )
  {
    c[0] = 0;
    return;
  }
  unsigned i,j;
  c[0] = a[0] + b[0] - 1;
  for (i=1;i<=c[0];i+=1)
    c[i] = 0;

  // pre-computation of log[b[j]] (to reduce table lookups)
  unsigned lb[MAXK+1];
  unsigned b0 = b[0];
  for (j=1;j<=b0;j+=1)
    lb[j] = log[b[j]];

  for (i=1;i<=a[0];i+=1)
  {
    unsigned lai = log[a[i]];
    if ( lai != BASE_M1 )
    {
      for (j=1;j<=b0;j+=1)
      {
        unsigned lbj = lb[j];
        if ( lbj != BASE_M1 )
        {
          unsigned x = lai + lbj;
          if ( x >= BASE_M1 ) x -= BASE_M1;
          c[i+j-1] ^= alog[x];
        }
      }
    }
  }
}

void field::square( const poly a, poly c )
{
  if ( a[0]==0 )
  {
    c[0] = 0;
    return;
  }
  unsigned i;
  c[0] = 2*a[0] - 1;
  for (i=1;;i+=1)
  {
    unsigned x = 0;
    if ( a[i] )
    {
      x = 2*log[a[i]];
      if ( x >= BASE_M1 ) x -= BASE_M1;
      x = alog[x];
    }
    c[2*i-1] = x;
    if (i==a[0]) break;
    c[2*i] = 0;
  }
}

void field::addmul( poly a, unsigned alpha, unsigned j, const poly b )
{
  unsigned la = log[alpha];
  for (unsigned i=1;i<=b[0];i+=1)
  {
    while ( a[0] < j+i ){ a[0] += 1; a[a[0]] = 0; }
    unsigned x = b[i];
    if (x)
    {
      x = log[x] + la;
      if ( x >= BASE_M1 ) x -= BASE_M1;
      a[i+j] ^= alog[x];
    }
  }
  while ( a[0] && a[a[0]]==0 ) a[0] -= 1;
}

void field::div( poly a, unsigned b )
{
  unsigned lb = log[b];
  for (unsigned i=1;i<=a[0];i+=1)
  {
    unsigned ix = a[i];
    if (ix)
    {
      ix = log[ix] + BASE_M1 - lb;
      if ( ix >= BASE - 1 ) ix -= BASE_M1;
      a[i] = alog[ix];
    }
  }
}

void field::inverse( const poly a, poly b )
{
  field::poly c,f,g;
  b[0] = 1; b[1] = 1;
  c[0] = 0;
  g[0] = K+1;
  copy( a, f );
  for (unsigned i=2;i<=K;i+=1)
    g[i] = 0;
  g[1] = 1; g[T+1] = 1; g[K+1] = 1;
  while (1)
  {
    if ( f[0]==1 )
    {
      div(b,f[1]);
      return;
    }
    if ( f[0]<g[0] ) // basically same code with b,c,f,g swapped
    {
      while (1)
      {
        unsigned j = g[0] - f[0];
        unsigned ix = log[g[g[0]]] - log[f[f[0]]] + BASE_M1;
        if ( ix >= BASE_M1 ) ix -= BASE_M1;
        unsigned alpha = alog[ix];
        addmul( g, alpha, j, f );
        addmul( c, alpha, j, b );
        if ( g[0]==1 )
        {
          div( c, g[1] );
          copy( c, b );
          return;
        }
        if ( g[0]<f[0] ) break;
      }
    }
    unsigned j = f[0] - g[0];
    unsigned ix = log[f[f[0]]] - log[g[g[0]]] + BASE_M1;
    if ( ix >= BASE_M1 ) ix -= BASE_M1;
    unsigned alpha = alog[ix];
    addmul( f, alpha, j, g );
    addmul( b, alpha, j, c );