/*
  Author:  Pate Williams (c) 1997

  Exercise
  "5.6 The field GF(2 ^ 5) can be constructed as
  Z_2[x]/(x ^ 5 + x ^ 2 + 1). Perform the following
  computations in this field.
  (a) Compute (x ^ 4 + x ^ 2) * (x ^ 3 + x + 1).
  (b) Using the Extended Euclidean algorithm
      compute (x ^ 3 + x ^ 2) ^ - 1.
  (c) Using the square and multiply algorithm,
      compute x ^ 25." -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by
  Douglas R. Stinson page 201.
*/

#include <math.h>
#include <stdio.h>

#define SIZE 32l

void poly_mul(long m, long n, long *a, long *b, long *c, long *p)
{
  long ai, bj, i, j, k, sum;

  *p = m + n;
  for (k = 0; k <= *p; k++) {
    sum = 0;
    for (i = 0; i <= k; i++) {
      j = k - i;
      if (i > m) ai = 0; else ai = a[i];
      if (j > n) bj = 0; else bj = b[j];
      sum += ai * bj;
    }
    c[k] = sum;
  }
}

void poly_div(long m, long n, long *u, long *v,
              long *q, long *r, long *p, long *s)
{
  long j, jk, k, nk, vn = v[n];

  for (j = 0; j <= m; j++) r[j] = u[j];
  if (m < n) {
    *p = 0, *s = m;
    q[0] = 0;
  }
  else {
    *p = m - n, *s = n - 1;
    for (k = *p; k >= 0; k--) {
      nk = n + k;
      q[k] = r[nk] * pow(vn, k);
      for (j = nk - 1; j >= 0; j--) {
        jk = j - k;
        if (jk >= 0)
          r[j] = vn * r[j] - r[nk] * v[j - k];
        else
          r[j] = vn * r[j];
      }
    }
    while (*p > 0 && q[*p] == 0) *p = *p - 1;
    while (*s > 0 && r[*s] == 0) *s = *s - 1;
  }
}

void poly_exp_mod(long degreeA, long degreem, long n,
                  long modulus, long *A, long *m, long *s,
                  long *ds)
{
  int zero;
  long dp, dq, dx = degreeA, i;
  long p[SIZE], q[SIZE], x[SIZE];

  *ds = 0, s[0] = 1;
  for (i = 0; i <= dx; i++) x[i] = A[i];
  while (n > 0) {
    if ((n & 1) == 1) {
      /* s = (s * x) % m; */
      poly_mul(*ds, dx, s, x, p, &dp);
      poly_div(dp, degreem, p, m, q, s, &dq, ds);
      for (i = 0; i <= *ds; i++) s[i] %= modulus;
      zero = s[*ds] == 0, i = *ds;
      while (i > 0 && zero) {
        if (zero) *ds = *ds - 1;
        zero = s[--i] == 0;
      }
    }
    n >>= 1;
    /* x = (x * x) % m; */
    poly_mul(dx, dx, x, x, p, &dp);
    poly_div(dp, degreem, p, m, q, x, &dq, &dx);
    for (i = 0; i <= dx; i++) x[i] %= modulus;
    zero = x[dx] == 0, i = dx;
    while (i > 0 && zero) {
      if (zero) dx--;
      zero = x[--i] == 0;
    }
  }
}

void poly_copy(long db, long *a, long *b, long *da)
/* a = b */
{
  long i;

  *da = db;
  for (i = 0; i <= db; i++) a[i] = b[i];
}

int poly_Extended_Euclidean(long db, long dn,
                            long *b, long *n,
                            long *t, long *dt)
{
  int nonzero;
  long db0, dn0, dq, dr, dt0 = 0, dtemp, du, i;
  long b0[SIZE], n0[SIZE], q[SIZE];
  long r[SIZE], t0[SIZE], temp[SIZE], u[SIZE];

  *dt = 0;
  poly_copy(dn, n0, n, &dn0);
  poly_copy(db, b0, b, &db0);
  t0[0] = 0;
  t[0] = 1;
  poly_div(dn0, db0, n0, b0, q, r, &dq, &dr);
  nonzero = r[0] != 0;
  for (i = 1; !nonzero && i <= dr; i++)
    nonzero = r[i] != 0;
  while (nonzero) {
    poly_mul(dq, *dt, q, t, u, &du);
    if (dt0 < du)
      for (i = dt0 + 1; i <= du; i++) t0[i] = 0;
    for (i = 0; i <= du; i++)
      temp[i] = t0[i] - u[i];
    dtemp = du;
    poly_copy(*dt, t0, t, &dt0);
    poly_copy(dtemp, t, temp, dt);
    poly_copy(db0, n0, b0, &dn0);
    poly_copy(dr, b0, r, &db0);
    poly_div(dn0, db0, n0, b0, q, r, &dq, &dr);
    nonzero = r[0] != 0;
    for (i = 1; !nonzero && i <= dr; i++)
      nonzero = r[i] != 0;
  }
  if (db0 != 0 && b0[0] != 1) return 0;
  return 1;
}

void poly_mod(long da, long p, long *a, long *new_da)
{
  int zero;
  long i;

  for (i = 0; i <= da; i++) {
    a[i] %= p;
    if (a[i] < 0) a[i] += p;
  }
  zero = a[da] == 0;
  for (i = da - 1; zero && i >= 0; i--) {
    da--;
    zero = a[i] == 0;
  }
  *new_da = da;
}

void poly_write(char *label, long da, long *a)
{
  long i;

  printf("%s", label);
  for (i = da; i >= 0; i--)
    printf("%ld ", a[i]);
  printf("\n");
}

int main(void)
{
  long da = 4, db = 3, dc = 3, dd = 5;
  long de, dq, dr, ds, p = 2;
  long a[5] = {0, 0, 1, 0, 1};
  long b[4] = {1, 1, 0, 1};
  long c[4] = {0, 0, 1, 1};
  long d[6] = {1, 0, 1, 0, 0, 1};
  long e[SIZE], q[SIZE], r[SIZE], s[SIZE];

  poly_write("A = ", da, a);
  poly_write("B = ", db, b);
  poly_write("C = ", dc, c);
  poly_write("D = ", dd, d);
  poly_mul(da, db, a, b, e, &de);
  poly_mod(de, p, e, &de);
  poly_div(de, dd, e, d, q, r, &dq, &dr);
  poly_mod(dq, p, q, &dq);
  poly_mod(dr, p, r, &dr);
  poly_write("A * B mod 2 = ", de, e);
  poly_write("A * B / D mod 2 = ", dq, q);
  poly_write("A * B mod D, 2  = ", dr, r);
  if (!poly_Extended_Euclidean(dc, dd, c, d, e, &de))
    printf("*error*\in poly_Extended_Euclidean\n");
  poly_mod(de, p, e, &de);
  poly_write("C ^ - 1 mod D = ", de, e);
  poly_mul(dc, de, c, e, s, &ds);
  poly_mod(ds, p, s, &ds);
  poly_div(ds, dd, s, d, q, r, &dq, &dr);
  poly_mod(dr, p, r, &dr);
  poly_write("C * C ^ - 1 mod D, 2 = ", dr, r);
  dq = 1;
  q[0] = 0;
  q[1] = 1;
  poly_exp_mod(dq, dd, 7, p, q, d, s, &ds);
  poly_mod(ds, p, s, &ds);
  poly_write("x ^  7 mod D, 2 = ", ds, s);
  poly_exp_mod(dq, dd, 9, p, q, d, s, &ds);
  poly_mod(ds, p, s, &ds);
  poly_write("x ^  9 mod D, 2 = ", ds, s);
  poly_exp_mod(dq, dd, 10, p, q, d, s, &ds);
  poly_mod(ds, p, s, &ds);
  poly_write("x ^ 10 mod D, 2 = ", ds, s);
  poly_exp_mod(dq, dd, 25, p, q, d, s, &ds);
  poly_mod(ds, p, s, &ds);
  poly_write("x ^ 25 mod D, 2 = ", ds, s);
  return 0;
}