vlong.cpp

ecppki-cpp curve

{
  return value->bit(i);
}

void vlong::setbit(unsigned i)
{
  docopy();
  value->setbit(i);
}

void vlong::clearbit(unsigned i)
{
  docopy();
  value->clearbit(i);
}

int vlong::cf( const vlong & x ) const
{
  int neg = negative && !value->is_zero();
  if ( neg == (x.negative && !x.value->is_zero()) )
    return value->cf( *x.value );
  else if ( neg ) return -1;
  else return +1;
}

vlong::vlong (unsigned x)
{
  value = new vlong_value;
  negative = 0;
  value->init(x);
}

vlong::vlong ( const vlong& x ) // copy constructor
{
  negative = x.negative;
  value = x.value;
  value->share += 1;
}

vlong& vlong::operator =(const vlong& x)
{
  if ( value->share ) value->share -=1; else delete value;
  value = x.value;
  value->share += 1;
  negative = x.negative;
  return *this;
}

vlong::~vlong()
{
  if ( value->share ) value->share -=1; else delete value;
}

unsigned to_unsigned ( const vlong & x ) // conversion to unsigned
{
  return x.value->to_unsigned();
}

vlong& vlong::operator ^=(const vlong& x)
{
  docopy();
  value->xor( *x.value );
  return *this;
}

vlong& vlong::operator &=(const vlong& x)
{
  docopy();
  value->and( *x.value );
  return *this;
}

vlong& vlong::ror(unsigned n)
{
  docopy();
  int carry = value->shr();
  if (carry)
    value->setbit(n-1);
  return *this;
}

vlong& vlong::rol(unsigned n)
{
  docopy();
  int carry = value->bit(n-1);
  if (carry) value->clearbit(n-1);
  value->shl();
  if (carry) value->setbit(0);
  return *this;
}

vlong& vlong::operator >>=( unsigned n ) // divide by 2**n
{
  docopy();
  value->shr(n);
  return *this;
}

vlong& vlong::operator +=(const vlong& x)
{
  if ( negative == x.negative )
  {
    docopy();
    value->add( *x.value );
  }
  else if ( value->cf( *x.value ) >= 0 )
  {
    docopy();
    value->subtract( *x.value );
  }
  else
  {
    vlong tmp = *this;
    *this = x;
    *this += tmp;
  }
  return *this;
}

vlong& vlong::operator -=(const vlong& x)
{
  if ( negative != x.negative )
  {
    docopy();
    value->add( *x.value );
  }
  else if ( value->cf( *x.value ) >= 0 )
  {
    docopy();
    value->subtract( *x.value );
  }
  else
  {
    vlong tmp = *this;
    *this = x;
    *this -= tmp;
    negative = 1 - negative;
  }
  return *this;
}

vlong operator +( const vlong& x, const vlong& y )
{
  vlong result = x;
  result += y;
  return result;
}

vlong operator -( const vlong& x, const vlong& y )
{
  vlong result = x;
  result -= y;
  return result;
}

vlong operator *( const vlong& x, const vlong& y )
{
  vlong result;
  result.value->mul( *x.value, *y.value );
  result.negative = x.negative ^ y.negative;
  return result;
}

vlong operator /( const vlong& x, const vlong& y )
{
  vlong result;
  vlong_value rem;
  result.value->divide( *x.value, *y.value, rem );
  result.negative = x.negative ^ y.negative;
  return result;
}

vlong operator %( const vlong& x, const vlong& y )
{
  vlong result;
  vlong_value divide;
  divide.divide( *x.value, *y.value, *result.value );
  result.negative = x.negative; // not sure about this?
  return result;
}

vlong operator ^( const vlong& x, const vlong& y )
{
  vlong result = x;
  result ^= y;
  return result;
}

vlong operator &( const vlong& x, const vlong& y )
{
  vlong result = x;
  result &= y;
  return result;
}

vlong operator <<( const vlong& x, unsigned n ) // multiply by 2**n
{
  vlong result = x;
  while (n)
  {
    n -= 1;
    result += result;
  }
  return result;
}


vlong abs( const vlong & x )
{
  vlong result = x;
  result.negative = 0;
  return result;
}

int product ( const vlong &X, const vlong &Y )
{
  return X.value->product( *Y.value );
}

vlong pow2( unsigned n )
{
  vlong result;
  result.setbit(n);
  return result;
}

vlong gcd( const vlong &X, const vlong &Y )
{
  vlong x=X, y=Y;
  while (1)
  {
    if ( y == 0 ) return x;
    x = x % y;
    if ( x == 0 ) return y;
    y = y % x;
  }
}

vlong modinv( const vlong &a, const vlong &m ) // modular inverse
// returns i in range 1..m-1 such that i*a = 1 mod m
// a must be in range 1..m-1
{
  vlong j=1,i=0,b=m,c=a,x,y;
  while ( c != 0 )
  {
    x = b / c;
    y = b - x*c;
    b = c;
    c = y;
    y = j;
    j = i - j*x;
    i = y;
  }
  if ( i < 0 )
    i += m;
  return i;
}

class monty // class for montgomery modular exponentiation
{
  vlong m,n1;
  vlong T,k;   // work registers
  unsigned N;  // bits for R
  void mul( vlong &x, const vlong &y );
public:
  vlong R,R1;
  vlong exp( const vlong &x, const vlong &e );
  vlong monty_exp( const vlong &x, const vlong &e );
  monty( const vlong &M );
};

monty::monty( const vlong &M )
{
  m = M;
  N = 0; R = 1; while ( R < M ) { R += R; N += 1; }
  R1 = modinv( R-m, m );
  n1 = R - modinv( m, R );
}

void monty::mul( vlong &x, const vlong &y )
{
  // T = x*y;
  T.value->fast_mul( *x.value, *y.value, N*2 );

  // k = ( T * n1 ) % R;
  k.value->fast_mul( *T.value, *n1.value, N );

  // x = ( T + k*m ) / R;
  x.value->fast_mul( *k.value, *m.value, N*2 );
  x += T;
  x.value->shr( N );

  if (x>=m) x -= m;
}

vlong monty::monty_exp( const vlong &x, const vlong &e )
{
  vlong result = R-m, t = x; // don't convert input
  t.docopy(); // careful not to modify input
  unsigned bits = e.value->bits();
  unsigned i = 0;
  while (1)
  {
    if ( e.value->bit(i) )
      mul( result, t);
    i += 1;
    if ( i == bits ) break;
    mul( t, t );
  }
  return result; // don't convert output
}

vlong monty::exp( const vlong &x, const vlong &e )
{
  return ( monty_exp( (x*R)%m, e ) * R1 ) % m;
}

vlong modexp( const vlong & x, const vlong & e, const vlong & m )
{
  monty me(m);
  return me.exp( x,e );
}

vlong monty_exp( const vlong & x, const vlong & e, const vlong & m )
{
  monty me(m);
  return me.monty_exp( x,e );
}

vlong monty_exp( const vlong & x, const vlong & d, const vlong & m, const vlong &p, const vlong &q )
{
  monty me(m);
  vlong x1 = ( x * me.R1 ) % m; // Transform input

  vlong u = modinv( p, q );
  vlong dp = d % (p-1);
  vlong dq = d % (q-1);

  // Apply chinese remainder theorem
  vlong a = modexp( x1 % p, dp, p );
  vlong b = modexp( x1 % q, dq, q );
  if ( b < a ) b += q;
  vlong result = a + p * ( ((b-a)*u) % q );

  // Transform result
  return ( result * me.R ) % m;
}

static vlong half(const vlong &a, vlong p)
{
  if (a.bit(0))
    return (a+p)/2;
  return a/2;
}

vlong lucas ( vlong P, vlong Z, vlong k, vlong p ) // P^2 - 4Z != 0
{
  vlong D = P*P - 4*Z, U = 1, V = P, U1,V1;
  unsigned i=k.bits()-1;
  while (i)
  {
    i -= 1;
    U1 = U*V; V1 = V*V + D*U*U; U = U1%p; V = half(V1%p,p);
    if ( k.bit(i) )
    {
      U1 = P*U+V; V1 = P*V+D*U; U = half(U1%p,p); V = half(V1%p,p);
    }
  }
  return V;
}

vlong sqrt( vlong g, vlong p )
{
  vlong result;
  if (p%4==3)
    result = modexp( g, p/4+1, p );
  else if (p%8==5)
  {
    vlong gamma = modexp( 2*g, p/8, p);
    vlong i = 2*g*gamma*gamma - 1;
    result = g*gamma*i;
  }
  else
  {
    vlong Z = g;
    vlong P = 1;
    while (1)
    {
      vlong D = (P*P-4*g)%p; if (D<0) D += p;
      if ( D==0 )
      {
        result = half(P,p);
        break;
      }
      if ( modexp( D, (p-1)/2, p ) != 1 )
      {
        result = half( lucas( P, Z, (p+1)/2, p ), p );
        break;
      }
      P += 1; // Is this ok (efficient?)
    }
  }
  result = result % p;
  if ( result < 0 ) result += p;

  return result;
}