/*
  Author:  Pate Williams (c) 1997

  Exercise V.3.5.
  "Let n = 2701. Use the B-numbers 52 ^ 2, 53 ^ 2
  mod n for a suitable factor-base B to factor
  2701. What are the epsilon's corresponding to
  52 and 53?" -Neal Koblitz-
  See "A Course in Number Theory and Cryptography"
  by Neal Koblitz second edition page 153.
*/

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

struct factor {long expon, prime;};

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

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

void trial_divide(long n, long Bn, long *B,
                  struct factor *f, long *count)
{
  long e, i, p;

  *count = 0;
  for (i = 0; i < Bn; i++) {
    p = B[i];
    if (n % p == 0) {
      e = 0;
      do {
        e++;
        n /= p;
      } while (n % p == 0);
      f[*count].expon = e;
      f[*count].prime = p;
      *count = *count + 1;
      if (n == 1) return;
    }
  }
}

void factor(long n, long Bnn, long *Bnumber)
{
  int found;
  long Bn = 7, a, b = 1, bi, c = 1, i, j, k, p, s, t;
  long count[32], epsilon[7] = {0};
  long B[7] = {2, 3, 5, 7, 11, 13, 17};
  struct factor f[32][7];

  for (i = 0; i < Bnn; i++) {
    bi = Bnumber[i];
    b *= bi;
    bi = (bi * bi) % n;
    trial_divide(bi, Bn, B, f[i], &count[i]);
    printf("%ld\n", bi);
    for (j = 0; j < count[i]; j++) {
      a = f[i][j].expon;
      p = f[i][j].prime;
      printf("%ld", p);
      if (a != 1) printf(" ^ %ld\n", a);
      else printf("\n");
      found = 0;
      for (k = 0; !found && k < Bn; k++) {
        found = p == B[k];
        if (found) epsilon[k] += a;
      }
    }
  }
  b %= n;
  for (i = 0; i < Bn; i++)
    c *= pow(B[i], epsilon[i] / 2);
  c %= n;
  s = gcd(b + c, n);
  t = n / s;
  if (s > t)
    b = s, s = t, t = b;
  printf("b = %ld c = %ld\n", b, c);
  printf("n = %ld = %ld * %ld\n", n, s, t);
}

int main(void)
{
  long Bnn = 2, Bnumber[2] = {52, 53}, n = 2701;

  factor(n, Bnn, Bnumber);
  return 0;
}