/*  FREUDENTHAL GENERATOR
 *
 * Author: Brett Berger
 * 28 March 2016
 *
 * Included is a very efficient method of generating solutions to the sam polly
 * problem, which I have learned is more properly known as the Freudenthal
 * problem.  The statement is reproduced below, generalized.
 *
 * ~~~~~~~ The Freudenthal problem ~~~~~~~~~
 *
 * Sam hears the sum of two numbers, Polly the product.  The numbers
 * are known to both be at least NMIN.
 *
 * S: "You don't know the numbers"
 * P: "That was true, but now I do"
 * S: "Now I do too"
 *
 * Classify possible pairs of numbers.
 *
 * In the formalism of https://www.win.tue.nl/~gwoegi/papers/freudenthal1.pdf
 * this is the problem FREUDENTHAL^=(NMIN,*), searching for stable solutions
 * with a given NMIN and without the constraint that the two numbers differ.
 *
 * ~~~~~ Overview of a naive approach: ~~~~~
 *
 * First, generate all possible products and sums of numbers greater than NMIN
 * with the product bounded above by some nmax, and order the pairs by product.
 *
 * Then, for products with only one possible sum, (i.e., only one factorization
 * given the NMIN constraint), tag the resulting sum as invalid, because Sam,
 * if given that sum, wouldn't be able to say Polly doesn't know the numbers.
 * (example: if NMIN is 3, 52 can only be written as 4 x 13.  Thus, the sum 17
 * is invalid, because there is a possible pair of numbers that give Polly the
 * answer right away.)
 *
 * Next, for products with multiple possible sums, look for those with exactly
 * one sum which has not been tagged as invalid.  (example: if NMIN is 3, then
 * 42 can be written as 3 x 14 or 6 x 7, so possible sums are 13 and 17. We
 * already showed that 17 is invalid, but 13 can be shown valid.  Thus 42 would
 * be accepted at this stage).  Save the product and corresponding valid sum,
 * now ordered by sum.
 *
 * For all sums that only appear with one viable product, we get a solution.
 * (example: if NMIN is 3, 13 could correspond to 42, but also to 30, 36, and
 * 40, so there is no solution there.  Sam can't learn the two numbers from
 * the fact that Polly figured them out.  However, 29 is a valid number that
 * corresponds to only the product 208.  Thus, 29, 208 gives us the solution
 * 13, 16, which is the smallest solution with NMIN set to 3)
 *
 * ~~~~~~ Why is that so inefficient? ~~~~~~
 *
 * Unfortunately, the upper bound nmax means the range over which we can be
 * certain solutions are stable is very small.  For products up to nmax, we
 * only exhaust the options for sums up to 2 * sqrt(nmax).  However, these
 * sums correspond to products as low as (2 * sqrt(nmax) - NMIN)  * NMIN.
 * So out of our nmax products, we only can say that they are uniquely fixed
 * by Sam's info for O(sqrt(nmax)) of them.
 *
 * Next, the sums.  We also need to be certain of uniqueness, but any sum can
 * correspond to a product on the order of its square.  The list of products
 * which are uniquely fixed is O(sqrt(nmax)), so the list of sums which we can
 * confidently say are unique is O(sqrt(sqrt(nmax))).
 *
 * This large amount of space overhead becomes a problem.  To reach confidence
 * up to a given sum, we need to compute and store O(sum^4) values.  It turns
 * out that limits us to sums around ~350 before running into system memory
 * constraints, at which point we have found the following solution pairs:
 *
 * {13, 16}
 * {16, 73}
 * {16, 133}
 * {16, 163}
 * {64, 127}
 * {16, 193}
 * {16, 223}
 * {13, 256}
 * {64, 241}
 *
 * ~~~~~~~~~~~~~ Altenatives: ~~~~~~~~~~~~~~
 *
 * I was initially drawn to the above method because it was simple, products
 * were generated from the factors instead of the factors from the products
 * (factoring gets slow).  But not allowing ourselves to factor large numbers
 * means that we have to wait until our program has reached all possible
 * factorizations of the number before we can say anything about it.
 *
 * Looking for alternatives, the next tempting idea is to start with the sums:
 * Once we know what the sum can be, we can generate a large collection of
 * corresponding products.  We remove duplicates, which has the same effect
 * as finding factorizations with exactly one matching sum (example, if
 * NMIN is 3, 13 is a possible sum, generating 30, 36, 40, and 42, none of
 * which show up elsewhere.  29 is a possible sum, but 100 -> 25, 120 -> 43,
 * 138 -> 49, 154 -> 25, 168 -> 31, 180 -> 49, 190 -> 43, 198 -> 31, 204 -> 55,
 * 210 -> 73, 37 and 31... leaving only 208.  For a given sum, that still means
 * we need to have this set known to a length of sum^2 / (4 * NMIN) + NMIN,
 * which for 29 is just 73 but for 350 is 10211.  However, this is still more
 * space efficient since we are generating the products from the sums instead
 * of a gargantuan lookup table of factorizations.
 *
 * ~~~~~ What can Sam's sum look like? ~~~~~
 *
 * This discussion almost exclusively looks at NMIN = 3, my favorite case.
 *
 * Our collection of valid sums looks basically like every odd composite number
 * plus 4: 13, 19, 25, 29, 31, 37, 43, 49, 53, 55, 59, 61, 67, 73, 79, etc..
 * except it excludes some: 39 and 69 are missing from the above.  This is
 * because, for 39, the sum 13 + 26 corresponds to a product 13 x 26 which is
 * unique, since the 2 cannot be moved onto the other factor.  In general,
 * exclude (p + 1) x q where p and q are prime and p is less than NMIN.
 *
 * However, we know for p > 2 that p + 1 is even, and so from the list of /odd/
 * composite numbers plus 4, we only actually exclude 3 x p for prime p.
 *
 * The fact that no even numbers appear is due to the expression of even
 * numbers as the sum of two primes.  However, this relies on the Goldbach
 * conjecture for NMIN < 4, and if NMIN is increased, there might be valid even
 * sums if all Goldbach partitions have one prime less than NMIN.  For NMIN = 3
 * the only assumption is the Goldbach conjecture, which has been verified
 * up to sums of about 10^18, far larger than we will run into.
 *
 * At NMIN = 3, sums can be generated simply by taking the products of all odd
 * integers, adding 4, and verifying that the result is not of the form 3 x p.
 * This is implemented via a sieve on a large bitset.
 *
 * For compactness, we store in s[k] the validity of 2k + 5.
 * Then (2a + 1) * (2b + 1) + 4, which is an odd composite plus 4,
 * is equal to 4ab + 2a + 2b + 5, and lives at index 2ab + a + b.
 * Using integer division, this is (2a + 1) * (2b + 1) / 2.  Take
 * products of odd numbers, divide by two, mark that index as true.
 *
 * Then, we can check each true number for full validity by checking its
 * remainder mod 3: 2k + 5 is divisible by 3 if k is congruent to 2,
 * and (2 * (3a + 2) + 5) / 3 = 2a + 3 = 2 * (k / 3) + 3 (integer division)
 * = 2 * (k / 3 + 1) + 5 - 4 is not composite - known from whether
 * s[k / 3 + 1] is true.
 *
 * Let's verify this works.  Row 1 is the index, row 2 its implied value,
 * row 3 a T if tagged as composite plus 4, and row 4 the index to check
 * to see whether s/3 is composite.
 *
 * 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 * 5  7  9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45
 *             T        T        T     T  T        T  T     T
 *       1        2        3        4        5        6        7
 *
 * Row 4 aligns with a tagged value at the 6th index, so we set s[17] = s[6],
 * rendering it false.  As we know, 39 is invalid. Looks good!
 *
 * But wait, we are going through and modifying data we might rely on
 * later - indeed, if we set false s[62] because 129, despite being
 * 125 + 4, is also 3 x 43, then we later reach 375, which is 7 x 53 + 4,
 * but also 125 x 3.  By that point, 125 "looks" prime to our list, since
 * we marked 129 false, and so 375 gets wrongly marked false.  However,
 * if we do the second step of the sieve in the reverse direction, every
 * odd prime can be looked up accurately.
 *
 * As a consequence, this method generates the sum list in O(S) space and
 * S log S time, where the log comes from multiplying all odd pairs that
 * have a product less than S: that's S / a operations for a up to sqrt(S)
 * which gives a harmonic series.  Also, space is fantastically optimized;
 * to reach sums of N takes at most N/2 bits.
 *
 * ~~~~~~~~~ Beyond the Sum list ~~~~~~~~~~~
 *
 * But generating sums doesn't solve the problem, just makes it easier.
 * We still need to make relevant products.  We can use a second bitset
 * to keep track of products and their uniqueness.
 *
 * For each sum, we tag all a * b with a + b = s.  We also need to tag whether
 * that product has already been seen before.  In the interest of time
 * efficiency, it would be nice to save for each product the corresponding
 * sum.  However, that turns the space requirement from 2 bits per product
 * into 32 or more.  Given that we want to monopolize system memory, we get
 * more out of it if we don't store the sum.
 *
 * Observe that since all the sums are odd, all pairs of numbers that add
 * to a given sum include one even number, so the products of interest are
 * all even.  Since we are storing 2 bits apiece, we hold information for
 * the product 2k in indices 2k and 2k + 1 as follows:
 *
 * Unseen: 00
 * Seen:   01
 * Reseen: 11
 * updated by ps[2k + 1] = ps[2k], ps[2k] = true;
 *
 * Regarding efficiency of generating products:
 * We can make a * (s - a) without multiplying every time, as follows:
 * 3 * (s - 3) is 3 * s - 9, then add s - 7, s - 9, etc until adding 0.
 * Since we compute with s = 2 * a + 5, we can translate the procedure into
 * functions of a instead of s.
 *
 * We are generating a fixed-size lookup table, so, we only insert products
 * if we know we will get full info on them, at most NMIN * (S - NMIN), where
 * S is the largest sum value we allow.
 *
 * ~~~~~ "Now I do," "Now I do too" ~~~~~~~~
 *
 * We want to use these giant lookup tables to quickly generate solutions.
 * For a given sum, we iterate over all products as before.  If the product
 * is tagged with 01, then we know:
 *
 * It has exactly one factorization for which the factors sum to an element of
 * the sum list - the element that brought us there.
 *
 * Because it is generated from the sum list, we know it has at least one other
 * factorization that does not sum to an element of the sum list.
 *
 * Thus, for this product, Polly can say "That was true, but now I do."
 *
 * There are three options when iterating over all products for a sum: either
 * we find 01 never, once, or more than once.  If we find it exactly once, then
 * Sam can say "Now I do too."  If we find it a second time, we can break early
 * from the loop.  If we reach the end with exactly one 01, we have a pair of
 * numbers satisfying the problem statement, as desired.  We save it to a set
 * for later processing of its properties.
 *
 * For further analysis, we provide a second operating mode that doesn't
 * display the solution pairs but instead displays for each sum the number of
 * times 01 is encountered.  We call this the Freudenthal partition function
 * and its behavior is detailed in comethelper.cpp.
 *
 * ~~~~~~~~~~ Once space runs out ~~~~~~~~~~
 *
 * The above has an upper limit: encoutering products not saved in our lookup
 * table.  This may happen at sums above 2 * sqrt(NMIN * (S - NMIN)).  It
 * seems like a bit of a waste then, to have such a vast sum list but only
 * scan through O(sqrt(S)) of it.
 *
 * Fortunately, there are memory-efficient time-inefficient ways to check
 * pairs once we've exceeded this limit.  This will be reasonably fast
 * at first, since it can reject any sum which finds two uniquely expressed
 * products, and it begins with the lowest products, which are still in the
 * lookup table.  But in general, once we exceed our limit, we need to do
 * explicit factoring.  We implement a slow p check and s check that will work
 * for any number.  We scan factors of p.  Every sum: if it is within our sum
 * lookup table, we check the table, otherwise, we factor s - 4 and s / 3 if
 * applicable.  If we find exactly one good sum, we return as though p was
 * tagged 01.  It just potentially takes O(p) time instead of constant.
 *
 * We use this p check to fill out the rest of the sum list.  But since we
 * make our S list humongous and we don't want to wait for completion, we
 * allow SIGINT to kill this function without killing the program.  After
 * SIGINT, the program dumps the sorted list of all solutions to cout, clears
 * memory, and cleanly exits.
 *
 * ~~~~~~~~ Complexity analysis ~~~~~~~~~~~~
 *
 * At the present version, we have a worst-case complexity analysis below:
 *
 * space: O(S), allocated in advance (no mid-program segfaults please)
 *
 * Generating our s list:
 *             sieve1: every odd N < S contributes S / N            -> S log S
 *             sieve2: every 3rd element contributes constant       -> S
 * Generating products:
 *             genprod: every N < S contributes min(N, S / N)       -> S log S
 *             gensols: every N < sqrt(S) contributes at most N     -> S
 *             gensolsslow: N > sqrt(S) may contribute S / N constant
 *                          and (N - sqrt(S))^2 calls to slowPcheck
 *                          with argument between S and N^2
 *                  slowPcheck: may take arg / S calls to slowScheck
 *                              with argument between S and arg
 *                  slowScheck: may take sqrt(arg) time
 *
 *  Thus, we generate solutions with sum up to order sqrt(S) in S log S time,
 *   the next N take roughly N^3 sqrt(S) time (valid only for N << sqrt(S)),
 *   and to generate all solutions with sum up to O(S), the form becomes S^5.
 *
 * Overnight, with space as large as my laptop could handle, this generated
 * 3141 solutions in FREUDENTHAL^=(3,*), with sum up to roughly 2 ^ 19.
 *
 * ~~~~~~~~~~~~~~~ NMIN > 3 ~~~~~~~~~~~~~~~~
 * It is not yet known whether solutions exist with NMIN > 3.  The paper I
 * referenced above found none with sum less than 50000 and conjectured none
 * exist.  We can follow the same procedure in general, but the S list must be
 * modified.
 *
 * Call a number pseudoprime in our case if it has no factorizations with at
 * least one factor greater than or equal to NMIN.  For NMIN = 3, that makes
 * 4 the only nonprime pseudoprime.  For NMIN = 4, the list is 4, 6, 9.  For
 * NMIN = 5, we get 6, 8, 9.  The largest even pseudoprime for a given NMIN can
 * be shown to be 2 * (NMIN - 1).  For a pseudoprime p and a prime q, p * q has
 * only one allowable factorization and so p + q is excluded from S.  The sum
 * of two pseudoprimes is a weird case.  It is disallowed if there is no way
 * to legally exchange factors. For NMIN = 6 we have 6 and 8 as pseudoprimes,
 * and 48 factors uniquely, but 15 is also a pseudoprime and 8 * 15 = 120 is
 * expressible as 6 * 20.  I don't see a particularly systematic way to find
 * this, so we can just check every option.
 *
 * We only need to look at even primes and pseudoprimes in the NMIN = 2 and 3
 * case because we were safe to assume that every even number at least 6 is the
 * sum of two odd primes.  But it could be the case that a given number is not
 * the sum of two odd primes both of which are at least NMIN.  For example,
 * if NMIN = 6, then 16 is out because both its partitions rely on smaller
 * primes.
 *
 * A small list of largeish even numbers made valid by raising NMIN:
 *     NMIN = 14, 44 becomes valid.  With NMIN = 14, pseudoprimes include 169.
 *     NMIN = 32, 92 becomes valid.  (pseudoprimes up to 961)
 *     NMIN = 104, 242 becomes valid.
 *
 * From muddling about a bit with goldbach partitions, it seems like a very
 * safe assumption that if we check all possible sums of pseudoprimes then we
 * can assume any even number beyond that range is the sum of two valid primes.
 *
 * ~~~~~ Respecialization to NMIN = 4 ~~~~~~
 *
 * With NMIN = 4, we have pseudoprimes 4, 6, and 9.  By hand we check that
 * 16 = 4 x 4 eliminates 8, 25 = 5 x 5 eliminates 10, 35 eliminates 12,
 * 49 eliminates 14, 55 eliminates 16, and 65 eliminates 18.  Our sum list
 * consists of odd numbers for which s - 4 and s - 6 are both composite
 * and which are not thrice a prime.  This is just one further condition
 * than we encountered for NMIN = 3, and can be implemented by altering
 * the way sieve2 works to also window every number with its neighbor.
 *
 * If we want to be able to use our sieve trick to eliminate numbers of the
 * form 3 * p, then we need to work backwards on our S list.  I'm debating
 * whether we want to have s[k] hold 2k + 5 or 2 * (k + NMIN) - 1 in general.
 * If NMIN = 2 or 3 the latter is ideal.  If NMIN = 4 or 5, then after the
 * sieve1 we need to have s[k]_post = (s[k] & s[k + 1])_pre.  This is difficult
 * to accomplish if we work backwards, but I guess we just store a value before
 * overwriting.
 *
 * The benefit of keeping 2k + 5 in the specialized case of NMIN = 4 is that
 * we window in the appropriate direction for a reverse scan.  All unscanned
 * indices hold relevant information about primality.  But generalizing with
 * this goal would mean that for an arbitrary NMIN we would have s[k] storing
 * 2k + ((NMIN + 1) & ~1) + 1.  That makes prime lookups awkward, but
 * if we specialize above NMIN = 5 it means we don't need to think about
 * storing arbitrary overhead information for our sieve.
 *
 * Let's try this latter goal, and look at a few sum lists in the making:
 *
 * NMIN = 2: tag composites plus 2 at index (c / 2).
 * 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 * 3  5  7  9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43
 *             T        T        T     T  T        T  T     T
 *
 * NMIN = 3: tag composites plus 4 at index (c / 2)
 * 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 * 5  7  9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45
 *             T        T        T     T  T        T  T     T
 * Then, as above, for any index congruent to 2 mod 3, look at (k + 1) / 3 for
 * information about the primality of s / 3.
 *
 * NMIN = 4: tag composites plus 4 at index (c / 2)
 * 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 * 5  7  9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45
 *             T        T        T     T  T        T  T     T
 * Then, look at k - 1 and, for any index congruent to 2 mod 3, (k + 1) / 3.
 * Perform bitwise and with both.
 *
 * NMIN = 5: tag composites plus 6 at index (c / 2)
 * 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 * 7  9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47
 *             T        T        T     T  T        T  T     T
 * Look at k - 1 and, for any index congruent to 1 mod 3, (k + 2) / 3.
 * Perform bitwise and with both.
 *
 * NMIN = 6: tag composites plus 6 at index (c / 2)
 * 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 * 7  9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47
 *             T        T        T     T  T        T  T     T
 * Look at k - 1 and k - 2, and for index congruent to 1 mod 3, (k + 2) / 3.
 * Perform bitwise and with both.
 *
 * NMIN = 7: tag composites plus 8 at index (c / 2)
 * 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 * 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
 *             T        T        T     T  T        T  T     T
 * Look at k - 1 and k - 2, and for index congruent to 0 mod 3, (k + 3) / 3.
 * etc..
 *
 * General: tag composites plus (NMIN + 1) & ~1 at index (c / 2)
 * Look at (distances to other even pseudoprimes, all backward), and for index
 * with k + (NMIN - 1) / 2 divisible by three, look at k / 3 + (NMIN - 1) / 6.
 *
 * Then when creating products, we need to begin with NMIN * (s - NMIN), which
 * is NMIN * (2 * k + 1 + ((NMIN + 1) & ~1) - NMIN), which is given most simply
 * by NMIN * (2 * k + (NMIN % 2) + 1).  Example, with NMIN = 7, k = 14, we
 * have 7 * (2 * 14 + 2) = 210, which is the smallest number with a factor sum
 * of 37, which is what s[14] holds for NMIN = 7.
 *
 * Incrementing from there we need (NMIN + 1) * (s - NMIN - 1), which is the
 * prior value plus s - 2 * NMIN - 1.  Again using our expression for s from k,
 * we get 2 * k - NMIN - ((NMIN % 2) ? 1 : 0).  These ternaries are pretty ugly
 * but fortunately the precompiler is the one evaluating them so our code will
 * still perform well.  Actually since this expression is always even, we can
 * use integer division to write 2 * (k - NMIN / 2).  This value on the right
 * decreases by one each iteration to zero to list off the products.  Example:
 * 2 * (14 - 3) = 22 adds to 210 to give us 232, which is 8 x 29.
 *
 * We ignore for the moment the potential for even numbers and sums of
 * pseudoprimes to be allowed sums and just focus on implementing this general
 * algorithm.  Presently, the code functions by assuming the only concerning
 * pseudoprimes are even numbers between NMIN and 2*(NMIN - 1).  I have yet to
 * find a case where this assumption affects the outcome.
 */

#include <iostream>
#include <cstdint> // For definite bit widths
#include <set>
#include <bitset>
#include <cmath>
#include <csignal>

// We have found S can be around 16 billion / NMIN before memory restrictions
#define NMIN 3
#define S 409ull

// The VERBOSE option allows progress updates to be sent to std::cerr.
#define VERBOSE true

/* We provide two modes.  If COMET is set to false, then std::cout recieves
 * solution pairs.  If COMET is set to true, then std::cout instead receives
 * ordered pairs of the form s, f(s) where f(s) is the Freudenthal partition
 * function (analogous to the Goldbach partition function), counting the
 * number of products which are uniquely fixed by Sam's claim.
 */
#define COMET true

/* Function nums : generates the original numbers given their sum and product
 * arguments: sum and product
 * return   : packaged form of number pairs.  If a + b = s, ab = p, a is even,
 *            then bits 0 - 31 hold b and bits 32 - 63 hold a.  This allows us
 *            to sort the outputs by their even number, breaking ties by the
 *            odd number.
 */
uint64_t nums(uint64_t s, uint64_t p)
{
    int64_t b(sqrt(s * s - 4 * p));
    if ((s + b) % 4)
        b = -b;
    return ((s + b) << 31) + ((s - b) >> 1);
}

/* Section headers of the top comment are reproduced below to give reference
 * for a more detailed discussion of functions in the section.
 */

// ~~~~~ What can Sam's sum look like? ~~~~~

/* Function sieve1 : implements an erastosthenes sieve on odd numbers.
 * argument: pointer to bitset with all bits off
 * post    : bits for which 2 * a + 1 is composite are set
 */
void sieve1(std::bitset<S/2> * s)
{
    uint64_t b;
    for (uint64_t a = 3;; a += 2)
    {
        b = a * a / 2;
        if (b > s->size())
            break;
        for (; b <= s->size(); b += a)
            (*s)[b] = true;
    }
}

/* Function window : A helper function for sieve2 with NMIN > 3.
 * note : I'm really hoping -O3 optimizes this whole function out and just
 * replaces it with an unrolled loop with short-circuited booleans.
 * Particularly for NMIN = 3, the loop doesn't acually run so the relevant
 * machine code is far simpler if optimized.
 */

 inline void window(std::bitset<S/2> * s, uint64_t index)
 {
    for (int a = 1; a < NMIN / 2; ++a)
        (*s)[index] = (*s)[index] && (*s)[index - a];
 }

/* Function sieve2 : Performs a reverse scan to generate sum lookup table
 * argument: pointer to bitset with composite numbers set
 * post    : bitset contains all possible sums for which no pair of summands
 *           may be uniquely determined from their product.
 */
void sieve2(std::bitset<S/2> * s)
{
    // We have two goals: window every index and take care of the 3p case.
    // So that our loop may have definite structure, we enforce that b always
    // holds an index encoding a multiple of 3 and a holds the index holding
    // the primality of the encoded sum divided by 3.

    // The highest b in our array is one of the top 3 indices, satisfying
    // b + (NMIN - 1) / 2 == 3 * a.
    uint64_t b = S / 2 - 1;
    if ((b + (NMIN - 1) / 2) % 3 > 1)
        window(s, b--);
    if ((b + (NMIN - 1) / 2) % 3 > 0)
        window(s, b--);

    uint64_t a = (b + (NMIN - 1) / 2) / 3;
    while (b > 2)
    {
        (*s)[b] = (*s)[b] && (*s)[a];
        window(s, b--);
        window(s, b--);
        window(s, b--);
        --a;
    }

/* Displays the sum list to cout: a nice sanity check when debugging.
 */
}
void writeS(std::bitset<S/2> * s)
{
    for (int i = 0; i < S/2; ++i)
        if ((*s)[i])
            std::cout << (2 * i + ((NMIN + 1) | 1)) << ", ";
    std::cout << std::endl;
}

// ~~~~~~~~~ Beyond the Sum list ~~~~~~~~~~~

/* Function genprod : tags products according to how many sums generate them
 * arguments: pointer to filled sum bitset and empty product bitset
 * post     : product bitset is tagged appropriately
 */
void genprod(std::bitset<S/2> * s, std::bitset<NMIN*S> * ps)
{
    uint64_t c, d;
    uint64_t cmax = NMIN * (S - NMIN);
    for (uint64_t a = 0; a < S / 2;)
        if ((*s)[a++])
            for (d = a - NMIN / 2, c = NMIN * (2 * a - 1 + NMIN % 2);
                 d && c <= cmax; c += (--d << 1))
            {
                (*ps)[c + 1] = (*ps)[c];
                (*ps)[c] = true;
            }
}

// ~~~~~ "Now I do," "Now I do too" ~~~~~~~~

/* Function gensols : scans sum and product tables for solutions
 * arguments: pointers to sum and product tables, reference to solution set
 * post     : solution set is filled to sums at most 2 * sqrt(NMIN * (S - NMIN))
 */
void gensols(std::bitset<S/2> * s, std::bitset<NMIN*S> * ps,
             std::set<uint64_t> &solutions)
{
    uint64_t c, d;
    const uint64_t smax = sqrt(NMIN * (S - NMIN));
    uint64_t uniqc;
    for (uint64_t a = 0; a < smax - NMIN;)
        if ((*s)[a++])
        {
            uniqc = 0;  // The solution product, or f(s) in COMET mode.
            for (d = a - NMIN / 2, c = NMIN * (2 * a - 1 + NMIN % 2);
                 d; c += (--d << 1))
                if (!(*ps)[c + 1])
                {
                    if (COMET)
                        ++uniqc;
                    else
                    {
                        if (uniqc)
                            break;
                        else
                            uniqc = c;
                    }
                }
            if (COMET)
                std::cout << 2 * a + ((NMIN - 1) | 1) << ", "
                          << uniqc << std::endl;
            else if (!d && uniqc)
                solutions.insert(nums(2 * a + ((NMIN - 1) | 1), uniqc));
        }
}

//~~~~~~~~~ Once space runs out ~~~~~~~~~~~

/* Function isOddPrime : helper to slowScheck
 * argument: an odd number
 * return  : true if prime, false otherwise
 */
bool isOddPrime(uint64_t oddo)
{
    uint64_t amax = sqrt(oddo);
    for (uint64_t a = 3; a <= amax; ++a)
        if (!(oddo % a))
            return false;
    return true;
}

/* Function slowScheck : determines if a number is a valid sum
 * argument: an odd number
 * return  : true if valid (emulates (*s)[sum / 2 - 2])
 */
bool slowScheck(uint64_t sum)
{
    if (NMIN == 2)
        return !isOddPrime(sum - 2);
    if (sum % 3 == 0 && isOddPrime(sum / 3))
        return false;
    for (int a = 2 * (NMIN - 1); a >= NMIN; a -= 2)
        if (isOddPrime(sum - a))
            return false;
    return true;
}

/* Function slowPcheck: Does p factor into any valid sum besides the one given?
 * arguments: a product and sum pair, and pointer to the known valid sums.
 * return   : true if another pair is found, false if p is uniquely determined
 *            (emulates (*ps)[p + 1])
 * note     : first considers factors with minimum sum to minimize chance of
 *            needing to call slowScheck.
 */
bool slowPcheck(uint64_t p, uint64_t s1, std::bitset<S/2> * s)
{
    uint64_t sum;
    for (uint64_t a = sqrt(p); a >= NMIN; --a)
        // To proceed, p must divide a, the factor sum must be odd and not s1,
        // and if all those succeed, we can check whether the sum is a valid
        // s number.  If it is, then we know p isn't uniquely determined by
        // the provided info.
        if (!(p % a) && ((sum = a + p / a) % 2) && sum != s1
            && (sum < S && (*s)[(sum - NMIN - 1) / 2]
                              || sum >= S && slowScheck(sum)))
            return true;
    return false;
}

/* Function gensolsslow : continues where gensols left off but slower
 * arguments: same as gensols
 * note     : as it runs long, we permit it to be interrupted by SIGINT
 * return   : the largest sum that was checked to completion
 *
 * aside    : In solution mode, it quits after finding a single point of
 *            inconsistency with Sam's 3rd claim, which happens very quickly
 *            for NMIN > 3.  Thus, we are able to verify the absence of
 *            solutions up to a very significant bound on the sum.  In
 *            comet mode, we must check every product, so this is far slower.
 */
volatile sig_atomic_t flag = 1;
void exitSlowFunction(int sig)
{
    flag = 0;
}
uint64_t gensolsslow(std::bitset<S/2> * s, std::bitset<NMIN*S> * ps,
                     std::set<uint64_t> &solutions)
{
    uint64_t c, d;
    uint64_t cmax = NMIN * (S - NMIN);
    uint64_t a = sqrt(cmax) - NMIN;
    uint64_t uniqc;
    while (flag && a < S/2)
        if ((*s)[a++])
        {
            uniqc = 0;
            for (d = a - NMIN / 2, c = NMIN * (2 * a - 1 + NMIN % 2);
                 d; c += (--d << 1))
                if (c <= cmax && !(*ps)[c + 1]
            || (c > cmax && !slowPcheck(c, 2 * a + ((NMIN - 1) | 1), s)))
                {
                    if (COMET)
                        ++uniqc;
                    else
                    {
                        if (uniqc)
                            break;
                        else
                            uniqc = c;
                    }
                }
            if (COMET)
                std::cout << 2 * a + ((NMIN - 1) | 1) << ", "
                          << uniqc << std::endl;
            else if (!d && uniqc)
            {
                solutions.insert(nums(2 * a + ((NMIN - 1) | 1), uniqc));
                if (NMIN > 3)
                    std::cerr << "Counterexample found!!" << std::endl;
            }
        }
    return 2 * a + ((NMIN - 1) | 1);
}

/* Procedure that utilizes all the above.
 */
int main()
{
    signal(SIGINT, exitSlowFunction);

    std::bitset<S/2> * s = new std::bitset<S/2>;
    std::bitset<NMIN*S> * ps = new std::bitset<NMIN*S>;

    if (VERBOSE)
        std::cerr << "prepared memory successfully" << std::endl;

    sieve1(s);
    if (NMIN > 2)
        sieve2(s);

    if (VERBOSE)
        std::cerr << "sieving complete, S list generated" << std::endl;

    genprod(s, ps);

    if (VERBOSE)
        std::cerr << "tagged products up to " << ps->size() << std::endl;

    std::set<uint64_t> solutions;
    gensols(s, ps, solutions);

    if (VERBOSE)
        std::cerr << solutions.size()
                  << " solutions with sum at most "
                  << (uint64_t) (2 * sqrt(NMIN * (S - NMIN))) << std::endl;

    if (VERBOSE && !COMET)
        std::cerr << "Product table exhausted, further solutions will be "
                  << "generated far more slowly" << std::endl;

    uint64_t maxsum = gensolsslow(s, ps, solutions);

    if (VERBOSE && !COMET)
        std::cerr << ((flag) ? "Maxed out sum table" : "SIGINT received")
                  << " at a sum of " << maxsum
                  << ", writing " << solutions.size()
                  << " solutions" << std::endl;

    if (VERBOSE && COMET)
        std::cerr << "Reached sum of " << maxsum;


    // Clear memory and present solution pairs.
    delete s;
    delete ps;
    if (!COMET)
        for (const auto& i : solutions)
            std::cout << (i >> 32) << ", " << (i & 0xffffffff) << std::endl;

    return 0;
}