scale = 0

/*
 * bc functions for modular group operations.
 * Anton Stiglic, 2000, 2001
 * send bug reports to: stiglic@cs.mcgill.ca 
 */ 


/*
 * modexp ( m, e, n ) returns m^e mod n
 * Right to left modular exponentiation, as described for example in the 
 * Handbook of Applied Cryptography, algo 14.76,
 * http://www.cacr.math.uwaterloo.ca/hac/
 */

define modexp( m, e, n ) 
{
  auto b, a, i, s;
  a = 1;
  s = m;
  while ( e > 0 ) { 
    if ((e % 2) != 0 ) {
      a = a * s % n;
    }
    e /= 2;
    if ( e > 0 ) {
      s = s*s % n;
    }
  }
  return ( a );
}

/*
 * mpower( a, b, c ) returns a^b mod c
 * This was written by Keith Matthews (U. of Queensland)
 * http://www.maths.uq.edu.au/~krm/MP313/handouts.html
 */

define mpower( a, b, c )
{
  auto x, y, z;
  x = a; y = b; z = 1;
  while ( y > 0 ) {
    while ( y%2 == 0) {
      y /= 2;
      x = (x*x)%c;
    }
    y -= 1;
    z = (z*x)%c;
  }
  return (z);
}

/*
 * modinv( a, n ) returns the inverse of a (mod n) or 0 if it doesn't exist 
 */

define modinv( a, n ) 
{
  auto ret[];

  ret=extended_gcd (a, n);

  if (ret[0] > 1) {
    return (0);
  }

  if (ret[1] < 0) {
    return (ret[1] + n);
  }

  return (ret[1]);
}

/*
 * lcm( a, b ) returns the least common multiple  of a and b.
 */

define lcm( a, b )
{
  return (a*b/gcd(a,b));
}

/* 
 * gcd( a, b) returns the greatest common dividor of a and b.
 */

define gcd( a, b ) 
{
  auto ret[];	

  ret=extended_gcd(a, b);

  return (ret[0]);
}

/* 
 * returns the h such that d = a*h + b*k, for the smallest d > 0.
 * (see _extended_gcd( a, b ).
 */
define gcdh( a, b)
{  
  auto ret[];

  ret=extended_gcd(a, b);

  if ( a < b ) {
    return (ret[1]);
  } else {
    return (ret[2]);
  }
	
  return (0);
} 

/*
 * returns the k such that d = a*h + b*k, for the smallest d > 0.
 * (see extended_gcd( a, b ).
 */
define gcdk( a, b)
{
  auto ret[];

  ret=extended_gcd(a, b);

  if ( a < b ) {
    return (ret[2]);
  } else {
    return (ret[1]);
  }

  return (0);
}

/*
 * Internal functions.
 */

define bits( x ) {
  auto b;

  while ( x > 0 ) {
    x /= 2;
    b += 1;
  }
  return( b );
}

/*
 * extended_gcd( a, b ) computes the smallest d > 0 such that d = a*h + b*k,
 * it sets the resulting d into ret[0], h and k are stored in ret[1] and ret[2] 
 * respectively (ret[] is assumed to have been declared by the calling function;
 * if you declared ret[] for some other reason in the scope of the call of this 
 * function, you'll loose the values that where in the array).
 *
 * The idea behind the algo is to maintain, after each iteration, the relations,
 * d = a*x + b*y and d1 = a*x1 + b*y1,
 * while reducing the value of the d's using the well known fact that 
 * gcd(d, d1) = gcd(d1, d mod d1).
 */

define extended_gcd ( a, b ) 
{
  auto d, d1, d2, q, x, x1, x2, y, y1, y2;

  /* 
   * we start with d = b = a*0 + b*1 and d1 = a = a*1 + b*0 when
   * a < b, the "opposite" when a >= b
   */
  if ( a < b ) {
   d1 = a;
   d = b;
  } else {
   d = a;
   d1 = b;
  }
  x = 1;
  y = 0;
  x1 = 0;
  y1 = 1;

  while ( d1 > 0 ) {
    q = d / d1;

    /* (d2, x2, y2) <- (d, x, y) - q(d1, x1, y1) */
    d2 = d - q*d1;   
    x2 = x - q*x1;
    y2 = y - q*y1;

    /* (d, x, y) <- (d1, x1, y1) */
    d = d1;
    x = x1;
    y = y1;

    /* (d1, x1, y1) <- (d2, x2, y2) */
    d1 = d2;
    x1 = x2;
    y1 = y2;
  }

  ret[0] = d; 
  ret[1] = y; 
  ret[2] = x;

  return (ret[0]);
}