/*
  Author:  Pate Williams (c) 1997

  Exercise 8.2
  "Suppose the Blom Scheme with k = 2 is implemented
  for a set of five users, U, V, W, X and Y. Suppose
  that p = 97, r_u = 14, r_v = 38, r_w = 93, r_x = 69,
  and r_y = 70. The secret g polnomials are as follows:
  g_u(x) = 15 + 15x + 2x^2
  g_v(x) = 95 + 77x + 83x^2
  g_w(x) = 88 + 32x + 18x^2
  g_x(x) = 62 + 91x + 59x^2
  g_y(x) = 10 + 82x + 52x^2
  (a) Show how U and V will compute the key
      K_u,v = K_v,u.
  (b) Show how W, X, and Y together can compute K_u,v."
  -Douglas R. Stinson-
  (b) Show that W and X together can compute
      K_u,v." -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by
  Douglas R. Stinson pages 281-282. The polynomial
  is:
  f(x, y) = a + b(x + y) + c(x^2 + y^2) + dxy
          + e(xy^2 + x^2y) + fx^2y^2
  Any three users can collaborate to determine
  the coefficients of the polynomial by solving
  a nine by six system of overdetermined
  equations.
*/

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

#define SIZE 10

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 Blom_Scheme(long k, long m, long p,
                 long *r, long **g, long **K)
{
  long i, j, l, s, t;

  for (i = 0; i < m; i++) {
    for (j = i + 1; j < m; j++) {
      t = r[j];
      s = g[i][k];
      for (l = k - 1; l >= 0; l--)
        s = t * s + g[i][l];
      K[i][j] = s % p;
    }
  }
  for (i = 0; i < m; i++) {
    for (j = 0; j < i; j++) {
      t = r[j];
      s = g[i][k];
      for (l = k - 1; l >= 0; l--)
        s = t * s + g[i][l];
      s %= p;
      if (s != K[j][i]) {
        fprintf(stderr, "fatal error\nin key calculation\n");
        exit(1);
      }
      K[i][j] = s;
    }
  }
}

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;
      if (!found) i++;
    }
    if (!found) {
      fprintf(stderr, "fatal error\nnon-invertible matrix\n");
      fprintf(stderr, "from gaussian_elimination\n");
      for (k = 0; k < n; k++) {
        for (l = 0; l < n; l++)
          printf("%ld ", 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] = Extended_Euclidean(m[j][j], p);
    if (dj == 0) {
      fprintf(stderr, "fatal error\nnon-invertlible element\n");
      fprintf(stderr, "from gaussian elimination");
      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) % p;
        if (m[k][l] < 0) m[k][l] += p;
      }
      b[k] = (b[k] - (ck * b[j]) % p) % 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)) % p;
    if (x[i] < 0) x[i] += p;
  }
}

void print_matrix(long m, long n, long **matrix)
{
  long i, j;

  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++)
      printf("%2ld ", matrix[i][j]);
    printf("\n");
  }
}

int main(void)
{
  long i, j, k = 2, m = 9, l, mm = 5, n = 6, p = 97, s;
  long b[SIZE], c[SIZE], r[SIZE], x[SIZE];
  long **a, **g, **h, **K;

  a = calloc(m, sizeof(long *));
  g = calloc(m, sizeof(long *));
  h = calloc(m, sizeof(long *));
  K = calloc(m, sizeof(long *));
  if (!a || !g || !h || !K) {
    fprintf(stderr, "fatal error\ninsufficient memory\n");
    exit(1);
  }
  for (i = 0; i < m; i++) {
    a[i] = calloc(n, sizeof(long));
    g[i] = calloc(k + 1, sizeof(long));
    h[i] = calloc(n, sizeof(long));
    K[i] = calloc(m, sizeof(long));
    if (!a[i] || !g[i] || !h[i] || !K[i]) {
      fprintf(stderr, "fatal error\ninsufficient memory\n");
      exit(1);
    }
  }
  g[0][0] = 15, g[0][1] = 15, g[0][2] = 2;
  g[1][0] = 95, g[1][1] = 77, g[1][2] = 83;
  g[2][0] = 88, g[2][1] = 32, g[2][2] = 18;
  g[3][0] = 62, g[3][1] = 91, g[3][2] = 59;
  g[4][0] = 10, g[4][1] = 82, g[4][2] = 52;
  r[0] = 14;
  r[1] = 38;
  r[2] = 92;
  r[3] = 69;
  r[4] = 70;
  Blom_Scheme(k, mm, p, r, g, K);
  printf("the g matrix is:\n");
  print_matrix(mm, k + 1, g);
  printf("the r vector is:\n");
  for (i = 0; i < mm; i++)
    printf("%2ld ", r[i]);
  printf("\n");
  printf("the K matrix is:\n");
  print_matrix(mm, mm, K);
  a[0][0] = 1;
  a[0][1] = r[2];
  a[0][2] = r[2] * r[2];
  a[1][1] = 1;
  a[1][3] = r[2];
  a[1][4] = r[2] * r[2];
  a[2][2] = 1;
  a[2][4] = r[2];
  a[2][5] = r[2] * r[2];
  a[3][0] = 1;
  a[3][1] = r[3];
  a[3][2] = r[3] * r[3];
  a[4][1] = 1;
  a[4][3] = r[3];
  a[4][4] = r[3] * r[3];
  a[5][2] = 1;
  a[5][4] = r[3];
  a[5][5] = r[3] * r[3];
  a[6][0] = 1;
  a[6][1] = r[4];
  a[6][2] = r[4] * r[4];
  a[7][1] = 1;
  a[7][3] = r[4];
  a[7][4] = r[4] * r[4];
  a[8][2] = 1;
  a[8][4] = r[4];
  a[8][5] = r[4] * r[4];
  b[0] = g[2][0];
  b[1] = g[2][1];
  b[2] = g[2][2];
  b[3] = g[3][0];
  b[4] = g[3][1];
  b[5] = g[3][2];
  b[6] = g[4][0];
  b[7] = g[4][1];
  b[8] = g[4][2];
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      s = 0;
      for (l = 0; l < m; l++)
        s += (a[l][i] * a[l][j]) % p;
      s %= p;
      h[i][j] = s;
    }
  }
  printf("the A matrix:\n");
  print_matrix(m, n, a);
  printf("A_t * A:\n");
  print_matrix(n, n, h);
  for (i = 0; i < n; i++) {
    s = 0;
    for (j = 0; j < m; j++)
      s += (a[j][i] * b[j]) % p;
    s %= p;
    c[i] = s;
  }
  gaussian_elimination(n, p, c, x, h);
  printf("a, b, c, d, e, f calculated:\n");
  for (i = 0; i < n; i++)
    printf("%ld ", x[i]);
  printf("\n");
  h[0][0] = x[0] + x[1] * r[0] + x[2] * r[0] * r[0];
  h[0][1] = x[1] + x[3] * r[0] + x[4] * r[0] * r[0];
  h[0][2] = x[2] + x[4] * r[0] + x[5] * r[0] * r[0];
  printf("g_u given and calculated\n");
  for (i = 0; i <= k; i++)
    printf("%2ld %2ld\n", g[0][i], h[0][i] % p);
  for (i = 0; i < m; i++) {
    free(a[i]);
    free(g[i]);
    free(h[i]);
    free(K[i]);
  }
  free(a);
  free(g);
  free(h);
  free(K);
  return 0;
}