/*
  Author:  Pate Williams (c) 1997

  Exercise II.1.8
  "By computing g.c.d.(f, f'), find all the
  multiple roots of
  f(X) = X ^ 7 + X ^ 5 + X ^ 4 - X ^ 3 - X ^ 2
       - X + 1 contained in F_3[X] in its
  splitting field." -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 < 0) l += p;
  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_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");
}

int main(void)
{
  long df = 7, dg = 6, dh, i, p = 3;
  long f[SIZE] = {0}, g[SIZE] ={0}, h[SIZE];

  f[7] = f[5] = f[4] = 1, f[3] = f[2] = f[1] = 2;
  f[0] = 1;
  g[6] = 7, g[4] = 5, g[3] = 4, g[2] = 6;
  g[1] = 4, g[0] = 2;
  for (i = 0; i <= 7; i++) g[i] %= p;
  poly_write("f  = ", df, f);
  poly_write("f' = ", dg, g);
  poly_gcd(df, dg, p, f, g, h, &dh);
  poly_monic(dh, p, h);
  poly_write("gcd(f, f') = ", dh, h);
  return 0;
}