/*
  Author:  Pate Williams (c) 1997

  Exercise IV.3.10
  "Using the Silver-Pohlig-Hellman algorithm,
  find the discrete logarithm to the base 2
  in F_181*. (2 is a generator of F_181*.)"
  -Neal Koblitz-
  See "A Course in Number Theory and Cryptography"
  by Neal Koblitz second edition page 110.
*/

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

#define BITS_PER_LONG   32
#define BITS_PER_LONG_1 31
#define SIZE 32

struct factor {long expon, prime;};

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 get_bit(long i, long *sieve)
{
  long b = i % BITS_PER_LONG;
  long c = i / BITS_PER_LONG;

  return (sieve[c] >> (BITS_PER_LONG_1 - b)) & 1;
}

void set_bit(long i, long v, long *sieve)
{
  long b = i % BITS_PER_LONG;
  long c = i / BITS_PER_LONG;
  long mask = 1 << (BITS_PER_LONG_1 - b);

  if (v == 1)
    sieve[c] |= mask;
  else
    sieve[c] &= ~mask;
}

void Sieve(long n, long *sieve)
{
  long c, i, inc;

  set_bit(0l, 0l, sieve);
  set_bit(1l, 0l, sieve);
  set_bit(2l, 1l, sieve);
  for (i = 3; i <= n; i++)
    set_bit(i, i & 1, sieve);
  c = 3;
  do {
    i = c * c, inc = c + c;
    while (i <= n) {
      set_bit(i, 0l, sieve);
      i += inc;
    }
    c += 2;
    while (!get_bit(c, sieve)) c++;
  } while (c * c <= n);
}

void factor(long n, long *sieve, struct factor *f,
            long *count)
/* factor using trial division */
{
  int one = 0;
  long e, p = 2, s = sqrt(n);

  *count = 0;
  while (!one && p <= s) {
    while (!get_bit(p, sieve)) p++;
    if (n % p == 0) {
      e = 0;
      do {
        e++;
        n = n / p;
      } while (n % p == 0);
      f[*count].expon = e;
      f[*count].prime = p;
      *count = *count + 1;
      one = n == 1;
    }
    p++;
  }
}

long Chinese_Remainder(long r, long *a, long *m)
{
  long i, N = 1, M[SIZE], y[SIZE], x = 0;

  for (i = 0; i < r; i++) N *= m[i];
  for (i = 0; i < r; i++) {
    M[i] = N / m[i];
    y[i] = Extended_Euclidean(M[i], m[i]);
    x += (a[i] * M[i] * y[i]) % N;
  }
  return x;
}

long Pohlig_Hellman(long alpha, long beta,
                    long c, long p, long q)
{
  long ai, delta, e1, s = 0;
  long i, j, p1 = p - 1, p2 = p1 / q, q1 = q;
  long *a = calloc(c, sizeof(long));
  long *betaj = calloc(c, sizeof(long));
  long *gamma = calloc(q, sizeof(long));

  if (!betaj || !gamma) {
    fprintf(stderr, "*error*\ninsufficient memory\n");
    fprintf(stderr, "from Pohlig_Hellman\n");
    exit(1);
  }
  ai = Extended_Euclidean(alpha, p);
  for (i = 0; i < q; i++)
    gamma[i] = exp_mod(alpha, i * p2, p);
  j = 0;
  betaj[j] = beta;
  while (j < c) {
    delta = exp_mod(betaj[j], p1 / q1, p);
    i = 0;
    while (gamma[i] != delta) i++;
    a[j] = i;
    e1 = exp_mod(ai, i * q1 / q, p);
    betaj[j + 1] = (betaj[j] * e1) % p;
    q1 *= q;
    j++;
  }
  q1 = 1;
  for (i = 0; i < c; i++) {
    s += a[i] * q1;
    q1 *= q;
  }
  return s % exp_mod(q, c, p);
}

int main(void)
{
  long alpha = 2;
  long beta = 153;
  long p = 181;
  long a[32], m[32];
  long count, i, log, q, sieve[10000];
  struct factor f[32];

  Sieve(p, sieve);
  q = p - 1;
  factor(q, sieve, f, &count);
  printf("the factorization of %ld is:\n", q);
  for (i = 0; i < count; i++) {
    printf("%ld", f[i].prime);
    if (f[i].expon > 1) printf(" ^ %ld\n", f[i].expon);
    else printf("\n");
  }
  for (i = 0; i < count; i++) {
    a[i] = Pohlig_Hellman(alpha, beta,
                          f[i].expon, p,
                          f[i].prime);
    m[i] = exp_mod(f[i].prime, f[i].expon, p);
  }
  log = Chinese_Remainder(count, a, m) % q;
  printf("log(%ld, %ld) = %ld in F_%ld*\n",
         alpha, beta, log, p);
  if (exp_mod(alpha, log, p) != beta)
    printf("*error*\nin calculation\n");
  return 0;
}