/*
  Author:  Pate Williams (c) 1997

  Exercise V.4.5
  "Following Examples 2 and 3, use the continued
  fraction algorithm to factor the following
  numbers: (a) 9509; (b) 13561; (c) 8777;
  (d) 14429; (e) 12403; (f) 14527; (g) 10123;
  (h) 12449; (i) 9353; (j) 25511; (k) 17873."
  -Neal Koblitz-
  See "A Course in Number Theory and Cryptography"
  by Neal koblitz second edition pages 159-160.
*/

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

long gcd(long a, long b)
{
  long r;

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

int trial_divide(long n, long Bn, long *B, long *vector)
{
  long e, i, imin, p;

  if (n < 0) {
    if (B[0] != - 1) return 0;
    vector[0] = 1;
    n = - n;
    imin = 1;
  }
  else if (B[0] == - 1) {
    vector[0] = 0;
    imin = 1;
  }
  else imin = 0;
  for (i = imin; i < Bn; i++) {
    p = B[i];
    if (n % p == 0) {
      e = 0;
      do {
        e++;
        n /= p;
      } while (n % p == 0);
      vector[i] = e;
    }
    else vector[i] = 0;
  }
  return n == 1;
}

int routine(long bi, long n, long Bn, long *B,
            long *epsilon, long *b)
{
  int found, value;
  long A = (bi * bi) % n, c = 1, d, e, i;
  long vector[32];

  value = trial_divide(A, Bn, B, vector);
  if (!value) {
    A = A - n;
    value = trial_divide(A, Bn, B, vector);
  }
  if (!value) return 0;
  *b *= bi;
  *b %= n;
  for (i = 0; i < Bn; i++)
    epsilon[i] += vector[i];
  found = (epsilon[0] & 1) == 0;
  for (i = 1; found && i < Bn; i++)
    found = (epsilon[i] & 1) == 0;
  if (!found) return 0;
  for (i = 0; i < Bn; i++)
    c *= pow(B[i], epsilon[i] / 2);
  if (c < 0) c = - c;
  c %= n;
  d = gcd(*b + c, n);
  if (d == 1 || d == n) return 0;
  e = n / d;
  if (d > e) c = d, d = e, e = c;
  printf("b = %5ld c = %5ld ", *b, c);
  printf("n = %5ld = %ld * %ld\n", n, d, e);
  return 1;
}

void factor(long n, long imax, long Bn, long *B)
{
  double x0, x1;
  long a0, a1, b = 1, b0, b1, bm2, i;
  long epsilon[12] = {0};

  bm2 = 1;
  a0 = b0 = sqrt(n);
  x0 = sqrt(n) - a0;
  routine(b0, n, Bn, B, epsilon, &b);
  for (i = 0; i < imax; i++) {
    a1 = 1.0 / x0;
    x1 = 1.0 / x0 - a1;
    b1 = (a1 * b0 + bm2) % n;
    if (routine(b1, n, Bn, B, epsilon, &b)) return;
    bm2 = b0;
    b0 = b1;
    x0 = x1;
  }
}

int main(void)
{
  long i;
  long Bn[11] = {4, 3, 2, 3, 3, 4, 6, 3, 4, 5, 5};
  long B[11][6] = {{- 1, 2, 5, 11}, {2, 3, 5}, {- 1, 2},
                   {- 1, 2, 29}, {- 1, 3, 13},
                   {- 1, 2, 3, 7}, {- 1, 2, 3, 7, 11, 13},
                   {- 1, 2, 5}, {- 1, 2, 7, 11},
                   {- 1, 2, 5, 23, 29}, {- 1, 2, 7, 11, 23}};
  long imax[11] = {4, 4, 3, 3, 7, 6, 10, 9, 10, 9};
  long n[11] = {9509, 13561, 8777, 14429, 12403,
                14527, 10123, 12449, 9353, 25511, 17873};

  for (i = 0; i < 11; i++)
    factor(n[i], imax[i], Bn[i], B[i]);
  return 0;
}