/*
  Author:  Pate Williams (c) 1997

  Exercise VI.2.11
  "Use the elliptic curve analog of ElGamal to send
  the message in Exercise 3(a) with E and p as in
  Exercise 3 and B = (0, 0). Suppose that your
  correspondent's public key is the point (201, 380)
  and your sequence of k's (one used to send each
  message unit) is 386, 209, 118, 589, 312, 483,
  335. What sequence of 7 pairs of points do you
  send?" -Neal Koblitz-
  See "A Course in Number Theory and Cryptography"
  by Neal Koblitz second edition page 186.
*/

#include <stdio.h>
#include <string.h>

struct point {long x, y;};
struct ciphertext {struct point y0; long y1, y2;};

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;
  }
  if (a < 0) a += 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;
}

struct point add(long a, long p, struct point P,
                 struct point Q)
/* elliptic curve point partial addition */
{
  long i, lambda;
  struct point R;

  if (P.x == Q.x && P.y == 0 && Q.y == 0) {
    R.x = 0;
    R.y = 1;
    return R;
  }
  if (P.x == Q.x && P.y == p - Q.y) {
    R.x = 0;
    R.y = 1;
    return R;
  }
  if (P.x == 0 && P.y == 1) return Q;
  if (Q.x == 0 && Q.y == 1) return P;
  if (P.x != Q.x) {
    i = Q.x - P.x;
    if (i < 0) i += p;
    i = Extended_Euclidean(i, p);
    lambda = ((Q.y - P.y) * i) % p;
  }
  else {
    i = Extended_Euclidean((2 * P.y) % p, p);
    lambda = ((3 * P.x * P.x + a) * i) % p;
  }
  if (lambda < 0) lambda += p;
  R.x = (lambda * lambda - P.x - Q.x) % p;
  R.y = (lambda * (P.x - R.x) - P.y) % p;
  if (R.x < 0) R.x += p;
  if (R.y < 0) R.y += p;
  return R;
}

struct point multiply(long a, long k, long p,
                      struct point P)
{
  struct point R, S;

  R.x = 0;
  R.y = 1;
  S = P;
  while (k != 0) {
    if (k & 1) R = add(a, p, R, S);
    k >>= 1;
    if (k != 0) S = add(a, p, S, S);
  }
  return R;
}

struct point plain_translate(char c)
{
  long a = - 1, b = 188, kappa = 20, m, p = 751;
  long x, y;
  struct point pt;

  if (c < '9') m = c - '0';
  else m = c - 'A' + 10;
  x = kappa * m;
  m = (p + 1) / 4;
  for (;;) {
    x++;
    if (x < 0) x += p;
    y = (x * x * x + a * x + b) % p;
    if (JACOBI(y, p) != - 1) {
      y = exp_mod(y, m, p);
      /*y = (y - 376) % p;*/
      if (y < 0) y += p;
      pt.x = x, pt.y = y;
      return pt;
    }
  }
}

int main(void)
{
  char plaintext[8] = "STOP007";
  long a = - 1, p = 751;
  long i, k[7] = {386, 209, 118, 589, 312, 483, 335};
  struct point P, P_m, Q, B = {0, 376}, a_B = {201, 380};

  a_B.y = (a_B.y + 376) % p;
  for (i = 0; i < strlen(plaintext); i++) {
    P_m = plain_translate(plaintext[i]);
    P = multiply(a, k[i], p, B);
    Q = multiply(a, k[i], p, a_B);
    Q = add(a, p, P_m, Q);
    P_m.y = (P_m.y - 376) % p;
    if (P_m.y < 0) P_m.y += p;
    P.y = (P.y - 376) % p;
    if (P.y < 0) P.y += p;
    Q.y = (Q.y - 376) % p;
    if (Q.y < 0) Q.y += p;
    printf("(%3ld, %3ld), (%3ld, %3ld), (%3ld, %3ld)\n",
           P_m.x, P_m.y, P.x, P.y, Q.x, Q.y);
  }
  return 0;
}