/*
  Author:  Pate Williams (c) 1997

  "Exercise 6.3. Controlling the factor base.
  Utilizing the function jacobi on p. 283 write
  a computer program helping you, for any given
  N to find the most efficient (square-free)
  multiplier k <= 1000, such that a maximal
  number of primes <= 100 will belong to the
  factor base." -Hans Riesel-
  See "Prime Numbers and Computer Methods for
  Factorization" by Hans Riesel second edition
  page 202.
*/

#include <stdio.h>

#define BITS_PER_LONG   32
#define BITS_PER_LONG_1 31

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 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(0l, 0l, sieve);
  set_bit(1l, 0l, sieve);
  set_bit(2l, 1l, 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, 0l, sieve);
      i += inc;
    }
    c += 2;
    while (!get_bit(c, sieve)) c++;
  } while (c * c <= n);
}

long get_k(long N, long f_max, long k_max,
           long *f_sieve, long *k_sieve,
           long *fb, long *fb_size)
{
  long K, c_max = 0, k = 1;
  long count, f_count = 0, kN, p;

  p = 2;
  while (p <= f_max) {
    while (!get_bit(p, f_sieve)) p++;
    f_count++;
    p++;
  }
  while (k <= k_max) {
    count = 0;
    kN = k * N;
    p = 1;
    do {
      while (!get_bit(p, f_sieve)) p++;
      if (JACOBI(kN, p) != - 1)
        fb[count++] = p;
      p++;
    } while (p <= f_max);
    if (count > c_max) {
      K = k;
      c_max = count;
      if (count == f_count) break;
    }
    k++;
    if (k <= k_max)
      while (!get_bit(k, k_sieve)) k++;
  }
  count = 0;
  kN = K * N;
  p = 1;
  do {
    while (!get_bit(p, f_sieve)) p++;
    if (JACOBI(kN, p) != - 1)
      fb[count++] = p;
    p++;
  } while (p <= f_max);
  *fb_size = count;
  return K;
}

int main(void)
{
  long f_max, fb_size, i, k, k_max, N;
  long f_sieve[1024], k_sieve[1024];
  long fb[1024];

  Sieve(100, f_sieve);
  Sieve(1000, k_sieve);
  f_max = 100;
  while (!get_bit(f_max, f_sieve)) f_max--;
  k_max = 1000;
  while (!get_bit(k_max, k_sieve)) k_max--;
  for (;;) {
    printf("N = ");
    scanf("%ld", &N);
    if (N == 0) return 0;
    k = get_k(N, f_max, k_max, f_sieve, k_sieve, fb, &fb_size);
    printf("k = %ld\n", k);
    printf("the factor base is:\n");
    for (i = 0; i < fb_size; i++) {
      printf("%2ld ", fb[i]);
      if ((i + 1) % 5 == 0) printf("\n");
    }
    if (fb_size % 5 != 0) printf("\n");
    printf("total number of factors is: %ld\n", fb_size);
  }
}