The Tor Browser / file revision

 /* This Source Code Form is subject to the terms of the Mozilla Public

  * License, v. 2.0. If a copy of the MPL was not distributed with this

  * file, You can obtain one at http://mozilla.org/MPL/2.0/. */

 #include "mp_gf2m.h"

 #include "mp_gf2m-priv.h"

 #include "mplogic.h"

 #include "mpi-priv.h"

 const mp_digit mp_gf2m_sqr_tb[16] =

 {

       0,     1,     4,     5,    16,    17,    20,    21,

      64,    65,    68,    69,    80,    81,    84,    85

 };

 /* Multiply two binary polynomials mp_digits a, b.

  * Result is a polynomial with degree < 2 * MP_DIGIT_BITS - 1.

  * Output in two mp_digits rh, rl.

  */

 #if MP_DIGIT_BITS == 32

 void 

 s_bmul_1x1(mp_digit *rh, mp_digit *rl, const mp_digit a, const mp_digit b)

 {

     register mp_digit h, l, s;

     mp_digit tab[8], top2b = a >> 30; 

     register mp_digit a1, a2, a4;

     a1 = a & (0x3FFFFFFF); a2 = a1 << 1; a4 = a2 << 1;

     tab[0] =  0; tab[1] = a1;    tab[2] = a2;    tab[3] = a1^a2;

     tab[4] = a4; tab[5] = a1^a4; tab[6] = a2^a4; tab[7] = a1^a2^a4;

     s = tab[b       & 0x7]; l  = s;

     s = tab[b >>  3 & 0x7]; l ^= s <<  3; h  = s >> 29;

     s = tab[b >>  6 & 0x7]; l ^= s <<  6; h ^= s >> 26;

     s = tab[b >>  9 & 0x7]; l ^= s <<  9; h ^= s >> 23;

     s = tab[b >> 12 & 0x7]; l ^= s << 12; h ^= s >> 20;

     s = tab[b >> 15 & 0x7]; l ^= s << 15; h ^= s >> 17;

     s = tab[b >> 18 & 0x7]; l ^= s << 18; h ^= s >> 14;

     s = tab[b >> 21 & 0x7]; l ^= s << 21; h ^= s >> 11;

     s = tab[b >> 24 & 0x7]; l ^= s << 24; h ^= s >>  8;

     s = tab[b >> 27 & 0x7]; l ^= s << 27; h ^= s >>  5;

     s = tab[b >> 30      ]; l ^= s << 30; h ^= s >>  2;

     /* compensate for the top two bits of a */

     if (top2b & 01) { l ^= b << 30; h ^= b >> 2; } 

     if (top2b & 02) { l ^= b << 31; h ^= b >> 1; } 

     *rh = h; *rl = l;

 } 

 #else

 void 

 s_bmul_1x1(mp_digit *rh, mp_digit *rl, const mp_digit a, const mp_digit b)

 {

     register mp_digit h, l, s;

     mp_digit tab[16], top3b = a >> 61;

     register mp_digit a1, a2, a4, a8;

     a1 = a & (0x1FFFFFFFFFFFFFFFULL); a2 = a1 << 1; 

     a4 = a2 << 1; a8 = a4 << 1;

     tab[ 0] = 0;     tab[ 1] = a1;       tab[ 2] = a2;       tab[ 3] = a1^a2;

     tab[ 4] = a4;    tab[ 5] = a1^a4;    tab[ 6] = a2^a4;    tab[ 7] = a1^a2^a4;

     tab[ 8] = a8;    tab[ 9] = a1^a8;    tab[10] = a2^a8;    tab[11] = a1^a2^a8;

     tab[12] = a4^a8; tab[13] = a1^a4^a8; tab[14] = a2^a4^a8; tab[15] = a1^a2^a4^a8;

     s = tab[b       & 0xF]; l  = s;

     s = tab[b >>  4 & 0xF]; l ^= s <<  4; h  = s >> 60;

     s = tab[b >>  8 & 0xF]; l ^= s <<  8; h ^= s >> 56;

     s = tab[b >> 12 & 0xF]; l ^= s << 12; h ^= s >> 52;

     s = tab[b >> 16 & 0xF]; l ^= s << 16; h ^= s >> 48;

     s = tab[b >> 20 & 0xF]; l ^= s << 20; h ^= s >> 44;

     s = tab[b >> 24 & 0xF]; l ^= s << 24; h ^= s >> 40;

     s = tab[b >> 28 & 0xF]; l ^= s << 28; h ^= s >> 36;

     s = tab[b >> 32 & 0xF]; l ^= s << 32; h ^= s >> 32;

     s = tab[b >> 36 & 0xF]; l ^= s << 36; h ^= s >> 28;

     s = tab[b >> 40 & 0xF]; l ^= s << 40; h ^= s >> 24;

     s = tab[b >> 44 & 0xF]; l ^= s << 44; h ^= s >> 20;

     s = tab[b >> 48 & 0xF]; l ^= s << 48; h ^= s >> 16;

     s = tab[b >> 52 & 0xF]; l ^= s << 52; h ^= s >> 12;

     s = tab[b >> 56 & 0xF]; l ^= s << 56; h ^= s >>  8;

     s = tab[b >> 60      ]; l ^= s << 60; h ^= s >>  4;

     /* compensate for the top three bits of a */

     if (top3b & 01) { l ^= b << 61; h ^= b >> 3; } 

     if (top3b & 02) { l ^= b << 62; h ^= b >> 2; } 

     if (top3b & 04) { l ^= b << 63; h ^= b >> 1; } 

     *rh = h; *rl = l;

 } 

 #endif

 /* Compute xor-multiply of two binary polynomials  (a1, a0) x (b1, b0)  

  * result is a binary polynomial in 4 mp_digits r[4].

  * The caller MUST ensure that r has the right amount of space allocated.

  */

 void 

 s_bmul_2x2(mp_digit *r, const mp_digit a1, const mp_digit a0, const mp_digit b1,

            const mp_digit b0)

 {

     mp_digit m1, m0;

     /* r[3] = h1, r[2] = h0; r[1] = l1; r[0] = l0 */

     s_bmul_1x1(r+3, r+2, a1, b1);

     s_bmul_1x1(r+1, r, a0, b0);

     s_bmul_1x1(&m1, &m0, a0 ^ a1, b0 ^ b1);

     /* Correction on m1 ^= l1 ^ h1; m0 ^= l0 ^ h0; */

     r[2] ^= m1 ^ r[1] ^ r[3];  /* h0 ^= m1 ^ l1 ^ h1; */

     r[1]  = r[3] ^ r[2] ^ r[0] ^ m1 ^ m0;  /* l1 ^= l0 ^ h0 ^ m0; */

 }

 /* Compute xor-multiply of two binary polynomials  (a2, a1, a0) x (b2, b1, b0)  

  * result is a binary polynomial in 6 mp_digits r[6].

  * The caller MUST ensure that r has the right amount of space allocated.

  */

 void 

 s_bmul_3x3(mp_digit *r, const mp_digit a2, const mp_digit a1, const mp_digit a0, 

 	const mp_digit b2, const mp_digit b1, const mp_digit b0)

 {

 	mp_digit zm[4];

 	s_bmul_1x1(r+5, r+4, a2, b2);         /* fill top 2 words */

 	s_bmul_2x2(zm, a1, a2^a0, b1, b2^b0); /* fill middle 4 words */

 	s_bmul_2x2(r, a1, a0, b1, b0);        /* fill bottom 4 words */

 	zm[3] ^= r[3];

 	zm[2] ^= r[2]; 

 	zm[1] ^= r[1] ^ r[5];

 	zm[0] ^= r[0] ^ r[4];

 	r[5]  ^= zm[3];

 	r[4]  ^= zm[2];

 	r[3]  ^= zm[1];

 	r[2]  ^= zm[0];

 }

 /* Compute xor-multiply of two binary polynomials  (a3, a2, a1, a0) x (b3, b2, b1, b0)  

  * result is a binary polynomial in 8 mp_digits r[8].

  * The caller MUST ensure that r has the right amount of space allocated.

  */

 void s_bmul_4x4(mp_digit *r, const mp_digit a3, const mp_digit a2, const mp_digit a1, 

 	const mp_digit a0, const mp_digit b3, const mp_digit b2, const mp_digit b1, 

 	const mp_digit b0)

 {

 	mp_digit zm[4];

 	s_bmul_2x2(r+4, a3, a2, b3, b2);            /* fill top 4 words */

 	s_bmul_2x2(zm, a3^a1, a2^a0, b3^b1, b2^b0); /* fill middle 4 words */

 	s_bmul_2x2(r, a1, a0, b1, b0);              /* fill bottom 4 words */

 	zm[3] ^= r[3] ^ r[7]; 

 	zm[2] ^= r[2] ^ r[6]; 

 	zm[1] ^= r[1] ^ r[5]; 

 	zm[0] ^= r[0] ^ r[4]; 

 	r[5]  ^= zm[3];    

 	r[4]  ^= zm[2];

 	r[3]  ^= zm[1];    

 	r[2]  ^= zm[0];

 }

 /* Compute addition of two binary polynomials a and b,

  * store result in c; c could be a or b, a and b could be equal; 

  * c is the bitwise XOR of a and b.

  */

 mp_err

 mp_badd(const mp_int *a, const mp_int *b, mp_int *c)

 {

     mp_digit *pa, *pb, *pc;

     mp_size ix;

     mp_size used_pa, used_pb;

     mp_err res = MP_OKAY;

     /* Add all digits up to the precision of b.  If b had more

      * precision than a initially, swap a, b first

      */

     if (MP_USED(a) >= MP_USED(b)) {

         pa = MP_DIGITS(a);

         pb = MP_DIGITS(b);

         used_pa = MP_USED(a);

         used_pb = MP_USED(b);

     } else {

         pa = MP_DIGITS(b);

         pb = MP_DIGITS(a);

         used_pa = MP_USED(b);

         used_pb = MP_USED(a);

     }

     /* Make sure c has enough precision for the output value */

     MP_CHECKOK( s_mp_pad(c, used_pa) );

     /* Do word-by-word xor */

     pc = MP_DIGITS(c);

     for (ix = 0; ix < used_pb; ix++) {

         (*pc++) = (*pa++) ^ (*pb++);

     }

     /* Finish the rest of digits until we're actually done */

     for (; ix < used_pa; ++ix) {

         *pc++ = *pa++;

     }

     MP_USED(c) = used_pa;

     MP_SIGN(c) = ZPOS;

     s_mp_clamp(c);

 CLEANUP:

     return res;

 } 

 #define s_mp_div2(a) MP_CHECKOK( mpl_rsh((a), (a), 1) );

 /* Compute binary polynomial multiply d = a * b */

 static void 

 s_bmul_d(const mp_digit *a, mp_size a_len, mp_digit b, mp_digit *d)

 {

     mp_digit a_i, a0b0, a1b1, carry = 0;

     while (a_len--) {

         a_i = *a++;

         s_bmul_1x1(&a1b1, &a0b0, a_i, b);

         *d++ = a0b0 ^ carry;

         carry = a1b1;

     }

     *d = carry;

 }

 /* Compute binary polynomial xor multiply accumulate d ^= a * b */

 static void 

 s_bmul_d_add(const mp_digit *a, mp_size a_len, mp_digit b, mp_digit *d)

 {

     mp_digit a_i, a0b0, a1b1, carry = 0;

     while (a_len--) {

         a_i = *a++;

         s_bmul_1x1(&a1b1, &a0b0, a_i, b);

         *d++ ^= a0b0 ^ carry;

         carry = a1b1;

     }

     *d ^= carry;

 }

 /* Compute binary polynomial xor multiply c = a * b.  

  * All parameters may be identical.

  */

 mp_err 

 mp_bmul(const mp_int *a, const mp_int *b, mp_int *c)

 {

     mp_digit *pb, b_i;

     mp_int tmp;

     mp_size ib, a_used, b_used;

     mp_err res = MP_OKAY;

     MP_DIGITS(&tmp) = 0;

     ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

     if (a == c) {

         MP_CHECKOK( mp_init_copy(&tmp, a) );

         if (a == b)

             b = &tmp;

         a = &tmp;

     } else if (b == c) {

         MP_CHECKOK( mp_init_copy(&tmp, b) );

         b = &tmp;

     }

     if (MP_USED(a) < MP_USED(b)) {

         const mp_int *xch = b;      /* switch a and b if b longer */

         b = a;

         a = xch;

     }

     MP_USED(c) = 1; MP_DIGIT(c, 0) = 0;

     MP_CHECKOK( s_mp_pad(c, USED(a) + USED(b)) );

     pb = MP_DIGITS(b);

     s_bmul_d(MP_DIGITS(a), MP_USED(a), *pb++, MP_DIGITS(c));

     /* Outer loop:  Digits of b */

     a_used = MP_USED(a);

     b_used = MP_USED(b);

 	MP_USED(c) = a_used + b_used;

     for (ib = 1; ib < b_used; ib++) {

         b_i = *pb++;

         /* Inner product:  Digits of a */

         if (b_i)

             s_bmul_d_add(MP_DIGITS(a), a_used, b_i, MP_DIGITS(c) + ib);

         else

             MP_DIGIT(c, ib + a_used) = b_i;

     }

     s_mp_clamp(c);

     SIGN(c) = ZPOS;

 CLEANUP:

     mp_clear(&tmp);

     return res;

 }

 /* Compute modular reduction of a and store result in r.  

  * r could be a. 

  * For modular arithmetic, the irreducible polynomial f(t) is represented 

  * as an array of int[], where f(t) is of the form: 

  *     f(t) = t^p[0] + t^p[1] + ... + t^p[k]

  * where m = p[0] > p[1] > ... > p[k] = 0.

  */

 mp_err

 mp_bmod(const mp_int *a, const unsigned int p[], mp_int *r)

 {

     int j, k;

     int n, dN, d0, d1;

     mp_digit zz, *z, tmp;

     mp_size used;

     mp_err res = MP_OKAY;

     /* The algorithm does the reduction in place in r, 

      * if a != r, copy a into r first so reduction can be done in r

      */

     if (a != r) {

         MP_CHECKOK( mp_copy(a, r) );

     }

     z = MP_DIGITS(r);

     /* start reduction */

     /*dN = p[0] / MP_DIGIT_BITS; */

     dN = p[0] >> MP_DIGIT_BITS_LOG_2;

     used = MP_USED(r);

     for (j = used - 1; j > dN;) {

         zz = z[j];

         if (zz == 0) {

             j--; continue;

         }

         z[j] = 0;

         for (k = 1; p[k] > 0; k++) {

             /* reducing component t^p[k] */

             n = p[0] - p[k];

             /*d0 = n % MP_DIGIT_BITS;   */

             d0 = n & MP_DIGIT_BITS_MASK;

             d1 = MP_DIGIT_BITS - d0;

             /*n /= MP_DIGIT_BITS; */

             n >>= MP_DIGIT_BITS_LOG_2;

             z[j-n] ^= (zz>>d0);

             if (d0) 

                 z[j-n-1] ^= (zz<<d1);

         }

         /* reducing component t^0 */

         n = dN;  

         /*d0 = p[0] % MP_DIGIT_BITS;*/

         d0 = p[0] & MP_DIGIT_BITS_MASK;

         d1 = MP_DIGIT_BITS - d0;

         z[j-n] ^= (zz >> d0);

         if (d0) 

             z[j-n-1] ^= (zz << d1);

     }

     /* final round of reduction */

     while (j == dN) {

         /* d0 = p[0] % MP_DIGIT_BITS; */

         d0 = p[0] & MP_DIGIT_BITS_MASK;

         zz = z[dN] >> d0;  

         if (zz == 0) break;

         d1 = MP_DIGIT_BITS - d0;

         /* clear up the top d1 bits */

         if (d0) {

 	    z[dN] = (z[dN] << d1) >> d1; 

 	} else {

 	    z[dN] = 0;

 	}

         *z ^= zz; /* reduction t^0 component */

         for (k = 1; p[k] > 0; k++) {

             /* reducing component t^p[k]*/

             /* n = p[k] / MP_DIGIT_BITS; */

             n = p[k] >> MP_DIGIT_BITS_LOG_2;

             /* d0 = p[k] % MP_DIGIT_BITS; */

             d0 = p[k] & MP_DIGIT_BITS_MASK;

             d1 = MP_DIGIT_BITS - d0;

             z[n] ^= (zz << d0);

             tmp = zz >> d1;

             if (d0 && tmp)

                 z[n+1] ^= tmp;

         }

     }

     s_mp_clamp(r);

 CLEANUP:

     return res;

 }

 /* Compute the product of two polynomials a and b, reduce modulo p, 

  * Store the result in r.  r could be a or b; a could be b.

  */

 mp_err 

 mp_bmulmod(const mp_int *a, const mp_int *b, const unsigned int p[], mp_int *r)

 {

     mp_err res;

     if (a == b) return mp_bsqrmod(a, p, r);

     if ((res = mp_bmul(a, b, r) ) != MP_OKAY)

 	return res;

     return mp_bmod(r, p, r);

 }

 /* Compute binary polynomial squaring c = a*a mod p .  

  * Parameter r and a can be identical.

  */

 mp_err 

 mp_bsqrmod(const mp_int *a, const unsigned int p[], mp_int *r)

 {

     mp_digit *pa, *pr, a_i;

     mp_int tmp;

     mp_size ia, a_used;

     mp_err res;

     ARGCHK(a != NULL && r != NULL, MP_BADARG);

     MP_DIGITS(&tmp) = 0;

     if (a == r) {

         MP_CHECKOK( mp_init_copy(&tmp, a) );

         a = &tmp;

     }

     MP_USED(r) = 1; MP_DIGIT(r, 0) = 0;

     MP_CHECKOK( s_mp_pad(r, 2*USED(a)) );

     pa = MP_DIGITS(a);

     pr = MP_DIGITS(r);

     a_used = MP_USED(a);

 	MP_USED(r) = 2 * a_used;

     for (ia = 0; ia < a_used; ia++) {

         a_i = *pa++;

         *pr++ = gf2m_SQR0(a_i);

         *pr++ = gf2m_SQR1(a_i);

     }

     MP_CHECKOK( mp_bmod(r, p, r) );

     s_mp_clamp(r);

     SIGN(r) = ZPOS;

 CLEANUP:

     mp_clear(&tmp);

     return res;

 }

 /* Compute binary polynomial y/x mod p, y divided by x, reduce modulo p.

  * Store the result in r. r could be x or y, and x could equal y.

  * Uses algorithm Modular_Division_GF(2^m) from 

  *     Chang-Shantz, S.  "From Euclid's GCD to Montgomery Multiplication to 

  *     the Great Divide".

  */

 int 

 mp_bdivmod(const mp_int *y, const mp_int *x, const mp_int *pp, 

     const unsigned int p[], mp_int *r)

 {

     mp_int aa, bb, uu;

     mp_int *a, *b, *u, *v;

     mp_err res = MP_OKAY;

     MP_DIGITS(&aa) = 0;

     MP_DIGITS(&bb) = 0;

     MP_DIGITS(&uu) = 0;

     MP_CHECKOK( mp_init_copy(&aa, x) );

     MP_CHECKOK( mp_init_copy(&uu, y) );

     MP_CHECKOK( mp_init_copy(&bb, pp) );

     MP_CHECKOK( s_mp_pad(r, USED(pp)) );

     MP_USED(r) = 1; MP_DIGIT(r, 0) = 0;

     a = &aa; b= &bb; u=&uu; v=r;

     /* reduce x and y mod p */

     MP_CHECKOK( mp_bmod(a, p, a) );

     MP_CHECKOK( mp_bmod(u, p, u) );

     while (!mp_isodd(a)) {

         s_mp_div2(a);

         if (mp_isodd(u)) {

             MP_CHECKOK( mp_badd(u, pp, u) );

         }

         s_mp_div2(u);

     }

     do {

         if (mp_cmp_mag(b, a) > 0) {

             MP_CHECKOK( mp_badd(b, a, b) );

             MP_CHECKOK( mp_badd(v, u, v) );

             do {

                 s_mp_div2(b);

                 if (mp_isodd(v)) {

                     MP_CHECKOK( mp_badd(v, pp, v) );

                 }

                 s_mp_div2(v);

             } while (!mp_isodd(b));

         }

         else if ((MP_DIGIT(a,0) == 1) && (MP_USED(a) == 1))

             break;

         else {

             MP_CHECKOK( mp_badd(a, b, a) );

             MP_CHECKOK( mp_badd(u, v, u) );

             do {

                 s_mp_div2(a);

                 if (mp_isodd(u)) {

                     MP_CHECKOK( mp_badd(u, pp, u) );

                 }

                 s_mp_div2(u);

             } while (!mp_isodd(a));

         }

     } while (1);

     MP_CHECKOK( mp_copy(u, r) );

 CLEANUP:

     mp_clear(&aa);

     mp_clear(&bb);

     mp_clear(&uu);

     return res;

 }

 /* Convert the bit-string representation of a polynomial a into an array

  * of integers corresponding to the bits with non-zero coefficient.

  * Up to max elements of the array will be filled.  Return value is total

  * number of coefficients that would be extracted if array was large enough.

  */

 int

 mp_bpoly2arr(const mp_int *a, unsigned int p[], int max)

 {

     int i, j, k;

     mp_digit top_bit, mask;

     top_bit = 1;

     top_bit <<= MP_DIGIT_BIT - 1;

     for (k = 0; k < max; k++) p[k] = 0;

     k = 0;

     for (i = MP_USED(a) - 1; i >= 0; i--) {

         mask = top_bit;

         for (j = MP_DIGIT_BIT - 1; j >= 0; j--) {

             if (MP_DIGITS(a)[i] & mask) {

                 if (k < max) p[k] = MP_DIGIT_BIT * i + j;

                 k++;

             }

             mask >>= 1;

         }

     }

     return k;

 }

 /* Convert the coefficient array representation of a polynomial to a 

  * bit-string.  The array must be terminated by 0.

  */

 mp_err

 mp_barr2poly(const unsigned int p[], mp_int *a)

 {

     mp_err res = MP_OKAY;

     int i;

     mp_zero(a);

     for (i = 0; p[i] > 0; i++) {

 	MP_CHECKOK( mpl_set_bit(a, p[i], 1) );

     }

     MP_CHECKOK( mpl_set_bit(a, 0, 1) );

 CLEANUP:

     return res;

 }