/*
  Author:  Pate Williams (c) 1997

  Exercise
  "5.3 Find log(5, 896) in Z_1103 using the algorithm
  presented in Figure 5.6, given that L_2(beta) = 1
  for beta = 25, 219, 841, and L_2(beta) = 0 for beta
  = 163, 532, 625, and 656." -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by Douglas
  R. Stinson page 200.
*/

#include <stdio.h>

#define BITS_PER_LONG 32l

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 = 1l, s = x;

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

long L_1(long beta, long p2, long p)
{
  return exp_mod(beta, p2, p) != 1;
}

long L_2(long beta)
{
  if (beta == 25 || beta == 219 || beta == 841)
    return 1;
  return 0;
}

long discrete_log(long alpha, long beta, long p)
{
  long ai = Extended_Euclidean(alpha, p);
  long gamma;
  long p2 = (p - 1) >> 1, p4 = (p + 1) >> 2;
  long s, t = 2, x;

  s = x = L_1(beta, p2, p);
  beta = (beta * exp_mod(ai, x, p)) % p;
  printf("beta = %4ld gamma =    0 x = %ld\n",
         beta, x);
  while (beta != 1) {
    x = L_2(beta);
    gamma = exp_mod(beta, p4, p);
    printf("beta = %4ld gamma = %4ld x = %ld\n",
           beta, gamma, x);
    if (L_1(gamma, p2, p) == x) beta = gamma;
    else beta = p - gamma;
    beta = (beta * exp_mod(ai, x, p)) % p;
    s += t * x;
    t *= 2;
  }
  return s;
}

int main(void)
{
  long alpha = 5l, beta = 896l, p = 1103l, log;

  log = discrete_log(alpha, beta, p);
  printf("log(%ld, %ld) = %ld in Z_%ld\n",
         alpha, beta, log, p);
  if (exp_mod(alpha, log, p) != beta)
    printf("*error*\nin computation\n");
  return 0;
}