/*
  Author:  Pate Williams (c) 1997

  Exercise 6.1
  "Suppose Bob is using the ElGamal Signature
  Scheme, and he signs two messages x1 and x2 with
  signatures (gamma_1, delta_1) and (gamma_2, delta_2),
  respectively. (The same value for gamma occurs in
  both signatures.) Suppose also that
  gcd(delta_1 - delta_2, p - 1) = 1.
  (c) Suppose p = 31847, alpha = 5 and beta = 25703.
  Perform the computations of k and a, given the
  signature (23972, 31396) for the message x = 8990
  and the signature (23972, 20481) for the message
  x = 31415." -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by Douglas
  R. Stinson page 231.
*/

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

struct data {long j, alpha_j;};
struct signature {long gamma, delta;};

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 gcd(long a, long b)
{
  long r;

  while (b != 0)
    r = a % b, a = b, b = r;
  return a;
}

int main(void)
{
  long A, B, a0;
  long a, alpha = 5, beta = 25703, p = 31847;
  long d, delta_p, epsilon, i, k, xp;
  long gamma, p1 = p - 1, x1 = 8990, x2 = 31415;
  struct signature s1 = {23972, 31396};
  struct signature s2 = {23972, 20481};

  delta_p = s1.delta - s2.delta;
  if (delta_p < 0) delta_p += p1;
  xp = x1 - x2;
  if (xp < 0) xp += p1;
  epsilon = Extended_Euclidean(delta_p, p1);
  k = (xp * epsilon) % p1;
  if (k < 0) k += p1;
  gamma = exp_mod(alpha, k, p);
  d = gcd(gamma, p1);
  A = gamma / d;
  B = (x1 - (k * s1.delta) % p1) / d;
  p1 /= d;
  if (B < 0) B += p1;
  a0 = (Extended_Euclidean(A, p1) * B) % p1;
  for (i = 0; i < d; i++) {
    a = a0 + i * p1;
    if (beta == exp_mod(alpha, a, p)) break;
  }
  printf("gcd(delta_1 - delta_2, p - 1) = %ld\n", d);
  printf("(x1 - x2) mod (p - 1) = %ld\n", xp);
  printf("epsilon = %ld\n", epsilon);
  printf("gamma = %ld\n", gamma);
  printf("d = %ld\n", d);
  printf("k = %ld\n", k);
  printf("a = %ld\n", a);
  return 0;
}