/*
  Author:  Pate Williams (c) 1997

  Exercise
  "4.18 Suppose p = 199, q = 211, and B = 1357
        in the Rabin Cryptosystem. Perform the
        following computations.
        (a) Determine the four square roots of 1
        modulo n, where n = pq.
        (b) Compute the encryption y = ek(32767).
        (c) Determine the four possible
            decryptions of this given ciphetext."
  -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by
  Douglas R. Stinson page 161.
*/

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

int JACOBI(long a, long n)
{
  int s;
  long a1, b = a, e = 0, m, n1;

  if (a == 0) return 0;
  if (a == 1) return 1;
  while ((b & 1) == 0)
    b >>= 1, e++;
  a1 = b;
  m = n % 8;
  if (!(e & 1)) s = 1;
  else if (m == 1 || m == 7) s = + 1;
  else if (m == 3 || m == 5) s = - 1;
  if (n % 4 == 3 && a1 % 4 == 3) s = - s;
  if (a1 != 1) n1 = n % a1; else n1 = 1;
  return s * JACOBI(n1, a1);
}

long exp_mod(long x, long b, long n)
/* returns x ^ b mod n */
{
  long a = 1l, s = x;

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

void extended_euclid(long a, long b, long *x,
                     long *y, long *d)
/* calculates a * *x + b * *y = gcd(a, b) = *d */
{
  long q, r, x1, x2, y1, y2;

  if (b == 0) {
    *d = a, *x = 1, *y = 0;
    return;
  }
  x2 = 1, x1 = 0, y2 = 0, y1 = 1;
  while (b > 0) {
    q = a / b, r = a - q * b;
    *x = x2 - q * x1, *y = y2 - q * y1;
    a = b, b = r;
    x2 = x1, x1 = *x, y2 = y1, y1 = *y;
  }
  *d = a, *x = x2, *y = y2;
}

long inverse(long a, long n)
/* computes the inverse of a modulo n */
{
  long d, x, y;

  extended_euclid(a, n, &x, &y, &d);
  if (x < 0) x += n;
  if (d == 1) return x;
  return 0;
}

void sqrt_mod(long a, long p, long *r)
{
  long ai, b, c, d, i, p1 = p - 1, s = 0, t = p1;
  long x, y;

  if (JACOBI(a, p) == - 1) {
    *r = 0;
    return;
  }
  do
    do b = rand() % p; while (b == 0);
  while (JACOBI(b, p) != - 1);
  while (!(t & 1)) t >>= 1, s++;
  ai = inverse(a, p);
  if ((ai * a) % p != 1) {
    printf("*error*\nin inverse\n");
    printf("a = %ld a ^ - 1 = %ld\n", a, ai);
  }
  c = exp_mod(b, t, p);
  *r = exp_mod(a, (t + 1) >> 1, p);
  x = pow(2l, s);
  for (i = 1; i < s; i++) {
    y = (((*r * *r) % p) * ai) % p;
    d = exp_mod(y, x, p);
    if (d < 0) d += p;
    if (d == p1) *r = (*r * c) % p;
    c = (c * c) % p;
    x *= 2;
  }
  if (*r < 0) *r += p;
}

void square_roots(long a, long p, long q,
                  long *x1, long *x2,
                  long *y1, long *y2)
{
  long b, c, d, e, n = p * q, r, s, x, y;

  if (a < p) b = a; else b = a % p;
  sqrt_mod(b, p, &r);
  if ((r * r) % p != b) {
    printf("*error*\nin sqrt_mod p\n");
    printf("(a / p) = %d\n", JACOBI(a, p));
    printf("%ld * %ld %% %ld = %ld\n", r, r,
           p, (r * r) % p);
  }
  if (a < q) b = a; else b = a % q;
  sqrt_mod(b, q, &s);
  if ((s * s) % q != b) {
    printf("*error*\nin sqrt_mod q\n");
    printf("(a / q) = %d\n", JACOBI(a, q));
    printf("%ld * %ld %% %ld = %ld\n", s, s,
           q, (s * s) % q);
  }
  extended_euclid(p, q, &c, &d, &e);
  x = (r * d * q + s * c * p) % n;
  y = (r * d * q - s * c * p) % n;
  *x1 = x;
  *x2 = - x % n;
  if (*x1 < 0) *x1 += n;
  if (*x2 < 0) *x2 += n;
  *y1 = y;
  *y2 = - y % n;
  if (*y1 < 0) *y1 += n;
  if (*y2 < 0) *y2 += n;
}

int main(void)
{
  long B = 1357l, B2, B4;
  long a = 1l, p = 199l, q = 211l;
  long i, n = p * q, s[4];
  long x = 32767l, y, z;

  square_roots(a, p, q, &s[0], &s[1], &s[2], &s[3]);
  printf("the four square roots of %ld mod %ld are:\n", a, n);
  for (i = 0; i < 4; i++) {
    printf("%5ld ", s[i]);
    if ((s[i] * s[i]) % n != a)
      printf("*error*\nin square_roots\n");
  }
  B2 = (B * inverse(2l, n)) % n;
  B4 = (((B * B) % n) * inverse(4l, n)) % n;
  y = (x * ((x + B) % n)) % n;
  printf("\nek(%5ld) = %5ld\n", x, y);
  z = (B4 + y) % n;
  square_roots(z, p, q, &s[0], &s[1], &s[2], &s[3]);
  printf("the four possible plaintexts are:\n");
  for (i = 0; i < 4; i++) {
    z = (s[i] - B2) % n;
    if (z < 0) z += n;
    printf("%5ld ", z);
    if ((z * (z + B)) % n != y)
      printf("*error*\nin square_roots\n");
  }
  return 0;
}