Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- /* 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;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement