/*
  Author:  Pate Williams (c) 1997

  Single precision Goldwasser-Kilian primality test.
  See "A Course in Computational Algebraic Number
  Theory" by Henri Cohen Algorithm 9.2.3 page 470.
*/

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

#define BITS_PER_LONG   32
#define BITS_PER_LONG_1 31
#define SIEVE_LIMIT 128
#define SIEVE_SIZE 128

struct point {long x, y;};

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

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

void set_bit(long i, long v, long *sieve)
{
  int b = (int)(i % BITS_PER_LONG);
  int c = (int)(i / BITS_PER_LONG);
  long mask = 1L << (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);
}

long gcd(long a, long b)
/* returns greatest common divisor of a and b */
{
  long r;

  if (a < 0) a = - a;
  if (b < 0) b = - b;
  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 inverse(long a, long b)
/* returns inverse of a modulo b or 0 if it does not exist */
{
  long d, u, v;

  Euclid_extended(a, b, &u, &v, &d);
  if (d == 1) return u;
  return 0;
}

int addition_1(long n,
               struct point P1,
               struct point P2,
               struct point *P3)
/* returns 1 if P1 = - P2 therefore P1 + P2 = O, 0 otherwise */
/* P1 != P2 */
{
  long delta_x = (P2.x - P1.x) % n;
  long delta_y = (P2.y - P1.y) % n;
  long m;

  if (P1.x == P2.x && P1.y == - P2.y) {
    P3->x = 0, P3->y = 1;
    return 1;
  }
  /* calculate m = (y2 - y1)(x2 - x1) ^ - 1 mod n */
  if (delta_x < 0) delta_x += n;
  if (delta_y < 0) delta_y += n;
  m = delta_y * inverse(delta_x, n) % n;
  /* calculate x3 = m ^ 2 - (x1 + x2) mod n */
  P3->x = (m * m - (P1.x + P2.x)) % n;
  /* calculate y3 = m(x1 - x3) - y1 mod n */
  P3->y = (m * (P1.x - P3->x) - P1.y) % n;
  return 0;
}

void addition_2(long a,
                long n,
                struct point P1,
                struct point *P3)
/* P1 == P2 */
{
  long m;

  /* calculate m = (3x1 ^ 2 + a)(2y1) ^ - 1 mod n */
  m = (3 * P1.x * P1.x + a) * inverse(2 * P1.y, n) % n;
  /* calculate x3 = m ^ 2 - 2x1 mod n */
  P3->x = (m * m - 2 * P1.x) % n;
  /* calculate y3 = m(x1 - x3) - y1 mod n */
  P3->y = (m * (P1.x - P3->x) - P1.y) % n;
}

int multiply(long a,
             long k,
             long n,
             struct point P,
             struct point *R,
             long *d)
/* returns - 1 if O encountered, 0 if divisor not found, 1 otherwise */
{
  int value = 1;
  struct point A, B, C;

  /*  A = P */
  A = P;
  /* B = O = (0, 1) the point at infinity */
  B.x = 0, B.y = 1;
  while (value && k > 0) {
    if (k & 1) {
      *d = gcd((B.x - A.x) % n, n);
      k--;
      value = *d == 1 || *d == n;
      if (A.x == 0 && A.y == 1);
      else if (B.x == 0 && B.y == 1) B = A;
      else if (value) {
        addition_1(n, A, B, &C);
        B = C;
      }
    }
    else {
      *d = gcd((2 * A.y) % n, n);
      k >>= 1;
      value = *d == 1 || *d == n;
      if (value) {
        addition_2(a, n, A, &C);
        A = C;
      }
    }
  }
  *R = B;
  if (R->x < 0) R->x += n;
  if (R->y < 0) R->y += n;
  if (R->x == 0 && R->y == 1) return - 1;
  return !value;
}

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;

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

int Rabin_Miller(long N)
/* given an integer N >= 3 returns 0 composite 1 probably prime */
{
  long N1 = N - 1, a, b, c = 20, e, q = N1, t = 0;

  if (N == 2 || N == 3) return 1;
  /* N - 1 = 2 ^ t * q with q odd */
  while (!(q & 1)) q >>= 1, t++;
  L2:
  do a = rand() % N; while (a == 0);
  e = 0;
  b = exp_mod(a, q, N);
  if (b == 1) goto L4;
  L3:
  while (b != 1 && b != N1 && e <= t - 2) {
    b = (b * b) % N;
    e++;
  }
  if (b != N1) return 0;
  L4:
  c--;
  if (c > 0) goto L2;
  return 1;
}

long cardinality(long a, long b, long p)
/* returns the cardinality of the elliptic curve
   E(F_p): y ^ 2 = x ^ 3 + ax + b using
   |E(F_p)| = p + 1 - a_p
   a_p = - sum(0 <= x < p, (y ^ 2 / p)) */
{
  long a_p = 0, x;

  for (x = 0; x < p; x++)
    a_p += JACOBI((((x * x) % p * x) % p + a * x + b) % p, p);
  return p + 1 + a_p;
}

long square_root_mod(long a, long p)
/* returns square root of a modulo p if it exists 0 otherwise */
{
  long b, e = 0, m, n, q = p - 1, r, t, x, y, z;

  /* p - 1 = 2 ^ e * q with q odd */
  while (!(q & 1)) q >>= 1, e++;
  /* find generator */
  do
    do n = rand() % p; while (p == 0);
  while (JACOBI(n, p) != - 1);
  z = exp_mod(n, q, p);
  y = z, r = e;
  x = exp_mod(a, (q - 1) / 2, p);
  b = ((a * x % p) * x) % p;
  x = (a * x) % p;
  L3:
  if (b == 1) return x;
  m = 1;
  while (exp_mod(b, pow(2, m), p) != 1) m++;
  if (m == r) return 0;
  t = exp_mod(y, pow(2, r - m - 1), p);
  y = t * t % p;
  r = m;
  x = x * t % p;
  b = b * y % p;
  goto L3;
}

int Goldwasser_Kilian(long N)
/* returns 0 if N is composite 1 if prime */
{
  int value;
  long Ni = N, a, b, d, i = 0, m, p, q, x, y;
  long sieve[SIEVE_SIZE];
  struct point P, P1, P2;

  Sieve(SIEVE_LIMIT, sieve);
  L2:
  L3:
  if (Ni <= SIEVE_LIMIT) {
    p = 2;
    while (p <= sqrt(Ni)) {
      if (Ni % p == 0) {
        printf("1 factor = %ld\n", p);
        return 0;
      }
      p++;
      while (!get_bit(p, sieve)) p++;
    }
    return 1;
  }
  if (!Rabin_Miller(Ni)) goto L9;
  for (;;) {
    do a = rand() % Ni; while (a == 0);
    do b = rand() % Ni; while (b == 0);
    m = (4 * a * a * a + 27 * b * b) % Ni;
    if (m != 0) break;
  }
  m = cardinality(a, b, Ni);
  q = m / 2;
  x = sqrt(sqrt(Ni)) + 1;
  if (m & 1 || !Rabin_Miller(q) || q <= x * x) goto L3;
  L6:
  do {
    do x = rand() % Ni; while (x == 0);
    y = (((x * x) % Ni * x) % Ni + a * x + b) % Ni;
    if (y < 0) y += p;
  } while (y == 0 || JACOBI(y, Ni) == - 1);
  y = square_root_mod(y, Ni);
  if (y == 0) goto L9;
  P.x = x, P.y = y;
  value = multiply(a, m, Ni, P, &P1, &d);
  if (value == 1) {
    printf("2 factor = %ld\n", d);
    goto L9;
  }
  if (value == 0) goto L6;
  value = multiply(a, 2, Ni, P, &P2, &d);
  if (value == 1) {
    printf("3 factor = %ld\n", d);
    goto L9;
  }
  if (value == - 1) goto L6;
  printf("N[%ld] = %ld\n", i, Ni);
  printf("a = %ld\n", a);
  printf("b = %ld\n", b);
  printf("m = %ld\n", m);
  printf("q = %ld\n", q);
  printf("P = (%ld, %ld)\n", P.x, P.y);
  printf("P1 = (%ld, %ld)\n", P1.x, P1.y);
  printf("P2 = (%ld, %ld)\n", P2.x, P2.y);
  i++;
  Ni = q;
  goto L2;
  L9:
  if (i == 0) return 0;
  i--;
  goto L3;
}

/*
void test_1(void)
{
  long a = 1, b = 6, d, p = 11;
  struct point P = {2, 7}, R;

  multiply(a, 11, p, P, &R, &d);
  printf("%ld(%ld, %ld) = (%ld, %ld)\n", 11, P.x, P.y, R.x, R.y);
  multiply(a, 12, p, P, &R, &d);
  printf("%ld(%ld, %ld) = (%ld, %ld)\n", 12, P.x, P.y, R.x, R.y);
  printf("|E(F_%ld)| = %ld\n", p, cardinality(a, b, p));
  printf("Rabin_Miller(5) = %d\n", Rabin_Miller(5));
}

void test_2(void)
{
  long a = 1, b = 2, d, p = 5;
  struct point P = {1, 3}, R;

  if (multiply(a, 2, p, P, &R, &d))
    printf("d = %ld\n", d);
  printf("%ld(%ld, %ld) = (%ld, %ld)\n", 2, P.x, P.y, R.x, R.y);
  if (multiply(a, 4, p, P, &R, &d))
    printf("d = %ld\n", d);
  printf("%ld(%ld, %ld) = (%ld, %ld)\n", 4, P.x, P.y, R.x, R.y);
  printf("|E(F_%ld)| = %ld\n", p, cardinality(a, b, p));
}
*/

int main(void)
{
  long N;

  /* test_1();
  test_2(); */
  for (;;) {
    printf("number to be tested or 0 to quit:\n");
    scanf("%ld", &N);
    if (N == 0) break;
    if (Goldwasser_Kilian(N))
      printf("proven prime\n");
    else
      printf("composite\n");
  }
  return 0;
}