/*
  Author:  Pate Williams (c) 1997

  Exercise
  "5.8 Let E be the elliptic curve
  y ^ 2 = x ^ 3 + x + 28 defined over Z_71.
  (a) Determine the number of points on E.
  (b) Show that E is not a cyclic group.
  (c) What is the maximum order of an element
      in E? Find an element having this order."
  -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 142l

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

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 = 1, b = 28, h, i, p = 71;
  long max_order, min_order, noncyclic, ord;
  struct point e[MAX_POINTS], P, Q;

  printf("E: y ^ 2 = x ^ 3 + %ld * x + %ld mod %ld\n",
         a, b, p);
  printf("the points (x, y) on E are:\n");
  h = E_points(a, b, p, e) + 1;
  printf("the number of points on E is: %ld\n", h);
  max_order = 0;
  min_order = h;
  printf("the points (x, y) on E and their orders are:\n");
  for (i = 0; i < h - 1; i++) {
    ord = order(a, p, e[i]);
    printf("(%2ld, %2ld) %2ld ", e[i].x, e[i].y, ord);
    if ((i + 1) % 3 == 0) printf("\n");
    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("\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);
  return 0;
}