// This program computes terms of A098550.  A later comment describes
// the rough approach.

// The program takes two optional arguments.  The first is when to
// stop; by default it runs forever.  The second is when to start
// printing; by default it starts at start.  The program skips the
// first three terms (1,2,3), so the output starts at n=4.

// So, for example, ("$ " is the prompt in these examples, to distinguish
// commands from output, and ./A098550 is the name of the compiled executable)

// $ ./A098550 10
// 4 4
// 5 9
// 6 8
// 7 15
// 8 14
// 9 5
// 10 6

// $ ./A098550 1000000 999995
// 999995 960002
// 999996 1094335
// 999997 959994
// 999998 1093931
// 999999 960000
// 1000000 1094537

// There are some additional outputs you can enable by editing the code
// and recompiling, controlled by some boolean flags at the start of the
// program.  For example, if you change print_factorizations and
// print_ratios to true, you get

// $ ./A098550 1000000 999995                             
// 999995 960002 2*37*12973 0.960007
// 999996 1094335 5*11*101*197 1.09434
// 999997 959994 2*3*3*7*19*401 0.959997
// 999998 1093931 101*10831 1.09393
// 999999 960000 2*2*2*2*2*2*2*2*2*3*5*5*5*5 0.960001
// 1000000 1094537 101*10837 1.09454


// variables controlling optional outputs

// print the first three terms (1,2,3)
const bool print_first3 = false;

// after a(n), print the factorization of a(n)
const bool print_factorizations = false;

// after that, print a(n)/n
const bool print_ratios = false;

// after that, print |i<=n : a(i)<=i|/n
const bool print_smallcnt = false;

// after that, print the sizes of the even border, odd border,
// ready primes, and large odd prime multiples
const bool print_program_stats = false;

// check the even/odd 2p/odd/p/even/k*p pattern and report all violations
// to stderr
const bool test_pattern = false;

#include <iostream>
#include <cstdlib>
#include <list>
#include <vector>
#include <set>
#include <map>
#include <iterator>

typedef std::list<size_t>::iterator LIter;
typedef std::set<size_t>::iterator SIter;

static inline
bool isprime(size_t x) {
  if (x % 2 == 0 || x < 2) return false;
  for (size_t d = 3; d*d <= x; d += 2) {
    if (x % d == 0) return false;
  }
  return true;
}

static inline
size_t nextprime(size_t x) {
  if (x % 2 == 0) x++;
  else x += 2;
  while (!isprime(x)) x += 2;
  return x;
}

static inline
size_t gcd (size_t x, size_t y) {
  if (x < y) return gcd(y,x);
  if (y == 0) return x;
  if (y == 1) return 1;
  return gcd (y, x%y);
}

void ifactor (size_t x, std::vector<size_t>& factors)
{
  factors.clear();
  while (x % 2 == 0) {
    factors.push_back (2);
    x /= 2;
  }
  for (size_t p=3; p*p <= x; p += 2) {
    while (x % p == 0) {
      factors.push_back (p);
      x /= p;
    }
  }
  if (x > 1) factors.push_back (x);
}

// APPROACH
// This program uses the fact that primes appear (as divisors) in increasing
// order (otherwise, the first occurrence of a prime could be replaced by
// a smaller prime that has not occurred), and that a prime cannot appear
// as itself until it has appeared as a divisor of another number.

// Beyond that, for efficiency it uses the observation that
// entries in the sequence take several forms:
// even numbers, a(n) ~ 0.96n
// odd primes, a(n) ~ 0.48n (a(n-2) = 2*a(n))
// "big" odd composites, a(n) = p*q, p prime, q small odd prime,
//    a(n-2)=p, a(n) ~ 0.48qn
// normal odd composites, a(n) ~ 1.09n
//
// so we track
// a2: a(n-2)
// a1: a(n-1)
// evenskip: the frontier of even numbers - the even numbers which
//    have not appeared up through the largest even number which has appeared
// oddskip: the frontier of "normal" odd numbers (does not include odd primes
//    or "big" odd composites
// bigodd: the "big" odd composites which have appeared and have not yet been
//    absorbed by the odd frontier
// next_prime: the smallest prime that hasn't occurred as a divisor
// primeskip: the primes that have occurred as divisors but not as themselves.
// a1prime,a2prime: whether or not a(n-1)/a(n-2) was an odd prime
//    (so a(n) is a "big" odd composite if it isn't 2*a(n-2))

// globals for list state
static std::list<size_t> evenskip;
static size_t evencnt;

static std::list<size_t> oddskip;
static size_t oddcnt;

static std::set<size_t> bigodd;

static std::list<size_t> primeskip;
static size_t primecnt;
static size_t next_prime;
static size_t next_oddprime;

static inline
void preincr_evenskip (LIter i)
{
  if (*i == evenskip.back()) {
    evenskip.push_back ((*i) + 2);
    evencnt++;
  }
}

static inline
void use_evenskip (LIter i)
{
  preincr_evenskip (i);
  evenskip.erase(i);
  evencnt--;
}

static inline
void use_evenskip_grow (size_t x)
{
  for (size_t y = evenskip.back()+2; y < x; y += 2) {
    evenskip.push_back (y);
    evencnt++;
  }
  evenskip.push_back (x+2);
  evencnt++;
}

static inline
void grow_oddskip ()
{
  for (size_t x = oddskip.back()+2; ; x += 2) {
    if (x == next_oddprime) {
      next_oddprime = nextprime(next_oddprime);
    } else {
      SIter i = bigodd.find(x);
      if (i != bigodd.end()) {
	bigodd.erase(i);
      } else {
	oddskip.push_back (x);
	oddcnt++;
	return;
      }
    }
  }
}

static inline
void preincr_oddskip (LIter i)
{
  if (*i == oddskip.back()) {
    grow_oddskip();
  }
}
static inline
void use_oddskip (LIter i)
{
  preincr_oddskip(i);
  oddskip.erase(i);
  oddcnt--;
}

int main (int ac, char **av)
{
  size_t nlim = (ac > 1) ? atol(av[1]) : 0;
  size_t nout = (ac > 2) ? atol(av[2]) : 0;

  size_t a2 = 2;
  bool a2prime = false;

  size_t a1 = 3;
  bool a1prime = true;

  size_t n = 3;
  size_t smallcnt = 0;

  if (print_first3) std::cout << "1 1\n2 2\n3 3\n";
  
  evenskip.push_back(4);
  evencnt = 1;

  oddskip.push_back (9);
  oddcnt = 1;
  next_oddprime = 11;

  primecnt = 0;
  next_prime = 5;

  while (nlim == 0 || n < nlim) {
    n++;
    // candidates for a(n):
    // if a2prime, then 2*a2, or odd multiples of a2
    // else if a2 is even
    // entries in evenskip
    // else if a2 is odd
    // entries in oddskip (not in bigodd)
    // entries in primeskip
    //
    // special boundary case, only relevant early on:
    //   if the odd frontier (oddskip) passes next_prime, we put the prime
    //   into primeskip and increment, even though we haven't seen the prime
    //   as a divisor yet

    size_t found = 0;
    bool found_prime = false;

    if (a2prime) {
      for (size_t m = 2; found == 0; m++) {
	if (gcd(m, a1) != 1) continue;
	size_t x = m * a2;
	if (m % 2 == 0) {
	  LIter i = std::find(evenskip.begin(), evenskip.end(), x);
	  if (i != evenskip.end()) {
	    found = x;
	    use_evenskip (i);
	    break;
	  } else if (x > evenskip.back()) {
	    found = x;
	    use_evenskip_grow (x);
	    break;
	  }
	} else {
	  LIter i = std::find(oddskip.begin(), oddskip.end(), x);
	  if (i != oddskip.end()) {
	    found = x;
	    use_oddskip (i);
	  } else if (x > oddskip.back()) {
	    SIter j = bigodd.find(x);
	    if (j == bigodd.end()) {
	      found = x;
	      bigodd.insert (x);
	      break;
	    }
	  }
	}
      }
    } else {
      LIter i_even = evenskip.begin();
      LIter i_odd = oddskip.begin();
      LIter i_prime = primeskip.begin();
      bool use_even = (a1%2 != 0);
      while (found == 0) {
	if (i_prime != primeskip.end() && (!use_even || *i_prime < *i_even) && *i_prime < *i_odd) {
	  size_t x = *i_prime;
	  /* if (gcd(x,a2) != 1 && gcd (x,a1) == 1) */
	  if (a2 % x == 0) {
	    found = x;
	    found_prime = true;
	    primeskip.erase(i_prime);
	    primecnt--;
	    break;
	  } else {
	    ++i_prime;
	  }
	} else if (use_even && *i_even < *i_odd) {
	  size_t x = *i_even;
	  if (gcd(x,a2) != 1 && gcd(x,a1) == 1) {
	    found = x;
	    use_evenskip(i_even);
	    break;
	  } else {
	    preincr_evenskip(i_even);
	    ++i_even;
	  }
	} else {
	  size_t x = *i_odd;
	  if (gcd(x,a2) != 1 && gcd(x,a1) == 1) {
	    found = x;
	    use_oddskip(i_odd);
	    break;
	  } else {
	    preincr_oddskip(i_odd);
	    ++i_odd;
	  }
	}
      }
    }

    if (found <= n) smallcnt++;

    if (found % next_prime == 0) {
      primeskip.push_back (next_prime);
      primecnt++;
      next_prime = nextprime(next_prime);
    }
    
    if (n >= nout) {
      std::cout << n << " " << found;
      if (print_factorizations) {
	std::vector<size_t> factors;
	ifactor (found, factors);
	std::cout << " ";
	std::copy (factors.begin(), factors.end()-1, std::ostream_iterator<size_t>(std::cout, "*"));
	std::cout << factors.back();
      }
      if (print_ratios) std::cout << " " << ((double) found)/n;
      if (print_smallcnt) std::cout << " " << ((double) smallcnt)/n;
      if (print_program_stats) {
	std::cout << " " << evencnt << " " << oddcnt << " " << primecnt << " "
		  << bigodd.size();
      }
      std::cout << "\n";
       
    }

    if (test_pattern) {
      if (found_prime && a2 != found*2) {
	std::cerr << n << " " << found << " prime, a(n-2) " << a2 << " not = 2*prime\n";
      }
      if (!found_prime && found % 2 != 0 && a1 % 2 != 0) {
	std::cerr << n << " " << found << " consecutive odd numbers, second not prime\n";
      }
    }

    a2 = a1;
    a2prime = a1prime;
    a1 = found;
    a1prime = found_prime;
  }

  return 0;
}