/*
  Author:  Pate Williams (c) 1997

  Exercise
  "4.8 This exercise illustrates another example of
  a protocol failure (due to Simmons) involving RSA;
  it is called the common modulus protocol failure.
  Suppose Bob has an RSA Cryptosystem with modulus
  n and encryption exponent b1, and Charlie has an
  RSA Cryptosystem with (the same) modulus n and
  encryption exponent b2. Suppose that
  gcd(b1, b2) = 1. Now, consisder the situation
  that arises if Alice encrypts the same plaintext
  x to send to both Bob and Charlie. Thus, she
  computes y1 = x ^ b1 mod n and y2 = x ^ b2 mod n,
  and sends y1 to Bob and y2 to Charlie. Suppose
  Oscar intercepts y1 and y2, and performs the
  following computations:
  1. compute c1 = b1 ^ - 1 mod b2
  2. compute c2 = (c1 * b1 - 1) / b2
  3. compute x1 = y1 ^ c1 * (y2 ^ c2) ^ - 1 mod n
  (b) Illustrate this attack by computing x by
      this method if n = 18721, b1 = 945,
      b2 = 7717, y1 = 10510, and y2 = 14702."
  -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by
  Douglas R. Stinson page 159.
*/

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

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 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 gcd(long a, long b)
/* Euclid's algorithm */
{
  long r;

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

int main(void)
{
  long b1 = 945l, b2 = 7717l, n = 18721l;
  long y1 = 10510l, y2 = 14702l;
  long c1, c2, x, x1, x2;

  if (gcd(b1, b2) != 1) {
    fprintf(stderr, "fatal error\n");
    fprintf(stderr, "gcd(b1, b2) != 1\n");
    exit(1);
  }
  c1 = Extended_Euclidean(b1, b2);
  c2 = (c1 * b1 - 1) / b2;
  x1 = exp_mod(y1, c1, n);
  x2 = exp_mod(y2, c2, n);
  x2 = Extended_Euclidean(x2, n);
  x = (x1 * x2) % n;
  x1 = exp_mod(x, b1, n);
  x2 = exp_mod(x, b2, n);
  printf("b1 = %5ld b2 = %5ld n = %ld\n", b1, b2, n);
  printf("c1 = %5ld c2 = %5ld\n", c1, c2);
  printf("y1 = %5ld y2 = %5ld x = %ld\n", y1, y2, x);
  printf("y1 = %5ld y2 = %5ld\n", x1, x2);
  if (x1 != y1 || x2 != y2) {
    fprintf(stderr, "fatal error\n");
    fprintf(stderr, "can't compute y1 and y2 from x\n");
    fprintf(stderr, "y1 %5ld y2 %5ld\n", x1, x2);
    exit(1);
  }
  return 0;
}