/*
  Author:  Pate Williams (c) 1997

  Exercise 1.6
  "Here is how we might cryptanalyze the Hill Cipher
  using a ciphertext-only attack. Suppose we know
  that m = 2. Break the ciphertext into blocks of
  length two letters (digrams). Each such digram
  is the encryption of a plaintext digram using the
  unknown encryption matrix. Pick out the most
  frequent ciphertext digram and assume it is the
  encryption of a common digram in the list follow-
  ing Table 1.1 (for example TH or ST). For each
  such guess, proceed as in the known-plaintext
  attack, until the correct encryption matrix is
  found. Here is a sample of ciphertext for you to
  decrypt using this method:
  LMQETXYEAGTXCTUIEWNCTXLZEWUAISPZYVAPEWLMGQWYA
  XFTCJMSQCADAGTXLMDXNXSNPJQSYVAPRIQSMHNOCVAXFV"
  -Douglas R. Stinson-
  See "Cryptography: Theory and Practice by
  Douglas R. Stinson page 42.
*/

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

struct digram {long alpha1, alpha2, count;};

long **create_square_matrix(long n)
{
  long i, **matrix = calloc(n, sizeof(long *));

  if (!matrix) {
    fprintf(stderr, "fatal error\ninsufficient memory\n");
    fprintf(stderr, "from create_matrix\n");
    exit(1);
  }
  for (i = 0; i < n; i++) {
    matrix[i] = calloc(n, sizeof(long));
    if (!matrix[i]) {
      fprintf(stderr, "fatal error\ninsufficient memory\n");
      fprintf(stderr, "from create_matrix\n");
      exit(1);
    }
  }
  return matrix;
}

void delete_square_matrix(long n, long **matrix)
{
  long i;

  for (i = 0; i < n; i++) free(matrix[i]);
  free(matrix);
}

void Euclid_extended(long a, long b, long *u,
                     long *v, long *d)
{
  long q, t1, t3, v1, v3;

  *u = 1, *d = a;
  if (b == 0) {
    *v = 0;
    return;
  }
  v1 = 0, v3 = b;
  #ifdef DEBUG
  printf("----------------------------------\n");
  printf(" q    t3   *u   *d   t1   v1   v3\n");
  printf("----------------------------------\n");
  #endif
  while (v3 != 0) {
    q = *d / v3;
    t3 = *d - q * v3;
    t1 = *u - q * v1;
    *u = v1, *d = v3;
    #ifdef DEBUG
    printf("%4ld %4ld %4ld ", q, t3, *u);
    printf("%4ld %4ld %4ld %4ld\n", *d, t1, v1, v3);
    #endif
    v1 = t1, v3 = t3;
  }
  *v = (*d - a * *u) / b;
  #ifdef DEBUG
  printf("----------------------------------\n");
  #endif
}

long inv(long number, long modulus)
{
  long d, u, v;

  Euclid_extended(number, modulus, &u, &v, &d);
  if (d == 1) return u;
  return 0;
}

void gaussian_elimination(long n, long p,
                          long *b, long *x,
                          long **m)
{
  int found;
  long *d = calloc(n, sizeof(long)), ck, dj;
  long i, j, k, l, sum, t;

  if (!d) {
    fprintf(stderr, "fatal error\ninsufficient memory\n");
    fprintf(stderr, "from gaussian_elimination\n");
    exit(1);
  }
  for (j = 0; j < n; j++) {
    found = 0, i = j;
    while (!found && i < n) {
      found = m[i][j] != 0 && inv(m[i][j], p) != 0;
      if (!found) i++;
    }
    if (!found) {
      fprintf(stderr, "fatal error\nnon-invertible matrix\n");
      fprintf(stderr, "from gaussian_elimination\n");
      fprintf(stderr, "j = %ld\n", j);
      for (k = 0; k < n; k++) {
        for (l = 0; l < n; l++)
          printf("%2ld ", m[k][l]);
        printf("\n");
      }
      exit(1);
    }
    if (i > j) {
      /* swap elements */
      for (l = j; l < n; l++)
        t = m[i][l], m[i][l] = m[j][l], m[j][l] = t;
      t = b[i], b[i] = b[j], b[j] = t;
    }
    dj = d[j] = inv(m[j][j], p);
    if (dj == 0) {
      fprintf(stderr, "fatal error\nnon-invertible element\n");
      fprintf(stderr, "from gaussian elimination\n");
      fprintf(stderr, "element %ld mod %ld\n", m[j][j], p);
      exit(1);
    }
    for (k = j + 1; k < n; k++) {
      ck = (dj * m[k][j]) % p;
      for (l = j + 1; l < n; l++) {
        m[k][l] = (m[k][l] - ck * m[j][l]) % p;
        if (m[k][l] < 0) m[k][l] += p;
      }
      b[k] = (b[k] - ck * b[j]) % p;
      if (b[k] < 0) b[k] += p;
    }
  }
  for (i = n - 1; i >= 0; i--) {
    sum = 0;
    for (j = i + 1; j < n; j++)
      sum += (m[i][j] * x[j]) % p;
    if (sum < 0) sum += p;
    x[i] = (d[i] * (b[i] - sum)) % p;
    if (x[i] < 0) x[i] += p;
  }
}

int main(void)
{
  char cipher[5][25] = {"LMQETXYEAGTXCTUIEWNC",
                       /*THEKINGWASINTHECOUNT*/
                        "TXLZEWUAISPZYVAPEWLM",
                       /*INGHOUSECOUNTINGOUTH*/
                        "GQWYAXFTCJMSQCADAGTX",
                       /*ISMONEYTHEQUEENWASIN*/
                        "LMDXNXSNPJQSYVAPRIQS",
                       /*THEPARLOUREATINGBREA*/
                        "MHNOCVAXFV"};
                       /*DANDHONEY*/
  char ciphertext[100], lt_ch, rt_ch;
  long count = 0, d_count = 0, i, j, k, n = 4, p = 26;
  long di_freq[26][26] = {{0}};
  long det, b[4], x[4];
  long **I = create_square_matrix(n);
  long **K = create_square_matrix(n);
  struct digram *d, temp;

  for (i = 0; i < 5; i++)
    for (j = 0; j < strlen(cipher[i]); j++)
      ciphertext[count++] = cipher[i][j];
  for (i = 0; i < count; i++)
    ciphertext[i] -= (char) 'A';
  for (i = 0; i < count; i += 2) {
    lt_ch = ciphertext[i];
    rt_ch = ciphertext[i + 1];
    di_freq[lt_ch][rt_ch]++;
  }
  for (i = 0; i < 26; i++)
    for (j = 0; j < 26; j++)
      if (di_freq[i][j] != 0) d_count++;
  d = calloc(d_count, sizeof(struct digram));
  if (!d) {
    fprintf(stderr, "*error*\ninsufficient memory\n");
    exit(1);
  }
  i = 0;
  for (j = 0; j < 26; j++) {
    for (k = 0; k < 26; k++) {
      if (di_freq[j][k] != 0) {
        d[i].alpha1 = j + 'A';
        d[i].alpha2 = k + 'A';
        d[i++].count = di_freq[j][k];
      }
    }
  }
  /* sort the digrams using the selection sort */
  for (i = 0; i < d_count - 1; i++)
    for (j = i + 1; j < d_count; j++)
      if (d[i].count < d[j].count)
        temp = d[i], d[i] = d[j], d[j] = temp;
  for (i = 0; i < d_count; i++)
    printf("%2ld %c %c %ld\n", i, d[i].alpha1, d[i].alpha2, d[i].count);
  K[0][0] = 'I' - 'A';
  K[0][2] = 'N' - 'A';
  K[1][1] = K[0][0];
  K[1][3] = K[0][2];
  K[2][0] = 'T' - 'A';
  K[2][2] = 'H' - 'A';
  K[3][1] = K[2][0];
  K[3][3] = K[2][2];
  b[0] = d[0].alpha1 - 'A';
  b[1] = d[0].alpha2 - 'A';
  b[2] = d[1].alpha1 - 'A';
  b[3] = d[1].alpha2 - 'A';
  gaussian_elimination(n, p, b, x, K);
  K[0][0] = x[0];
  K[0][1] = x[1];
  K[1][0] = x[2];
  K[1][1] = x[3];
  printf("the key matrix is:\n");
  for (i = 0; i < 2; i++) {
    for (j = 0; j < 2; j++)
      printf("%2ld ", K[i][j]);
    printf("\n");
  }
  det = (K[0][0] * K[1][1] - K[0][1] * K[1][0]) % p;
  if (det < 0) det += p;
  det = inv(det, p);
  if (det < 0) det += p;
  printf("det(K) = %ld\n", det);
  I[0][0] = (det * K[1][1]) % p;
  I[1][1] = (det * K[0][0]) % p;
  I[0][1] = - (det * K[0][1]) % p;
  I[1][0] = - (det * K[1][0]) % p;
  if (I[0][1] < 0) I[0][1] += p;
  if (I[1][0] < 0) I[1][0] += p;
  printf("the inverse key matrix is:\n");
  for (i = 0; i < 2; i++) {
    for (j = 0; j < 2; j++)
      printf("%2ld ", I[i][j]);
    printf("\n");
  }
  for (i = 0; i < count / 2; i++) {
    b[0] = I[0][0] * ciphertext[2 * i]
         + I[1][0] * ciphertext[2 * i + 1];
    b[1] = I[0][1] * ciphertext[2 * i]
         + I[1][1] * ciphertext[2 * i + 1];
    b[0] = b[0] % p + 'A';
    b[1] = b[1] % p + 'A';
    printf("%c%c", b[0], b[1]);
    if ((i + 1) % 10 == 0) printf("\n");
  }
  delete_square_matrix(n, I);
  delete_square_matrix(n, K);
  return 0;
}