/*
  Author:  Pate Williams (c) 1997

  Exercise 8.6
  "Consider the Girault Scheme where p = 167, q = 179,
  and hence n = 29893. Suppose alpha = 2 and e = 11101.
  (a) Compute d.
  (b) Given ID(U) = 10021 and a_u = 9843, compute b_u
      and p_u. Given that ID(V) = 10022 and a_v = 7692,
      compute b_v and p_v.
  (c) Show how b_u can be computed from p_u and ID(U)
      using the public exponent e. Similarly, show how
      b_v can be computed from p_v and ID(V).
  (d) Suppose that U chooses r_u = 15556 and V chooses
      r_v = 6420. Compute s_u and s_v, and show how U
      and V each compute a common key."
  -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by Douglas
  R. Stinson page 282.
*/

#include <stdio.h>

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;
}

int main(void)
{
  long alpha = 2, e = 11101, p = 167, q = 179;
  long a_u = 9843, a_v = 7692, c_u, c_v;
  long n = p * q, phi = (p - 1) * (q - 1);
  long d = Extended_Euclidean(e, phi);
  long ID_U = 10021, ID_V = 10022;
  long r_u = 15556, r_v = 6420;
  long b_u, b_v, p_u, p_v, s_u, s_v, K_uv, K_vu;

  b_u = exp_mod(alpha, a_u, n);
  b_v = exp_mod(alpha, a_v, n);
  c_u = b_u - ID_U;
  if (c_u < 0) c_u += n;
  c_v = b_v - ID_V;
  if (c_v < 0) c_v += n;
  p_u = exp_mod(c_u, d, n);
  p_v = exp_mod(c_v, d, n);
  s_u = exp_mod(alpha, r_u, n);
  s_v = exp_mod(alpha, r_v, n);
  c_u = (exp_mod(p_v, e, n) + ID_V) % n;
  c_v = (exp_mod(p_u, e, n) + ID_U) % n;
  K_uv = (exp_mod(s_v, a_u, n) * exp_mod(c_u, r_u, n)) % n;
  K_vu = (exp_mod(s_u, a_v, n) * exp_mod(c_v, r_v, n)) % n;
  c_u = (exp_mod(p_u, e, n) + ID_U) % n;
  c_v = (exp_mod(p_v, e, n) + ID_V) % n;
  printf("alpha = %ld\n", alpha);
  printf("p = %ld\n", p);
  printf("q = %ld\n", q);
  printf("n = %ld\n", n);
  printf("a_u = %ld\n", a_u);
  printf("a_v = %ld\n", a_v);
  printf("ID(U) = %ld\n", ID_U);
  printf("ID(V) = %ld\n", ID_V);
  printf("r_u = %ld\n", r_u);
  printf("r_v = %ld\n", r_v);
  printf("e = %ld\n", e);
  printf("d = %ld\n", d);
  printf("b_u = %ld\n", b_u);
  printf("b_v = %ld\n", b_v);
  printf("p_u = %ld\n", p_u);
  printf("p_v = %ld\n", p_v);
  printf("b_u another way = %ld\n", c_u);
  printf("b_v another way = %ld\n", c_v);
  printf("K_uv = %ld\n", K_uv);
  printf("K_vu = %ld\n", K_vu);
  return 0;
}