/*
  Author:  Pate Williams (c) 1997

  Exercise IV.3.1
  "If one has occasion to do a lot of arithmetic
  in a fixed finite field F_q which is not too
  large, it can save time first to compose a
  complete "table of logarithms." In other words,
  choose a generator g of F_q and make a 2-column
  list of all pairs n, g ^ n, as n  goes from
  1 to q - 1; then make third and fourth columns
  listing all pairs a, log_a(a). That is, list
  the elements a of F_q* in some convenient order
  in the third column, and then run down the first
  two columns, putting each n in the fourth column
  next to the a which is g ^ n.
  (a) Make a log table for F_31*, and use it to
  compute 16 * 17, 19 * 13, 1 / 17, 20 / 23.
  (b) Make a log table for F_8*; and use it to
  compute the following: (x + 1) * (x ^ 2 + x),
  (x ^ 2 + x + 1) * (x ^ 2 + 1), 1 / (x ^ 2 + 1),
  x / (x ^ 2 + x + 1)." -Neal Koblitz-
  See "A Course in  Number Theory and Cryptography"
  by Neal Koblitz second edition pages 106-107.
*/

#include <stdio.h>

struct table {long n, gn, a, log_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;
}

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

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

long look_up(long a, struct table *t)
{
  long p = 31, i;

  for (i = 0; i < p - 1; i++)
    if (a == t[i].a) return t[i].log_a;
  return 0;
}

long inverse_look_up(long log_a, struct table *t)
{
  long p = 31, i;

  for (i = 0; i < p - 1; i++)
    if (log_a == t[i].log_a) return t[i].a;
  return 0;
}

int main(void)
{
  int found;
  long a, p = 31, p1 = 30, log_a, i, j;
  struct table t[30];

  for (i = 1; i < p; i++) {
    t[i - 1].n = i;
    t[i - 1].gn = exp_mod(3, i, p);
    t[i - 1].a = i;
  }
  for (i = 0; i < p - 1; i++) {
    a = t[i].a;
    found = 0;
    for (j = 0; !found && j < p - 1; j++) {
      found = a == t[j].gn;
      if (found) t[i].log_a = t[j].n;
    }
  }
  for (i = 0; i < p - 1; i++)
    printf("%2ld %2ld %2ld %2ld\n", t[i].n,
           t[i].gn, t[i].a, t[i].log_a);
  log_a = (look_up(16, t) + look_up(17, t)) % p1;
  printf("16 * 17 = %2ld\n", inverse_look_up(log_a, t));
  printf("16 * 17 = %2ld\n", (16 * 17) % p);
  log_a = (look_up(19, t) + look_up(13, t)) % p1;
  printf("19 * 13 = %2ld\n", inverse_look_up(log_a, t));
  printf("19 * 13 = %2ld\n", (19 * 13) % p);
  log_a = (look_up(1, t) - look_up(17, t)) % p1;
  if (log_a < 0) log_a += p1;
  printf("1  / 17 = %2ld\n", inverse_look_up(log_a, t));
  printf("1  / 17 = %2ld\n", Extended_Euclidean(17, p));
  log_a = (look_up(20, t) - look_up(23, t)) % p1;
  if (log_a < 0) log_a += p1;
  printf("20 / 23 = %2ld\n", inverse_look_up(log_a, t));
  log_a = (20 * Extended_Euclidean(23, p)) % p;
  printf("20 / 23 = %2ld\n", log_a);
  return 0;
}