/*
  Author:  Pate Williams (c) 1997

  "Example 5.2. Suppose p = 809, and we wish
  to find log(3, 525)." -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by
  Douglas R. Stinson pages 165 - 166.
*/

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

struct data {long j, alpha_j;};

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;
  }
  for (i = 0; i < m; i++) {
    printf("(%3ld, %3ld) ", table1[i].j,
           table1[i].alpha_j);
    if ((i + 1) % 5 == 0) printf("\n");
  }
  printf("\n");
  for (i = 0; i < m; i++) {
    printf("(%3ld, %3ld) ", table2[i].j,
           table2[i].alpha_j);
    if ((i + 1) % 5 == 0) printf("\n");
  }
  printf("\n");
  /* 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);
      printf("log(%3ld, %3ld) = %2ld * %2ld + %ld\n",
             alpha, beta, m, key.j, d->j);
      printf("              = %3ld in Z_p%3ld\n",
             log, p);
      if (exp_mod(alpha, log, p) == beta) return log;
    }
  }
  return 0;
}

int main(void)
{
  long alpha = 3, beta = 525l, p = 809l, log;

  log = Shanks(alpha, beta, p);
  if (log == 0)
    printf("log(%3ld, %3ld) in Z_%3ld does not exist\n",
           alpha, beta, p);
  return 0;
}