/*
  Author:  Pate Williams (c) 1997

  Exercise 11.1
  "Write a computer program to compute the key for
  a Shamir (t, w)-threshold scheme implemented in
  Z_p. That is, given t public x-coordinates,
  x_1, x_2,..., x_t, and t y-coordinates
  y_1, y_2,..., y_t, compute the resulting key.
  Use the Lagrange interpolation method, as it is
  easier to program.
  (a) Test your program if p = 31847, t = 5, and
      w = 10, with the following shares:
      x_1  = 413 y_1 = 25439
      x_2  = 432 y_2 = 14847
      x_3  = 451 y_3 = 24780
      x_4  = 470 y_4 =  5910
      x_5  = 489 y_5 = 12734
      x_6  = 508 y_1 = 12492
      x_7  = 527 y_2 = 12555
      x_8  = 546 y_3 = 28578
      x_9  = 565 y_4 = 20806
      x_10 = 584 y_5 = 21462
      Verify that the same key is computed using
      several different subsets of the five shares.
  (b) Having determined the key, compute the
      share that would be given to a participant
      with x-coordinate 10000." -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by
  Douglas R. Stinson pages 359-360.
*/

#include <stdio.h>
#include <stdlib.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;
}

long Lagrange(long p, long t, long *x, long *y)
/* computes the key using the Lagrange
   interpolation formula, returns the key */
{
  long bj, i, j, k, s = 0;

  for (j = 0; j < t; j++) {
    bj = 1;
    for (k = 0; k < t; k++) {
      if (j != k) {
        i = (x[k] - x[j]) % p;
        if (i < 0) i += p;
        bj = (((bj * x[k]) % p) * Extended_Euclidean(i, p)) % p;
      }
    }
    s = (s + (bj * y[j]) % p) % p;
  }
  if (s < 0) s += p;
  return s;
}

long share(long p, long t, long x, long *a)
{
  long j, s;

  s = a[t - 1];
  for (j = t - 2; j >= 0; j--)
    s = ((s * x) % p + a[j]) % p;
  if (s < 0) s += p;
  return s;
}

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;
  while (v3 != 0) {
    q = *d / v3;
    t3 = *d - q * v3;
    t1 = *u - q * v1;
    *u = v1, *d = v3;
    v1 = t1, v3 = t3;
  }
  *v = (*d - a * *u) / b;
}

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;
      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] = inv(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;
        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;
  }
}

void solve(long p, long t, long *a, long *x, long *y)
{
  long i, j, k, q, z;
  long **A = create_square_matrix(t);

  for (i = 0; i < t; i++) {
    z = x[i];
    for (j = 0; j < t; j++) {
      q = 1;
      for (k = 1; k <= j; k++) q = (q * z) % p;
      A[i][j] = q;
    }
  }
  printf("the matrix is:\n");
  for (i = 0; i < t; i++) {
    for (j = 0; j < t; j++)
      printf("%5ld ", A[i][j]);
    printf("\n");
  }
  gaussian_elimination(t, p, y, a, A);
  delete_square_matrix(t, A);
}

int main(void)
{
  long i, p = 31847, r = 10000, t = 5;
  long a[5], u[5], v[5];
  long x[10] = {413, 432, 451, 470, 489,
                508, 527, 546, 565, 584};
  long y[5] = {25439, 14847, 24780, 5910, 12734};
  long z[5] = {12492, 12555, 28578, 20806, 21462};

  for (i = 0; i < 10; i++) {
    printf("%ld ", x[i]);
    if (i < t) printf("%5ld\n", y[i]);
    else printf("%5ld\n", z[i - t]);
  }
  u[0] = x[0], v[0] = y[0];
  u[1] = x[1], v[1] = y[1];
  u[2] = x[4], v[2] = y[4];
  u[3] = x[5], v[3] = z[0];
  u[4] = x[6], v[4] = z[1];
  printf("key = %ld\n", Lagrange(p, t, x, y));
  printf("key = %ld\n", Lagrange(p, t, x + 5, z));
  printf("key = %ld\n", Lagrange(p, t, u, v));
  for (i = 0; i < t; i++) v[i] = y[i];
  solve(p, t, a, x, v);
  printf("the polynomial coefficients are:\n");
  for (i = 0; i < t; i++)
    printf("a[%ld] = %ld\n", i, a[i]);
  for (i = 0; i < t; i++)
    printf("%ld %5ld %5ld\n", x[i], y[i], share(p, t, x[i], a));
  x[4] = r, y[4] = share(p, t, r, a);
  printf("x = %ld y = %ld\n", x[4], y[4]);
  printf("key = %ld\n", Lagrange(p, t, x, y));
  return 0;
}