/*
  Author:  Pate Williams (c) 1997

  Exercise VI.4.4
  "Each of the values of n in Exercise 5 of Section
  V.4 has a factor p < 100. In each case (a)-(k)
  find this factor by Lentra's elliptic curve
  method, choosing B = 5, C = 120, P = (1, 1),
  and E: y ^ 2 = x ^ 3 + a x - a with a = 1, 2,...
  (taking a's for which the discriminant is prime
  to n). In each case, what is the first value of
  a for which you find the factor, and what is the
  value of k_1 for which the factor appears as
  gcd(denominator, n) in your computation of
  k_1 P?" -Neal Koblitz-
  See "A Course in Number Theory and Cryptography"
  by Neal Koblitz second edition page 199.
*/

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

struct factor {long expon, prime;};
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 gcd(long a, long b)
{
  long r;

  while (b > 0)
    r = a % b, a = b, b = r;
  return a;
}

void Euclid_extended(long a, long b, long *u,
                     long *v, long *d)
{
  long q, t1, t3, v1, v3;

  *u = 1, *d = a;
  if (b == 0) {
    *v = 0;
    return;
  }
  v1 = 0, v3 = b;
  #ifdef DEBUG
  printf("----------------------------------\n");
  printf(" q    t3   *u   *d   t1   v1   v3\n");
  printf("----------------------------------\n");
  #endif
  while (v3 != 0) {
    q = *d / v3;
    t3 = *d - q * v3;
    t1 = *u - q * v1;
    *u = v1, *d = v3;
    #ifdef DEBUG
    printf("%4ld %4ld %4ld ", q, t3, *u);
    printf("%4ld %4ld %4ld %4ld\n", *d, t1, v1, v3);
    #endif
    v1 = t1, v3 = t3;
  }
  *v = (*d - a * *u) / b;
  #ifdef DEBUG
  printf("----------------------------------\n");
  #endif
}

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

  if (P.x == Q.x && P.y == 0 && Q.y == 0) {
    R->x = 0;
    R->y = 1;
    return 0;
  }
  if (P.x == Q.x && P.y == p - Q.y) {
    R->x = 0;
    R->y = 1;
    return 0;
  }
  if (P.x == 0 && P.y == 1) {
    *R = Q;
    return 0;
  }
  if (Q.x == 0 && Q.y == 1) {
    *R = P;
    return 0;
  }
  if (P.x != Q.x) {
    i = Q.x - P.x;
    if (i < 0) i += p;
    Euclid_extended(i, p, &u, &v, &d);
    if (d != 1) return d;
    lambda = ((Q.y - P.y) * u) % p;
  }
  else {
    Euclid_extended((2 * P.y) % p, p, &u, &v, &d);
    if (d != 1) return d;
    lambda = ((3 * P.x * P.x + a) * u) % 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 0;
}

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

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

int main(void)
{
  int found;
  long a, d, i, j, k, l, ni;
  long n[11] = {9509, 13561, 8777, 14429, 12403,
                14527, 10123, 12449, 9353, 25511, 17873};
  struct factor f[3] = {{6, 2}, {4, 3}, {2, 5}};
  struct point P = {1, 1}, Q;

  for (i = 0; i < 11; i++) {
    found = 0;
    a = 1;
    ni = n[i];
    while (!found) {
      for (;;) {
        d = 4 * a * a * a + 27 * a * a;
        if (gcd(d, ni) == 1) break;
        a++;
      }
      k = 1;
      for (j = 0; !found && j < 3; j++) {
        for (l = 0; !found && l < f[j].expon; l++) {
          k *= f[j].prime;
          d = multiply(a, k, ni, P, &Q);
          found = d != 0;
          if (found)
            printf("(%c) %ld %2ld %ld\n", i + 'a', a, d, k);
        }
      }
      a++;
    }
  }
  return 0;
}