/*
  Author:  Pate Williams (c) 1997

  Exercise VI.2.3
  "Let E be the elliptic curve
  y ^ 2 + y = x ^ 3 - x defined over the field of
  p = 751 elements. (A change of variables of the
  form y' = y + 376 will convert this equation to
  the form (1) of Section 1.) This curve contains
  N = 727 points. Suppose that the plaintext
  message units are the decimal digits 0-9 and the
  letters A-Z with numerical equivalents 10-35,
  respectively. Take kappa = 20.
  (a) Use the method in the text to write the
  message "STOP007" as a sequence of seven points
  on the curve.
  (b) Translate the sequence of points (361, 383),
  (241, 605), (201, 380), (461, 467), (581, 395)
  into a reply message." -Neal Koblitz-
  See "A Course in Number Theory and Cryptography"
  by Neal Koblitz second edition page 185.
*/

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

struct point {long x, y;};

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

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

char point_translate(struct point pt)
{
  long kappa = 20, m, x;

  for (x = pt.x; x >= 0; x--) {
    m = x / kappa;
    if (m <= 35) {
      if (m < 9) return (char) (m + '0');
      return (char) (m + 'A' - 10);
    }
  }
  return 0;
}

int main(void)
{
  char plaintext[8] = "STOP007";
  long a = - 1, b = 188, a_p = 0, i;
  struct point p[5] = {{361, 383}, {241, 605},
                       {201, 380}, {461, 467},
                       {581, 395}};
  struct point pt;

  for (i = 0; i < 751; i++)
    a_p += JACOBI((i * i * i + a * i + b) % 751, 751);
  printf("N = |E(F_751)| = %ld\n", 752 + a_p);
  for (i = 0; i < strlen(plaintext); i++) {
    pt = plain_translate(plaintext[i]);
    printf("(%3ld, %3ld)\n", pt.x, pt.y);
  }
  for (i = 0; i < 5; i++)
    printf("%c", point_translate(p[i]));
  printf("\n");
  return 0;
}