/* 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;

}