/*
  Author:  Pate Williams (c) 1997

  Exercise
  "6.2 Suppose I implement the ElGamal Signature
  Scheme with p = 31847, a = 5, and beta = 26379.
  Write a computer program does the following.
  (a) Verify the signature (20679, 11082) on the
      message x = 20543.
  (b) Determine my secret exponent, a, using
      Shanks time-memory tradeoff. Then determine
      the random value k used in signing the
      message x." -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by
  Douglas R. Stinson page 231.
*/

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

#define MAX_COUNT 256l

struct data {long j, alpha_j;};
struct factor {long expon, prime;};

long exp_mod(long x, long b, long n)
/* returns x ^ b mod n */
{
  long a = 1l, s = x;

  while (b != 0) {
    if (b & 1l) a = (a * s) % n;
    b >>= 1;
    if (b != 0) s = (s * s) % n;
  }
  if (a < 0) a += n;
  return a;
}

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;
}

int fcmp(const void *d1, const void *d2)
{
  struct data *data1 = (struct data *) d1;
  struct data *data2 = (struct data *) d2;

  if (data1->alpha_j < data2->alpha_j) return - 1;
  if (data1->alpha_j > data2->alpha_j) return + 1;
  return 0;
}

long Shanks(long alpha, long beta, long p)
/* returns log(alpha, beta) in Z_p where
   alpha is a generator and beta is in Z_p */
{
  long a, alpha_i, i, log, m = ceil(sqrt(p - 1));
  struct data *table1 = calloc(m, sizeof(struct data));
  struct data *table2 = calloc(m, sizeof(struct data));
  struct data *d, key;

  if (!table1 || !table2) {
    fprintf(stderr, "*error*\nin sufficient memory\n");
    fprintf(stderr, "from Shanks\n");
    exit(1);
  }
  /* create a logarithm table */
  alpha_i = exp_mod(alpha, m, p);
  for (i = 0; i < m; i++) {
    table1[i].j = i;
    table1[i].alpha_j = exp_mod(alpha_i, i, p);
  }
  alpha_i = Extended_Euclidean(alpha, p);
  if ((alpha_i * alpha) % p != 1) {
    fprintf(stderr, "*error*\nin Extended_Euclidean\n");
    exit(1);
  }
  for (i = 0; i < m; i++) {
    table2[i].j = i;
    a = exp_mod(alpha_i, i, p);
    table2[i].alpha_j = (beta * a) % p;
  }
  /* sort the tables by second data item */
  qsort(table1, m, sizeof(struct data), fcmp);
  qsort(table2, m, sizeof(struct data), fcmp);
  for (i = 0; i < m; i++) {
    key.j = table1[i].j;
    key.alpha_j = table1[i].alpha_j;
    d = bsearch(&key, table2, m, sizeof(struct data), fcmp);
    if (d) {
      log = (key.j * m + d->j) % (p - 1);
      if (exp_mod(alpha, log, p) == beta) return log;
    }
  }
  return 0;
}

int main(void)
{
  long alpha = 5, beta = 26379, gamma = 20679;
  long a, delta = 11082, k, p = 31847;
  long p1 = p - 1, x = 20543, y, z;

  y = exp_mod(beta, gamma, p);
  z = exp_mod(gamma, delta, p);
  y = (y * z) % p;
  if (y == exp_mod(alpha, x, p))
    printf("signature accepted\n");
  else
    printf("signature rejected\n");
  a = Shanks(alpha, beta, p);
  k = Shanks(alpha, gamma, p);
  printf("a = %ld\n", a);
  printf("k = %ld\n", k);
  y = Extended_Euclidean(k, p1);
  z = ((((x - a * gamma) % p1) * y) % p1) % p1;
  if (z < 0) z += p1;
  if (z != delta)
    printf("error in a, k calculations\n");
  return 0;
}