/*
  Author:  Pate Williams (c) 1997

  "Algorithm 8.7.2 (Shanks's SQUFOF). Given an odd
  integer N, this algorithm tries to find a
  non-trivial factor of N." -Henri Cohen-
  See "A Course in Computational Algebraic Number
  Theory" by Henri Cohen page 437.
*/

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

struct factor {long expon, prime;};

long exp_mod(long x, long b, long n)
/* returns x ^ b mod n */
{
  long a = 1l, s = x;

  while (b != 0) {
    if (b & 1l) a = (a * s) % n;
    b >>= 1;
    if (b != 0) s = (s * s) % n;
  }
  return a;
}

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

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

int RabinMiller(long N)
/* returns 0 if N is composite 1 otherwise */
{
  long N1 = N - 1, a, b, c = 20, e, q = N1, t = 0;

  if (N == 2) return 1;
  while (!(q & 1)) q >>= 1, t++;
  while (c > 0) {
    do a = rand() % N; while (a == 0);
    e = 0;
    b = exp_mod(a, q, N);
    if (b != 1) {
      while (b != 1 && b != N1 && e <= t - 2) {
        b = (b * b) % N;
         e++;
      }
      if (b != N1) return 0;
    }
    c--;
  }
  return 1;
}

int square_test(long n, long *s)
{
  int square = 1;
  long q11[11], q63[63], q64[64], q65[65];
  long k, r, t;

  for (k = 0; k < 11; k++) q11[k] = 0;
  for (k = 0; k < 6; k++) q11[(k * k) % 11] = 1;
  for (k = 0; k < 63; k++) q63[k] = 0;
  for (k = 0; k < 32; k++) q63[(k * k) % 63] = 1;
  for (k = 0; k < 64; k++) q64[k] = 0;
  for (k = 0; k < 32; k++) q64[(k * k) % 64] = 1;
  for (k = 0; k < 65; k++) q65[k] = 0;
  for (k = 0; k < 33; k++) q65[(k * k) % 65] = 1;
  t = n % 64;
  r = n % 45045l;
  if (q64[t] == 0) square = 0;
  if (square && q63[r % 63] == 0) square = 0;
  if (square && q65[r % 65] == 0) square = 0;
  if (square && q11[r % 11] == 0) square = 0;
  if (square) *s = sqrt(n); else *s = 0;
  return square;
}

void reduce(long A, long B, long C,
            long *a, long *b, long *c)
{
  long D, l, r, s;

  *a = A;
  *b = B;
  *c = C;
  D = B * B - 4 * A * C;
  s = sqrt(D);
  l = labs(*c);
  if (l < s / 2)
    r = - *b + 2 * l * ((*b + s) / (2 * l));
  else
    r = - *b + 2 * l;
  *a = *c;
  *b = r;
  *c = (r * r - D) / (4 * *c);
}

void reduction(long A, long B, long C,
               long *a, long *b, long *c)
{
  long D, l, r, s;

  *a = A;
  *b = B;
  *c = C;
  D = B * B - 4 * A * C;
  s = sqrt(D);
  L1:
  if (*b < 0) *b = - *b, *a = - *a, *c = - *c;
  if (*b > labs(s - 2 * labs(*a)) && *b <= s) return;
  l = labs(*c);
  if (l < s / 2)
    r = - *b + 2 * l * ((*b + s) / (2 * l));
  else
    r = - *b + 2 * l;
  *a = *c;
  *b = r;
  *c = (r * r - D) / (4 * *c);
  goto L1;
}

int SQUFOF(long *N, long *f, long *e)
{
  int found;
  long A, B, C, D, L, Q[1024], Q_size, a, b, b1, c;
  long ap, bp, cp, d, i, j, s, z;

  if (RabinMiller(*N)) return - 2;
  if (square_test(*N, f)) {
    *e = 2;
    *N = 1;
    return - 1;
  }
  if (*N % 4 == 1) {
    D = *N;
    d = sqrt(*N);
    b = (d - 1) / 2;
    b = (b << 1) + 1;
  }
  else {
    D = 4 * *N;
    d = sqrt(D);
    b = 2 * (d / 2);
  }
  a = 1;
  c = (b * b - D) / 4;
  Q_size = 0;
  i = 0;
  L = sqrt(d);
  ap = a;
  bp = b;
  cp = c;
  L4:
  reduce(ap, bp, cp, &A, &B, &C);
  i++;
  if (i & 1) goto L7;
  if (square_test(labs(A), &z)) {
    a = z;
    found = 0;
    for (j = 0; !found && j < Q_size; j++)
      found = a == Q[j];
    if (!found) goto L8;
  }
  if (A == 1) return 0;
  L7:
  if (labs(A) <= L) {
    Q[Q_size] = A;
    Q_size++;
  }
  ap = A;
  bp = B;
  cp = C;
  goto L4;
  L8:
  s = gcd(a, gcd(B, D));
  if (s > 1) {
    *f = s;
    *e = 2;
    *N /= s * s;
    return 1;
  }
  else
    reduction(a, - B, a * C, &a, &b, &c);
  L9:
  b1 = b;
  ap = a;
  bp = b;
  cp = c;
  reduce(a, b, c, &a, &b, &c);
  if (a == 1) {
    a = ap;
    b = bp;
    c = cp;
    goto L10;
  }
  if (b1 != b) goto L9;
  L10:
  if (a == 1) return 0;
  if (a & 1) {
    *e = 1;
    *f = labs(a);
    if (*f == 1 || *f == *N) return 0;
    *N /= *f;
    return 1;
  }
  *e = 1;
  *f = labs(a / 2);
  if (*f == 1 || *f == *N) return 0;
  *N /= *f;
  return 1;
}

int main(void)
{
  int flag;
  long N, count, i, j, n;
  struct factor f[32], t;

  for (;;) {
    printf("N = ");
    scanf("%ld", &N);
    if (N <= 0) break;
    if (N & 1) {
      n = N;
      flag = SQUFOF(&N, &f[0].prime, &f[0].expon);
      if (flag == - 2) printf("%ld is prime\n", N);
      else if (flag == - 1)
        printf("%ld = %ld * %ld\n", f[0].prime * f[0].prime,
               f[0].prime, f[0].prime);
      else if (flag == 0) printf("%ld not factored\n", N);
      else {
        count = 1;
        printf("%ld is composite\nfactors:\n", n);
        for (;;) {
          flag = SQUFOF(&N, &f[count].prime, &f[count].expon);
          if (flag == - 1) break;
          if (flag == 1) count++;
          if (flag == - 2) {
            f[count].expon = 1;
            f[count].prime = N;
            count++;
            break;
          }
          if (flag == 0) break;
          if (N == 1) break;
        }
        if (flag != 0) {
          for (i = 0; i < count - 1; i++)
            for (j = i + 1; j < count; j++)
              if (f[i].prime > f[j].prime)
                t = f[i], f[i] = f[j], f[j] = t;
          for (i = 0; i < count; i++) {
            printf("%ld", f[i].prime);
            if (f[i].expon > 1) printf(" ^ %ld\n", f[i].expon);
            else printf("\n");
          }
        }
        else
          printf("error - can't factor number\n");
      }
    }
    else
      printf("%ld must be odd!\n", N);
  }
  return 0;
}