/*
  Author:  Pate Williams (c) 1997

  "Algorithm 8.1.1 (Trial Division). We assume given
  a table of prime numbers p[1] = 2, p[2] = 3,..., p[k],
  with k > 3, an array t <- [6, 4, 2, 4, 4, 6, 2], and
  an index j such that if p[k] mod 30 is equal to
  1, 7, 11, 13, 17, 19, 23 or 29 then j is equal to
  0, 1, 2, 3, 4, 5, 6 or 7 respectively. Finally, we
  give ourselves an upper bound B such that B >= p[k],
  essentially to avoid spending too much time. Then
  given a positive integer N, this algorithm tries to
  factor (or split N) and if it fails, N will be free
  of prime factors less than or equal to B."
  -Henri Cohen-
  See "A Course in Computational Algebraic Number
  Theory" by Henri Cohen page 420.
*/

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

#define BITS_PER_LONG   32
#define BITS_PER_LONG_1 31
#define MAX_SIEVE 500000l
#define PRIME_SIZE 41538l
#define SIEVE_SIZE (MAX_SIEVE / BITS_PER_LONG + 1)

struct factor {long expon, prime;};

long get_bit(long i, long *sieve)
{
  long b = i % BITS_PER_LONG;
  long c = i / BITS_PER_LONG;

  return (sieve[c] >> (BITS_PER_LONG_1 - b)) & 1;
}

void set_bit(long i, long v, long *sieve)
{
  long b = i % BITS_PER_LONG;
  long c = i / BITS_PER_LONG;
  long mask = 1 << (BITS_PER_LONG_1 - b);

  if (v == 1)
    sieve[c] |= mask;
  else
    sieve[c] &= ~mask;
}

void Sieve(long n, long *sieve)
{
  long c, i, inc;

  set_bit(0, 0, sieve);
  set_bit(1, 0, sieve);
  set_bit(2, 1, sieve);
  for (i = 3; i <= n; i++)
    set_bit(i, i & 1, sieve);
  c = 3;
  do {
    i = c * c, inc = c + c;
    while (i <= n) {
      set_bit(i, 0, sieve);
      i += inc;
    }
    c += 2;
    while (!get_bit(c, sieve)) c++;
  } while (c * c <= n);
}

int trial_division(long *N, long *prime,
                   struct factor *f, long *count)
{
  int found;
  long B, d, e, i, j = 0, k, l, m, n, r;
  long t[8] = {6, 4, 2, 4, 2, 4, 6, 2};
  long table[8] = {1, 7, 11, 13, 17, 19, 23, 29};

  if (*N <= 5) {
    *count = 1;
    f[0].expon = 1;
    switch (*N) {
      case 1 : f[0].prime = 1; break;
      case 2 : f[0].prime = 2; break;
      case 3 : f[0].prime = 3; break;
      case 4 : f[0].expon = 2; f[0].prime = 2; break;
      case 5 : f[0].prime = 5; break;
    }
    return 1;
  }
  *count = 0;
  B = prime[PRIME_SIZE - 1];
  i = - 1, l = sqrt(*N), m = - 1;
  L2:
  m++;
  if (i == PRIME_SIZE) {
    i = j - 1;
    goto L5;
  }
  d = prime[m];
  k = d % 30;
  found = 0;
  for (n = 0; !found && n < 8; n++) {
    found = k == table[n];
    if (found) j = n;
  }
  L3:
  r = *N % d;
  if (r == 0) {
    e = 0;
    do {
      e++;
      *N /= d;
    } while (*N % d == 0);
    f[*count].expon = e;
    f[*count].prime = d;
    *count = *count + 1;
    if (*N == 1) return 1;
    goto L2;
  }
  if (d >= l) {
    if (*N > 1) {
      f[*count].expon = 1;
      f[*count].prime = *N;
      *count = *count + 1;
      return 1;
    }
  }
  else if (i < 0) goto L2;
  L5:
  i = (i + 1) % 8;
  d = d + t[i];
  if (d > B) return 0;
  goto L3;
}

int main(void)
{
  long N, count, i, p = 2;
  long prime[PRIME_SIZE], sieve[SIEVE_SIZE];
  struct factor f[32];

  Sieve(MAX_SIEVE, sieve);
  for (i = 0; i < PRIME_SIZE; i++) {
    while (!get_bit(p, sieve)) p++;
    prime[i] = p++;
  }
  for (;;) {
    printf("N = ");
    scanf("%ld", &N);
    if (N <= 0) break;
    trial_division(&N, prime, f, &count);
    printf("factors:\n");
    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");
    }
  }
  return 0;
}