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]);
}
