/*
  Author:  Pate Williams (c) 1997

  Exercise
  "5.9 Let E be the elliptic curve
  y ^ 2 = x ^ 3 + x + 13 defined over Z_31. It
  can be shown that #E = 34 and (9, 10) is an
  element of order 34 in E. The Menezes-Vanstone
  Cryptosystem defined on E will have as its
  plaintext space Z_34* * Z_34*. Suppose Bob's
  secret exponent is a = 25.
  (a) Compute beta = a alpha.
  (b) Decrypt the following string of ciphertext:
      ((4, 9), 28, 7), ((19, 28), 9, 13),
      ((5, 22), 20, 17), ((25, 16), 12, 27).
  (c) Assuming that each plaintext represents two
      alphabetic characters, convert the plaintext
      into an English word. (Here we will use the
      correspondence A = 1,..., Z = 26, since 0
      is not allowed in a (plaintext) ordered
      pair." -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by
  Douglas R. Stinson page 201.
*/

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

#define MAX_POINTS 64l

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

long E_points(long a, long b, long p,
              struct point *e)
/* returns the number of points on the elliptic
   curve y ^ 2 = x ^ 3 + ax + b mod p */
{
  long count = 0, m = (p + 1) / 4, x, y;

  if (p % 4 == 3) {
    for (x = 0; x < p; x++) {
      y = x * x * x + a * x + b;
      if (JACOBI(y, p) != - 1) {
        y = exp_mod(y, m, p);
        printf("(%2ld, %2ld) ", x, y);
        e[count].x = x;
        e[count].y = y;
        y = - y % p;
        if (y < 0) y += p;
        printf("(%2ld, %2ld) ", x, y);
        e[count + 1].x = x;
        e[count + 1].y = y;
        count += 2;
        if (count % 4 == 0) printf("\n");
      }
    }
    if (count % 4 != 0) printf("\n");
  }
  return count;
}

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

  if (P.x == Q.x && P.y == 0 && Q.y == 0) {
    R->x = 0;
    R->y = 1;
    return;
  }
  if (P.x == Q.x && P.y == p - Q.y) {
    R->x = 0;
    R->y = 1;
    return;
  }
  if (P.x == 0 && P.y == 1) {
    *R = Q;
    return;
  }
  if (Q.x == 0 && Q.y == 1) {
    *R = P;
    return;
  }
  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;
}

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

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

long order(long a, long p, struct point P)
{
  long ord = 1;
  struct point Q = P, R;

  do {
    ord++;
    add(a, p, P, Q, &R);
    Q = R;
  } while (R.x != 0 && R.y != 1);
  return ord;
}

int main(void)
{
  long A =  25, a = 1, b = 13, h, i, p = 31;
  long max_order, min_order, noncyclic, ord;
  long c1, c2, m1, m2;
  struct ciphertext c[4] ={{{4, 9}, 28, 7},
                           {{19, 28}, 9, 13},
                           {{5, 22}, 20, 17},
                           {{25, 16}, 12, 27}};
  struct point alpha = {9, 10}, beta;
  struct point e[MAX_POINTS], P, Q;

  printf("E: y ^ 2 = x ^ 3 + %ld * x + %ld mod %ld\n",
         a, b, p);
  printf("the the non-origin points (x, y) on E are:\n");
  h = E_points(a, b, p, e) + 1;
  printf("the number of non-origin points on E is: %ld\n", h - 1);
  max_order = 0;
  min_order = h + 1;
  printf("the non-origin points (x, y) on E ");
  printf("and their orders are:\n");
  for (i = 0; i < h - 1; i++) {
    ord = order(a, p, e[i]);
    if (ord < h) noncyclic = 1;
    if (ord < min_order) {
      P = e[i];
      min_order = ord;
    }
    if (ord > max_order) {
      Q = e[i];
      max_order = ord;
    }
    printf("(%2ld, %2ld) %2ld ", e[i].x, e[i].y, ord);
    if ((i + 1) % 3 == 0) printf("\n");
  }
  printf("\n");
  if (noncyclic) printf("E is not a cyclic group\n");
  printf("minimum order: %ld\n", min_order);
  printf("maximum order: %ld\n", max_order);
  printf("minimum (%2ld, %2ld)\n", P.x, P.y);
  printf("maximum (%2ld, %2ld)\n", Q.x, Q.y);
  printf("alpha = (%2ld, %2ld)\n", alpha.x, alpha.y);
  multiply(a, A, p, alpha, &beta);
  printf("beta = %ld alpha = (%2ld, %2ld)\n",
          A, beta.x, beta.y);
  printf("the ciphertext is:\n");
  for (i = 0; i < 4; i++)
    printf("((%2ld, %2ld), %2ld, %2ld)\n",
           c[i].y0.x, c[i].y0.y, c[i].y1, c[i].y2);
  printf("the corresponding plaintext is:\n");
  for (i = 0; i < 4; i++) {
    multiply(a, A, p, c[i].y0, &beta);
    c1 = beta.x;
    c2 = beta.y;
    m1 = (c[i].y1 * Extended_Euclidean(c1, p)) % p;
    m2 = (c[i].y2 * Extended_Euclidean(c2, p)) % p;
    printf("%c%c", m1 + 'A' - 1, m2 + 'A' - 1);
  }
  printf("\n");
  return 0;
}