/******************************************************************************** * * * The functions in this file use the basic large integer arithmetic sub- * * routines to perform such things as greatest common divisor, modular * * inversion and modular exponentiation. Note that many other math packages * * do all these things much better, but this is very easy to follow. * * * ********************************************************************************/ #include #include "bigint.h" /* divide a large integer by 2. A simple shift right operation. */ void int_div2( x) BIGINT *x; { INDEX j; ELEMENT mask; INTLOOP(j) { if (j) mask = ( x->hw[j-1] & 1) ? CARRY : 0; else mask = 0; x->hw[j] = (x->hw[j] | mask) >> 1; } } /* compute greatest common divisor using binary method. See [DEK, pg 321] Algorithm B for theoretical details. Enter with large integers u, v. Returns gcd(u, v) as large integer in w. */ void int_gcd(u, v, w) BIGINT *u, *v, *w; { INDEX k, i, flag; ELEMENT check, carry_bit; BIGINT t, U, V; int_copy( u, &U); int_copy( v, &V); /* find common powers of 2 and eliminate them */ k = 0; /* check that both u and v are even */ while ( !(U.hw[INTMAX] & 1 || V.hw[INTMAX] & 1)) { /* increment power of 2 and divide both u and v by 2 */ k++; int_div2( &U); int_div2( &V); } /* Now both u and v have been divided by 2^k. If u is odd, set t = v and flag, otherwise t = u and clear flag */ if (U.hw[INTMAX] & 1) { int_copy( &V, &t); flag = -1; } else { int_copy( &U, &t); flag = 1; } check = 0; INTLOOP (i) check |= t.hw[i]; while (check) { /* while t is even, divide by 2 */ while ( !(t.hw[INTMAX] & 1)) int_div2( &t); /* reset u or v to t depending on sign of flag */ if (flag > 0) int_copy( &t, &U); else int_copy( &t, &V); /* t = u - v; core reduction step, gcd remains unchanged */ int_sub( &U, &V, &t); if (t.hw[0] & MSB_HW) { flag = -1; int_neg( &t); } else flag = 1; check = 0; INTLOOP (i) check |= t.hw[i]; } /* reapply common powers of 2. First do words, then do bits.*/ int_copy( &U, w); while ( k > HALFSIZE ) { for (i=0; ihw[i] = w->hw[i+1]; k -= HALFSIZE; w->hw[INTMAX] = 0; } carry_bit = 0; while ( k > 0 ) { INTLOOP (i) { w->hw[i] = (w->hw[i] << 1) | carry_bit; carry_bit = w->hw[i] & CARRY ? 1 : 0; w->hw[i] &= LOMASK; } k--; } } /* Binary method for modular exponentiation. Taken from [DEK,pg 442] Algorithm A. Computes z = x^n mod q for x, n and q large integers. */ void mod_exp(x, n, q, z) BIGINT *x, *n, *q, *z; { BIGINT N, Y, Z, temp, dummy; ELEMENT check; INDEX i; /* initialize variables */ int_copy (n, &N); int_null( &Y); Y.hw[INTMAX] = 1; int_copy (x, &Z); /* Main loop divides N by 2 each step. Repeat until N = 0, and return Y as result. */ check = 0; INTLOOP (i) check |= N.hw[i]; while (check) { /* if N is odd, multiply by extra factor of Y */ if (N.hw[INTMAX] & 1) { /* Y = (Y * Z) % q; */ int_mul (&Y, &Z, &temp); int_div (&temp, q, &dummy, &Y); } int_div2( &N); /* divide N by 2 */ /* Z = (Z * Z) % q; square Z */ int_mul (&Z, &Z, &temp); int_div( &temp, q, &dummy, &Z); check = 0; INTLOOP (i) check |= N.hw[i]; } int_copy (&Y, z); } /* Inversion of numbers in a prime field is similar to solving the linear congruence ax = c mod b. Replace c with 1 and x is the inverse of a mod b. Taken from [HR, pg 309]. Inputs are large integer a and modulus b. Output is x, inverse of a mod b (ax = 1 mod b). */ void mod_inv(a, b, x) BIGINT *a, *b, *x; { BIGINT m, n, p0, p1, p2, q, r, temp, dummy; ELEMENT check; INDEX sw, i; /* initialize loop variables sw = 1; m = b; n = a; p0 = 1; p1 = m/n; q = p1; r = m % n; */ sw = 1; int_copy( b, &m); int_copy( a, &n); int_null ( &p0); p0.hw[INTMAX] = 1; int_div ( &m, &n, &p1, &r); int_copy ( &p1, &q); /* main loop, compute continued fraction intermediates */ check = 0; INTLOOP (i) check |= r.hw[i]; while (check) { sw = -sw; int_copy( &n, &m); int_copy( &r, &n); int_div( &m, &n, &q, &r); /* p2 = (q * p1 + p0) % b; core operation of routine */ int_mul( &q, &p1, &temp); int_add( &temp, &p0, &temp); int_div( &temp, b, &dummy, &p2); int_copy( &p1, &p0); int_copy( &p2, &p1); check = 0; INTLOOP (i) check |= r.hw[i]; } /* sw keeps track of sign. If sw < 0, add modulus to result */ if (sw < 0) int_sub( b, &p0, x); else int_copy( &p0, x); }