/*
  Author:  Pate Williams (c) 1997

  Exercise II.1.4
  "For each degree d <= 6, find the number of
  irreducible polynomials over F_2 of degree d
  and make a list of them." -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

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

void to_binary(long count, long n, long *digit)
{
  long i;

  for (i = 0; i < count; i++) {
    digit[i] = n & 1;
    n >>= 1;
  }
}

void Sieve(char *sieve, long n)
{
  long c, i, inc;

  sieve[2] = 1;
  for (i = 3; i <= n; i++)
    sieve[i] = (char)((i & 1) == 1);
  c = 3;
  do {
    i = c * c, inc = c + c;
    while (i <= n) {
      sieve[i] = 0;
      i += inc;
    }
    c += 2;
    while (!sieve[c]) c++;
  } while (c * c <= n);
}

int irreducible(long n, long p, long *A)
{
  char sieve[SIZE];
  long B[SIZE], C[SIZE], X[SIZE] = {0}, dB, dC;
  long e = pow(p, n), prime = 2;

  Sieve(sieve, n);
  X[1] = 1;
  poly_exp_mod(1, n, e, p, X, A, B, &dB);
  if (dB != 1) return 0;
  if (B[0] != 0) return 0;
  if (B[1] < 0) B[1] += p;
  if (B[1] != 1) return 0;
  while (prime <= n) {
    while (!sieve[prime]) prime++;
    if (n % prime == 0) {
      e = pow(p, n / prime);
      poly_exp_mod(1, n, e, p, X, A, B, &dB);
      if (B[1] < 0) B[1] += p;
      if (dB == 1 && B[1] == 0) dB = 0;
      B[1] -= 1;
      if (dB == 1 && B[1] == 0) dB = 0;
      poly_gcd(n, dB, p, A, B, C, &dC);
      poly_mod(dC, p, C, &dC);
      if (dC != 0) return 0;
    }
    prime++;
  }
  return 1;
}

int main(void)
{
  long A[SIZE], c, d, i, n, p = 2;

  printf("A = 1 0\n");
  printf("A = 1 1\n");
  for (d = 2; d <= 6; d++) {
    n = pow(2, d);
    c = n / 2;
    for (i = 0; i < n; i++) {
      to_binary(c, i, A);
      A[d] = 1;
      if (irreducible(d, p, A))
        poly_write("A = ", d, A);
    }
  }
  return 0;
}