/*
  Author:  Pate Williams (c) 1997

  Exercise
  "5.7 We give an example of ElGamal Cryptosystem
  implemented in GF(3 ^ 3). The polynomial
  x ^ 3 + 2x ^ 2 + 1 is irreducible over Z_3[x]
  and hence Z_3[x]/(x ^ 3 + 2x ^ 2 + 1) is the
  field GF(3 ^ 3). We can associate the 26 letters
  of the alphabet with the 26 nonzero field
  elements, and thus encrypt ordinary text in a
  convenient way. We will use the lexicographic
  ordering of the (nonzero) polynomials to set up
  the correspondence. The correspondence is as
  follows:
  A = 1, B = 2, C = x, D = x + 1, E = x + 2,
  F = 2x, G = 2x + 1, H = 2x + 2, I = x ^ 2,
  J = x ^ 2 + 1, K = x ^ 2 + 2, L = x ^ 2 + x,
  M = x ^ 2 + x + 1, N = x ^ 2 + x + 2,
  O = x ^ 2 + 2x, P = x ^ 2 + 2x + 1,
  Q = x ^ 2 + 2x + 2, R = 2x, S = 2x ^ 2 + 1,
  T = 2x ^ 2 + 2, U = 2x ^ x + x,
  V = 2x ^ 2 + x + 1, W = 2x ^ 2 + x + 2,
  X = 2x ^ 2 + 2x, Y = 2x ^ 2 + 2x + 1,
  Z = 2x ^ 2 + 2x + 2
  Suppose Bob uses alpha = x and a = 11 in an
  ElGamal system; then beta = x + 2. Show how
  Bob will decrypt the following string of
  ciphertext:
  (K, H), (P, X), (N, K), (H, R), (T, F), (V, Y),
  (E, H), (F, A), (T, W), (J, D), (U, J)"
  -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by
  Douglas R. Stinson page 201.
*/

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

#define SIZE 32l

struct ciphertext {long y1, y2;};

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;
      if (k != 0)
        q[k] = r[nk] * pow(vn, k);
      else if (vn != 0)
        q[k] = r[nk];
      else
        q[k] = 0;
      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_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_exp_mod(long degreeA, long degreem, long n,
                  long modulus, long *A, long *m, long *s,
                  long *ds)
{
  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_mod(dp, modulus, p, &dp);
      poly_div(dp, degreem, p, m, q, s, &dq, ds);
      poly_mod(*ds, modulus, s, ds);
    }
    n >>= 1;
    /* x = (x * x) % m; */
    poly_mul(dx, dx, x, x, p, &dp);
    poly_mod(dp, modulus, p, &dp);
    poly_div(dp, degreem, p, m, q, x, &dq, &dx);
    poly_mod(dx, modulus, x, &dx);
  }
}

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 poly_Extended_Euclidean(long db, long dn,
                            long *b, long *n,
                            long *t, long *dt)
{
  int count = 0, 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 && count++ < 8) {
    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);
    if (dn0 - db0 > 0)
      poly_div(dn0, db0, n0, b0, q, r, &dq, &dr);
    else
      return db0 == 0 && b[0] == 1;
    nonzero = r[0] != 0;
    for (i = 1; !nonzero && i <= dr; i++)
      nonzero = r[i] != 0;
    if (db0 == 0) return 1;
  }
  return 0;
}


int main(void)
{
  char message[4] = "A ";
  int found, plaintext;
  long a = 11, dm = 3, dy1, dy2, i, j, k;
  long dq, dr, dx, dy, dz;
  long p = 3, m[SIZE] = {1, 0, 2, 1};
  long q[SIZE], r[SIZE], y1[SIZE], y2[SIZE];
  long x[SIZE], y[SIZE], z[SIZE];
  long field[26][4] = {{1, 0, 0, 0}, {2, 0, 0, 0},
                       {0, 1, 0, 0}, {1, 1, 0, 0},
                       {2, 1, 0, 0}, {0, 2, 0, 0},
                       {1, 2, 0, 0}, {2, 2, 0, 0},
                       {0, 0, 1, 0}, {1, 0, 1, 0},
                       {2, 0, 1, 0}, {0, 1, 1, 0},
                       {1, 1, 1, 0}, {2, 1, 1, 0},
                       {0, 2, 1, 0}, {1, 2, 1, 0},
                       {2, 2, 1, 0}, {0, 0, 2, 0},
                       {1, 0, 2, 0}, {2, 0, 2, 0},
                       {0, 1, 2, 0}, {1, 1, 2, 0},
                       {2, 1, 2, 0}, {0, 2, 2, 0},
                       {1, 2, 2, 0}, {2, 2, 2, 0}};
  struct ciphertext c[11] = {{'K', 'H'}, {'P', 'X'},
                             {'N', 'K'}, {'H', 'R'},
                             {'T', 'F'}, {'V', 'Y'},
                             {'E', 'H'}, {'F', 'A'},
                             {'T', 'W'}, {'J', 'D'},
                             {'U', 'J'}};

  for (i = 0; i < 26; i++) {
    message[0] = (char) (i + 'A');
    poly_write(message, 2, field[i]);
  }
  for (i = 0; i < 11; i++) {
    printf("(%c, %c) ", c[i].y1, c[i].y2);
    if ((i + 1) % 5 == 0) printf("\n");
  }
  printf("\n");
  for (i = 0; i < 11; i++) {
    poly_copy(2, y1, field[c[i].y1 - 'A'], &dy1);
    poly_copy(2, y2, field[c[i].y2 - 'A'], &dy2);
    poly_mod(dy1, p, y1, &dy1);
    poly_mod(dy2, p, y2, &dy2);
    poly_exp_mod(dy1, dm, a, p, y1, m, x, &dx);
    poly_Extended_Euclidean(dx, dm, x, m, y, &dy);
    poly_mod(dy, p, y, &dy);
    poly_mul(dy2, dy, y2, y, z, &dz);
    poly_mod(dz, p, z, &dz);
    poly_div(dz, dm, z, m, q, r, &dq, &dr);
    poly_mod(dr, p, r, &dr);
    found = 0;
    for (j = 0; !found && j < 26; j++) {
      found = r[0] == field[j][0];
      for (k = 1; found && k <= dr; k++)
        found = r[k] == field[j][k];
      if (found) plaintext = 'A' + j;
    }
    printf("%c", plaintext);
  }
  printf("\n");
  return 0;
}