/*
  Author:  Pate Williams (c) 1997

  Exercise II.1.7
  "Use the polynomial version of the Euclidean
  algorithm to find g.c.d.(f, g) for f, g in
  F_p[X] in each of the following examples. In
  each case express the g.c.d. polynomial as a
  combination of f and g, i.e. in the form
  d(X) = u(X) f(X) + v(X) g(X).
  (a) f = X ^ 3 + X + 1, g = X ^ 2 + X + 1, p = 2;
  (b) f = X ^ 6 + X ^ 5 + X ^ 4 + X ^ 3
        + X ^ 2 + X + 1, g = X ^ 4 + X ^ 2 + X + 1,
      p = 2;
  (c) f = X ^ 3 - X + 1, g = X ^ 2 + 1, p = 3;
  (d) f = X ^ 5 + X ^ 4 + X ^ 3 - X ^ 2 - X + 1,
      g = X ^ 3 + X ^ 2 + X + 1, p = 3;
  (e) f = X ^ 5 + 88 X ^ 4 + 73 X ^ 3 + 83 X ^ 2
        + 51 X + 67, g = X ^ 3 + 97 X ^ 2 + 40 X
        + 38, p = 101." -Neal Koblitz-
  See "A Course in Number Theory and Cryptography"
  by Neal Koblitz second edition page 41.
*/

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

#define SIZE 256

long Extended_Euclidean(long b, long n)
{
  long b0 = b, n0 = n, t = 1, t0 = 0, temp, q, r;

  q = n0 / b0;
  r = n0 - q * b0;
  while (r > 0) {
    temp = t0 - q * t;
    if (temp >= 0) temp = temp % n;
    else temp = n - (- temp % n);
    t0 = t;
    t = temp;
    n0 = b0;
    b0 = r;
    q = n0 / b0;
    r = n0 - q * b0;
  }
  if (b0 != 1) return 0;
  else return t % n;
}

void poly_monic(long da, long p, long *a)
{
  long i, l = a[da];

  if (l != 1) {
    l = Extended_Euclidean(l, p);
    a[da] = 1;
    for (i = da - 1; i >= 0; i--)
      a[i] = (l * a[i]) % p;
  }
}

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

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_Extended_Euclidean(long dg, long dh, long p,
                             long *g, long *h,
                             long *d, long *s,
                             long *t, long *dd,
                             long *ds, long *dt)
{
  int nonzero = 1;
  long da, db, dq, dr, ds1, ds2, dt1, dt2, i;
  long a[SIZE], b[SIZE];
  long q[SIZE], r[SIZE];
  long s1[SIZE], s2[SIZE];
  long t1[SIZE], t2[SIZE];

  if (dh == 0 && h[0] == 0) {
    poly_copy(dg, d, g, dd);
    *ds = *dt = 0;
    s[0] = 1;
    t[0] = 0;
  }
  ds1 = ds2 = dt1 = dt2 = 0;
  s2[0] = 1, s1[0] = 0, t2[0] = 0, t1[0] = 1;
  while (nonzero) {
    poly_div(dg, dh, g, h, q, r, &dq, &dr);
    poly_mod(dq, p, q, &dq);
    poly_mod(dr, p, r, &dr);
    poly_mul(dq, ds1, q, s1, a, &da);
    poly_mod(da, p, a, &da);
    poly_mul(dq, dt1, q, t1, b, &db);
    poly_mod(db, p, b, &db);
    if (ds2 <= da) {
      for (i = 0; i <= ds2; i++)
        s[i] = s2[i] - a[i];
      for (i = ds2 + 1; i <= da; i++) s[i] = - a[i];
    }
    else {
      for (i = 0; i <= da; i++)
        s[i] = s2[i] - a[i];
      for (i = da + 1; i <= ds2; i++) s[i] = s2[i];
    }
    poly_mod(da, p, s, &da);
    if (dt2 <= db) {
      for (i = 0; i <= dt2; i++)
        t[i] = t2[i] - a[i];
      for (i = dt2 + 1; i <= db; i++) t[i] = - b[i];
    }
    else {
      for (i = 0; i <= db; i++)
        t[i] = t2[i] - b[i];
      for (i = db + 1; i <= dt2; i++) t[i] = t2[i];
    }
    poly_mod(db, p, t, &db);
    poly_copy(dh, g, h, &dg);
    poly_copy(dr, h, r, &dh);
    poly_copy(ds1, s2, s1, &ds2);
    poly_copy(da, s1, s, &ds1);
    poly_copy(dt1, t2, t1, &dt2);
    poly_copy(db, t1, t, &dt1);
    nonzero = h[dh] != 0;
    if (!nonzero) {
      for (i = dh - 1; !nonzero && i >= 0; i--) {
        dh--;
        nonzero = h[i] != 0;
      }
    }
  }
  poly_copy(dg, d, g, dd);
  poly_copy(ds2, s, s2, ds);
  poly_copy(dt2, t, t2, dt);
  if (*dd == 0) d[0] = 1;
  for (i = 0; i <= *ds; i++)
    if (s[i] < 0) s[i] += p;
  for (i = 0; i <= *dt; i++)
    if (t[i] < 0) t[i] += p;
  poly_monic(*dd, p, d);
  poly_monic(*ds, p, s);
  poly_monic(*dt, p, t);
}

void poly_gcd(long degreeA, long degreeB, long p,
              long *A, long *B, long *a, long *da)
{
  int nonzero = 0, zero;
  long b[SIZE], db, dq, dr, i, q[SIZE], r[SIZE];

  if (degreeA > degreeB) {
    *da = degreeA;
    db = degreeB;
    for (i = 0; i <= *da; i++) a[i] = A[i];
    for (i = 0; i <= db; i++) b[i] = B[i];
  }
  else {
    *da = degreeB;
    db = degreeA;
    for (i = 0; i <= *da; i++) a[i] = B[i];
    for (i = 0; i <= db; i++) b[i] = A[i];
  }
  for (i = 0; i <= db && !nonzero; i++) nonzero = b[i] != 0;
  while (nonzero) {
    poly_div(*da, db, a, b, q, r, &dq, &dr);
    for (i = 0; i <= db; i++) a[i] = b[i] % p;
    *da = db, zero = a[*da] == 0, i = *da;
    while (i > 0 && zero) {
      if (zero) (*da)--;
      zero = a[--i] == 0;
    }
    for (i = 0; i <= dr; i++) b[i] = r[i] % p;
    db = dr, zero = b[db] == 0, i = db;
    while (i > 0 && zero) {
      if (zero) db--;
      zero = b[--i] == 0;
    }
    nonzero = 0;
    for (i = 0; i <= db && !nonzero; i++) nonzero = b[i] != 0;
  }
}

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

void part(long df, long dg, long p, long *f, long *g)
{
  long dd, du, dv;
  long d[SIZE], u[SIZE], v[SIZE];

  printf("p = %ld\n", p);
  poly_write("f = ", df, f);
  poly_write("g = ", dg, g);
  poly_Extended_Euclidean(df, dg, p, f, g, d, u, v,
                          &dd, &du, &dv);
  poly_write("d = ", dd, d);
  poly_write("u = ", du, u);
  poly_write("v = ", dv, v);
}

int main(void)
{
  long df, dg, f[SIZE], g[SIZE], p;

  df = 3, dg = 2, p = 2;
  f[0] = 1, f[1] = 1, f[2] = 0, f[3] = 1;
  g[0] = 1, g[1] = 1, g[2] = 1;
  part(df, dg, p, f, g);
  df = 6, dg = 4, p = 2;
  f[6] = f[5] = f[4] = f[3] = f[2] = f[1] = f[0] = 1;
  g[4] = 1, g[3] = 0, g[2] = g[1] = g[0] = 1;
  part(df, dg, p, f, g);
  df = 3, dg = 2, p = 3;
  f[3] = 1, f[2] = 0, f[1] = 2, f[0] = 1;
  g[2] = 1, g[1] = 0, g[0] = 1;
  part(df, dg, p, f, g);
  df = 5, dg = 3, p = 3;
  f[5] = f[4] = f[3] = 1, f[2] = f[1] = 2, f[0] = 1;
  g[3] = g[2] = g[1] = g[0] = 1;
  part(df, dg, p, f, g);
  df = 5, dg = 3, p = 101;
  f[5] = 1, f[4] = 88, f[3] = 73, f[2] = 83;
  f[1] = 51, f[0] = 67;
  g[3] = 1, g[2] = 97, g[1] = 40, g[0] = 38;
  part(df, dg, p, f, g);
  return 0;
}