Advertisement
Guest User

Untitled

a guest
May 12th, 2018
633
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 32.10 KB | None | 0 0
  1. /* FREUDENTHAL GENERATOR
  2. *
  3. * Author: Brett Berger
  4. * 28 March 2016
  5. *
  6. * Included is a very efficient method of generating solutions to the sam polly
  7. * problem, which I have learned is more properly known as the Freudenthal
  8. * problem. The statement is reproduced below, generalized.
  9. *
  10. * ~~~~~~~ The Freudenthal problem ~~~~~~~~~
  11. *
  12. * Sam hears the sum of two numbers, Polly the product. The numbers
  13. * are known to both be at least NMIN.
  14. *
  15. * S: "You don't know the numbers"
  16. * P: "That was true, but now I do"
  17. * S: "Now I do too"
  18. *
  19. * Classify possible pairs of numbers.
  20. *
  21. * In the formalism of https://www.win.tue.nl/~gwoegi/papers/freudenthal1.pdf
  22. * this is the problem FREUDENTHAL^=(NMIN,*), searching for stable solutions
  23. * with a given NMIN and without the constraint that the two numbers differ.
  24. *
  25. * ~~~~~ Overview of a naive approach: ~~~~~
  26. *
  27. * First, generate all possible products and sums of numbers greater than NMIN
  28. * with the product bounded above by some nmax, and order the pairs by product.
  29. *
  30. * Then, for products with only one possible sum, (i.e., only one factorization
  31. * given the NMIN constraint), tag the resulting sum as invalid, because Sam,
  32. * if given that sum, wouldn't be able to say Polly doesn't know the numbers.
  33. * (example: if NMIN is 3, 52 can only be written as 4 x 13. Thus, the sum 17
  34. * is invalid, because there is a possible pair of numbers that give Polly the
  35. * answer right away.)
  36. *
  37. * Next, for products with multiple possible sums, look for those with exactly
  38. * one sum which has not been tagged as invalid. (example: if NMIN is 3, then
  39. * 42 can be written as 3 x 14 or 6 x 7, so possible sums are 13 and 17. We
  40. * already showed that 17 is invalid, but 13 can be shown valid. Thus 42 would
  41. * be accepted at this stage). Save the product and corresponding valid sum,
  42. * now ordered by sum.
  43. *
  44. * For all sums that only appear with one viable product, we get a solution.
  45. * (example: if NMIN is 3, 13 could correspond to 42, but also to 30, 36, and
  46. * 40, so there is no solution there. Sam can't learn the two numbers from
  47. * the fact that Polly figured them out. However, 29 is a valid number that
  48. * corresponds to only the product 208. Thus, 29, 208 gives us the solution
  49. * 13, 16, which is the smallest solution with NMIN set to 3)
  50. *
  51. * ~~~~~~ Why is that so inefficient? ~~~~~~
  52. *
  53. * Unfortunately, the upper bound nmax means the range over which we can be
  54. * certain solutions are stable is very small. For products up to nmax, we
  55. * only exhaust the options for sums up to 2 * sqrt(nmax). However, these
  56. * sums correspond to products as low as (2 * sqrt(nmax) - NMIN) * NMIN.
  57. * So out of our nmax products, we only can say that they are uniquely fixed
  58. * by Sam's info for O(sqrt(nmax)) of them.
  59. *
  60. * Next, the sums. We also need to be certain of uniqueness, but any sum can
  61. * correspond to a product on the order of its square. The list of products
  62. * which are uniquely fixed is O(sqrt(nmax)), so the list of sums which we can
  63. * confidently say are unique is O(sqrt(sqrt(nmax))).
  64. *
  65. * This large amount of space overhead becomes a problem. To reach confidence
  66. * up to a given sum, we need to compute and store O(sum^4) values. It turns
  67. * out that limits us to sums around ~350 before running into system memory
  68. * constraints, at which point we have found the following solution pairs:
  69. *
  70. * {13, 16}
  71. * {16, 73}
  72. * {16, 133}
  73. * {16, 163}
  74. * {64, 127}
  75. * {16, 193}
  76. * {16, 223}
  77. * {13, 256}
  78. * {64, 241}
  79. *
  80. * ~~~~~~~~~~~~~ Altenatives: ~~~~~~~~~~~~~~
  81. *
  82. * I was initially drawn to the above method because it was simple, products
  83. * were generated from the factors instead of the factors from the products
  84. * (factoring gets slow). But not allowing ourselves to factor large numbers
  85. * means that we have to wait until our program has reached all possible
  86. * factorizations of the number before we can say anything about it.
  87. *
  88. * Looking for alternatives, the next tempting idea is to start with the sums:
  89. * Once we know what the sum can be, we can generate a large collection of
  90. * corresponding products. We remove duplicates, which has the same effect
  91. * as finding factorizations with exactly one matching sum (example, if
  92. * NMIN is 3, 13 is a possible sum, generating 30, 36, 40, and 42, none of
  93. * which show up elsewhere. 29 is a possible sum, but 100 -> 25, 120 -> 43,
  94. * 138 -> 49, 154 -> 25, 168 -> 31, 180 -> 49, 190 -> 43, 198 -> 31, 204 -> 55,
  95. * 210 -> 73, 37 and 31... leaving only 208. For a given sum, that still means
  96. * we need to have this set known to a length of sum^2 / (4 * NMIN) + NMIN,
  97. * which for 29 is just 73 but for 350 is 10211. However, this is still more
  98. * space efficient since we are generating the products from the sums instead
  99. * of a gargantuan lookup table of factorizations.
  100. *
  101. * ~~~~~ What can Sam's sum look like? ~~~~~
  102. *
  103. * This discussion almost exclusively looks at NMIN = 3, my favorite case.
  104. *
  105. * Our collection of valid sums looks basically like every odd composite number
  106. * plus 4: 13, 19, 25, 29, 31, 37, 43, 49, 53, 55, 59, 61, 67, 73, 79, etc..
  107. * except it excludes some: 39 and 69 are missing from the above. This is
  108. * because, for 39, the sum 13 + 26 corresponds to a product 13 x 26 which is
  109. * unique, since the 2 cannot be moved onto the other factor. In general,
  110. * exclude (p + 1) x q where p and q are prime and p is less than NMIN.
  111. *
  112. * However, we know for p > 2 that p + 1 is even, and so from the list of /odd/
  113. * composite numbers plus 4, we only actually exclude 3 x p for prime p.
  114. *
  115. * The fact that no even numbers appear is due to the expression of even
  116. * numbers as the sum of two primes. However, this relies on the Goldbach
  117. * conjecture for NMIN < 4, and if NMIN is increased, there might be valid even
  118. * sums if all Goldbach partitions have one prime less than NMIN. For NMIN = 3
  119. * the only assumption is the Goldbach conjecture, which has been verified
  120. * up to sums of about 10^18, far larger than we will run into.
  121. *
  122. * At NMIN = 3, sums can be generated simply by taking the products of all odd
  123. * integers, adding 4, and verifying that the result is not of the form 3 x p.
  124. * This is implemented via a sieve on a large bitset.
  125. *
  126. * For compactness, we store in s[k] the validity of 2k + 5.
  127. * Then (2a + 1) * (2b + 1) + 4, which is an odd composite plus 4,
  128. * is equal to 4ab + 2a + 2b + 5, and lives at index 2ab + a + b.
  129. * Using integer division, this is (2a + 1) * (2b + 1) / 2. Take
  130. * products of odd numbers, divide by two, mark that index as true.
  131. *
  132. * Then, we can check each true number for full validity by checking its
  133. * remainder mod 3: 2k + 5 is divisible by 3 if k is congruent to 2,
  134. * and (2 * (3a + 2) + 5) / 3 = 2a + 3 = 2 * (k / 3) + 3 (integer division)
  135. * = 2 * (k / 3 + 1) + 5 - 4 is not composite - known from whether
  136. * s[k / 3 + 1] is true.
  137. *
  138. * Let's verify this works. Row 1 is the index, row 2 its implied value,
  139. * row 3 a T if tagged as composite plus 4, and row 4 the index to check
  140. * to see whether s/3 is composite.
  141. *
  142. * 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  143. * 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45
  144. * T T T T T T T T
  145. * 1 2 3 4 5 6 7
  146. *
  147. * Row 4 aligns with a tagged value at the 6th index, so we set s[17] = s[6],
  148. * rendering it false. As we know, 39 is invalid. Looks good!
  149. *
  150. * But wait, we are going through and modifying data we might rely on
  151. * later - indeed, if we set false s[62] because 129, despite being
  152. * 125 + 4, is also 3 x 43, then we later reach 375, which is 7 x 53 + 4,
  153. * but also 125 x 3. By that point, 125 "looks" prime to our list, since
  154. * we marked 129 false, and so 375 gets wrongly marked false. However,
  155. * if we do the second step of the sieve in the reverse direction, every
  156. * odd prime can be looked up accurately.
  157. *
  158. * As a consequence, this method generates the sum list in O(S) space and
  159. * S log S time, where the log comes from multiplying all odd pairs that
  160. * have a product less than S: that's S / a operations for a up to sqrt(S)
  161. * which gives a harmonic series. Also, space is fantastically optimized;
  162. * to reach sums of N takes at most N/2 bits.
  163. *
  164. * ~~~~~~~~~ Beyond the Sum list ~~~~~~~~~~~
  165. *
  166. * But generating sums doesn't solve the problem, just makes it easier.
  167. * We still need to make relevant products. We can use a second bitset
  168. * to keep track of products and their uniqueness.
  169. *
  170. * For each sum, we tag all a * b with a + b = s. We also need to tag whether
  171. * that product has already been seen before. In the interest of time
  172. * efficiency, it would be nice to save for each product the corresponding
  173. * sum. However, that turns the space requirement from 2 bits per product
  174. * into 32 or more. Given that we want to monopolize system memory, we get
  175. * more out of it if we don't store the sum.
  176. *
  177. * Observe that since all the sums are odd, all pairs of numbers that add
  178. * to a given sum include one even number, so the products of interest are
  179. * all even. Since we are storing 2 bits apiece, we hold information for
  180. * the product 2k in indices 2k and 2k + 1 as follows:
  181. *
  182. * Unseen: 00
  183. * Seen: 01
  184. * Reseen: 11
  185. * updated by ps[2k + 1] = ps[2k], ps[2k] = true;
  186. *
  187. * Regarding efficiency of generating products:
  188. * We can make a * (s - a) without multiplying every time, as follows:
  189. * 3 * (s - 3) is 3 * s - 9, then add s - 7, s - 9, etc until adding 0.
  190. * Since we compute with s = 2 * a + 5, we can translate the procedure into
  191. * functions of a instead of s.
  192. *
  193. * We are generating a fixed-size lookup table, so, we only insert products
  194. * if we know we will get full info on them, at most NMIN * (S - NMIN), where
  195. * S is the largest sum value we allow.
  196. *
  197. * ~~~~~ "Now I do," "Now I do too" ~~~~~~~~
  198. *
  199. * We want to use these giant lookup tables to quickly generate solutions.
  200. * For a given sum, we iterate over all products as before. If the product
  201. * is tagged with 01, then we know:
  202. *
  203. * It has exactly one factorization for which the factors sum to an element of
  204. * the sum list - the element that brought us there.
  205. *
  206. * Because it is generated from the sum list, we know it has at least one other
  207. * factorization that does not sum to an element of the sum list.
  208. *
  209. * Thus, for this product, Polly can say "That was true, but now I do."
  210. *
  211. * There are three options when iterating over all products for a sum: either
  212. * we find 01 never, once, or more than once. If we find it exactly once, then
  213. * Sam can say "Now I do too." If we find it a second time, we can break early
  214. * from the loop. If we reach the end with exactly one 01, we have a pair of
  215. * numbers satisfying the problem statement, as desired. We save it to a set
  216. * for later processing of its properties.
  217. *
  218. * For further analysis, we provide a second operating mode that doesn't
  219. * display the solution pairs but instead displays for each sum the number of
  220. * times 01 is encountered. We call this the Freudenthal partition function
  221. * and its behavior is detailed in comethelper.cpp.
  222. *
  223. * ~~~~~~~~~~ Once space runs out ~~~~~~~~~~
  224. *
  225. * The above has an upper limit: encoutering products not saved in our lookup
  226. * table. This may happen at sums above 2 * sqrt(NMIN * (S - NMIN)). It
  227. * seems like a bit of a waste then, to have such a vast sum list but only
  228. * scan through O(sqrt(S)) of it.
  229. *
  230. * Fortunately, there are memory-efficient time-inefficient ways to check
  231. * pairs once we've exceeded this limit. This will be reasonably fast
  232. * at first, since it can reject any sum which finds two uniquely expressed
  233. * products, and it begins with the lowest products, which are still in the
  234. * lookup table. But in general, once we exceed our limit, we need to do
  235. * explicit factoring. We implement a slow p check and s check that will work
  236. * for any number. We scan factors of p. Every sum: if it is within our sum
  237. * lookup table, we check the table, otherwise, we factor s - 4 and s / 3 if
  238. * applicable. If we find exactly one good sum, we return as though p was
  239. * tagged 01. It just potentially takes O(p) time instead of constant.
  240. *
  241. * We use this p check to fill out the rest of the sum list. But since we
  242. * make our S list humongous and we don't want to wait for completion, we
  243. * allow SIGINT to kill this function without killing the program. After
  244. * SIGINT, the program dumps the sorted list of all solutions to cout, clears
  245. * memory, and cleanly exits.
  246. *
  247. * ~~~~~~~~ Complexity analysis ~~~~~~~~~~~~
  248. *
  249. * At the present version, we have a worst-case complexity analysis below:
  250. *
  251. * space: O(S), allocated in advance (no mid-program segfaults please)
  252. *
  253. * Generating our s list:
  254. * sieve1: every odd N < S contributes S / N -> S log S
  255. * sieve2: every 3rd element contributes constant -> S
  256. * Generating products:
  257. * genprod: every N < S contributes min(N, S / N) -> S log S
  258. * gensols: every N < sqrt(S) contributes at most N -> S
  259. * gensolsslow: N > sqrt(S) may contribute S / N constant
  260. * and (N - sqrt(S))^2 calls to slowPcheck
  261. * with argument between S and N^2
  262. * slowPcheck: may take arg / S calls to slowScheck
  263. * with argument between S and arg
  264. * slowScheck: may take sqrt(arg) time
  265. *
  266. * Thus, we generate solutions with sum up to order sqrt(S) in S log S time,
  267. * the next N take roughly N^3 sqrt(S) time (valid only for N << sqrt(S)),
  268. * and to generate all solutions with sum up to O(S), the form becomes S^5.
  269. *
  270. * Overnight, with space as large as my laptop could handle, this generated
  271. * 3141 solutions in FREUDENTHAL^=(3,*), with sum up to roughly 2 ^ 19.
  272. *
  273. * ~~~~~~~~~~~~~~~ NMIN > 3 ~~~~~~~~~~~~~~~~
  274. * It is not yet known whether solutions exist with NMIN > 3. The paper I
  275. * referenced above found none with sum less than 50000 and conjectured none
  276. * exist. We can follow the same procedure in general, but the S list must be
  277. * modified.
  278. *
  279. * Call a number pseudoprime in our case if it has no factorizations with at
  280. * least one factor greater than or equal to NMIN. For NMIN = 3, that makes
  281. * 4 the only nonprime pseudoprime. For NMIN = 4, the list is 4, 6, 9. For
  282. * NMIN = 5, we get 6, 8, 9. The largest even pseudoprime for a given NMIN can
  283. * be shown to be 2 * (NMIN - 1). For a pseudoprime p and a prime q, p * q has
  284. * only one allowable factorization and so p + q is excluded from S. The sum
  285. * of two pseudoprimes is a weird case. It is disallowed if there is no way
  286. * to legally exchange factors. For NMIN = 6 we have 6 and 8 as pseudoprimes,
  287. * and 48 factors uniquely, but 15 is also a pseudoprime and 8 * 15 = 120 is
  288. * expressible as 6 * 20. I don't see a particularly systematic way to find
  289. * this, so we can just check every option.
  290. *
  291. * We only need to look at even primes and pseudoprimes in the NMIN = 2 and 3
  292. * case because we were safe to assume that every even number at least 6 is the
  293. * sum of two odd primes. But it could be the case that a given number is not
  294. * the sum of two odd primes both of which are at least NMIN. For example,
  295. * if NMIN = 6, then 16 is out because both its partitions rely on smaller
  296. * primes.
  297. *
  298. * A small list of largeish even numbers made valid by raising NMIN:
  299. * NMIN = 14, 44 becomes valid. With NMIN = 14, pseudoprimes include 169.
  300. * NMIN = 32, 92 becomes valid. (pseudoprimes up to 961)
  301. * NMIN = 104, 242 becomes valid.
  302. *
  303. * From muddling about a bit with goldbach partitions, it seems like a very
  304. * safe assumption that if we check all possible sums of pseudoprimes then we
  305. * can assume any even number beyond that range is the sum of two valid primes.
  306. *
  307. * ~~~~~ Respecialization to NMIN = 4 ~~~~~~
  308. *
  309. * With NMIN = 4, we have pseudoprimes 4, 6, and 9. By hand we check that
  310. * 16 = 4 x 4 eliminates 8, 25 = 5 x 5 eliminates 10, 35 eliminates 12,
  311. * 49 eliminates 14, 55 eliminates 16, and 65 eliminates 18. Our sum list
  312. * consists of odd numbers for which s - 4 and s - 6 are both composite
  313. * and which are not thrice a prime. This is just one further condition
  314. * than we encountered for NMIN = 3, and can be implemented by altering
  315. * the way sieve2 works to also window every number with its neighbor.
  316. *
  317. * If we want to be able to use our sieve trick to eliminate numbers of the
  318. * form 3 * p, then we need to work backwards on our S list. I'm debating
  319. * whether we want to have s[k] hold 2k + 5 or 2 * (k + NMIN) - 1 in general.
  320. * If NMIN = 2 or 3 the latter is ideal. If NMIN = 4 or 5, then after the
  321. * sieve1 we need to have s[k]_post = (s[k] & s[k + 1])_pre. This is difficult
  322. * to accomplish if we work backwards, but I guess we just store a value before
  323. * overwriting.
  324. *
  325. * The benefit of keeping 2k + 5 in the specialized case of NMIN = 4 is that
  326. * we window in the appropriate direction for a reverse scan. All unscanned
  327. * indices hold relevant information about primality. But generalizing with
  328. * this goal would mean that for an arbitrary NMIN we would have s[k] storing
  329. * 2k + ((NMIN + 1) & ~1) + 1. That makes prime lookups awkward, but
  330. * if we specialize above NMIN = 5 it means we don't need to think about
  331. * storing arbitrary overhead information for our sieve.
  332. *
  333. * Let's try this latter goal, and look at a few sum lists in the making:
  334. *
  335. * NMIN = 2: tag composites plus 2 at index (c / 2).
  336. * 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  337. * 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43
  338. * T T T T T T T T
  339. *
  340. * NMIN = 3: tag composites plus 4 at index (c / 2)
  341. * 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  342. * 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45
  343. * T T T T T T T T
  344. * Then, as above, for any index congruent to 2 mod 3, look at (k + 1) / 3 for
  345. * information about the primality of s / 3.
  346. *
  347. * NMIN = 4: tag composites plus 4 at index (c / 2)
  348. * 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  349. * 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45
  350. * T T T T T T T T
  351. * Then, look at k - 1 and, for any index congruent to 2 mod 3, (k + 1) / 3.
  352. * Perform bitwise and with both.
  353. *
  354. * NMIN = 5: tag composites plus 6 at index (c / 2)
  355. * 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  356. * 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47
  357. * T T T T T T T T
  358. * Look at k - 1 and, for any index congruent to 1 mod 3, (k + 2) / 3.
  359. * Perform bitwise and with both.
  360. *
  361. * NMIN = 6: tag composites plus 6 at index (c / 2)
  362. * 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  363. * 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47
  364. * T T T T T T T T
  365. * Look at k - 1 and k - 2, and for index congruent to 1 mod 3, (k + 2) / 3.
  366. * Perform bitwise and with both.
  367. *
  368. * NMIN = 7: tag composites plus 8 at index (c / 2)
  369. * 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  370. * 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
  371. * T T T T T T T T
  372. * Look at k - 1 and k - 2, and for index congruent to 0 mod 3, (k + 3) / 3.
  373. * etc..
  374. *
  375. * General: tag composites plus (NMIN + 1) & ~1 at index (c / 2)
  376. * Look at (distances to other even pseudoprimes, all backward), and for index
  377. * with k + (NMIN - 1) / 2 divisible by three, look at k / 3 + (NMIN - 1) / 6.
  378. *
  379. * Then when creating products, we need to begin with NMIN * (s - NMIN), which
  380. * is NMIN * (2 * k + 1 + ((NMIN + 1) & ~1) - NMIN), which is given most simply
  381. * by NMIN * (2 * k + (NMIN % 2) + 1). Example, with NMIN = 7, k = 14, we
  382. * have 7 * (2 * 14 + 2) = 210, which is the smallest number with a factor sum
  383. * of 37, which is what s[14] holds for NMIN = 7.
  384. *
  385. * Incrementing from there we need (NMIN + 1) * (s - NMIN - 1), which is the
  386. * prior value plus s - 2 * NMIN - 1. Again using our expression for s from k,
  387. * we get 2 * k - NMIN - ((NMIN % 2) ? 1 : 0). These ternaries are pretty ugly
  388. * but fortunately the precompiler is the one evaluating them so our code will
  389. * still perform well. Actually since this expression is always even, we can
  390. * use integer division to write 2 * (k - NMIN / 2). This value on the right
  391. * decreases by one each iteration to zero to list off the products. Example:
  392. * 2 * (14 - 3) = 22 adds to 210 to give us 232, which is 8 x 29.
  393. *
  394. * We ignore for the moment the potential for even numbers and sums of
  395. * pseudoprimes to be allowed sums and just focus on implementing this general
  396. * algorithm. Presently, the code functions by assuming the only concerning
  397. * pseudoprimes are even numbers between NMIN and 2*(NMIN - 1). I have yet to
  398. * find a case where this assumption affects the outcome.
  399. */
  400.  
  401. #include <iostream>
  402. #include <cstdint> // For definite bit widths
  403. #include <set>
  404. #include <bitset>
  405. #include <cmath>
  406. #include <csignal>
  407.  
  408. // We have found S can be around 16 billion / NMIN before memory restrictions
  409. #define NMIN 3
  410. #define S 409ull
  411.  
  412. // The VERBOSE option allows progress updates to be sent to std::cerr.
  413. #define VERBOSE true
  414.  
  415. /* We provide two modes. If COMET is set to false, then std::cout recieves
  416. * solution pairs. If COMET is set to true, then std::cout instead receives
  417. * ordered pairs of the form s, f(s) where f(s) is the Freudenthal partition
  418. * function (analogous to the Goldbach partition function), counting the
  419. * number of products which are uniquely fixed by Sam's claim.
  420. */
  421. #define COMET true
  422.  
  423. /* Function nums : generates the original numbers given their sum and product
  424. * arguments: sum and product
  425. * return : packaged form of number pairs. If a + b = s, ab = p, a is even,
  426. * then bits 0 - 31 hold b and bits 32 - 63 hold a. This allows us
  427. * to sort the outputs by their even number, breaking ties by the
  428. * odd number.
  429. */
  430. uint64_t nums(uint64_t s, uint64_t p)
  431. {
  432. int64_t b(sqrt(s * s - 4 * p));
  433. if ((s + b) % 4)
  434. b = -b;
  435. return ((s + b) << 31) + ((s - b) >> 1);
  436. }
  437.  
  438. /* Section headers of the top comment are reproduced below to give reference
  439. * for a more detailed discussion of functions in the section.
  440. */
  441.  
  442. // ~~~~~ What can Sam's sum look like? ~~~~~
  443.  
  444. /* Function sieve1 : implements an erastosthenes sieve on odd numbers.
  445. * argument: pointer to bitset with all bits off
  446. * post : bits for which 2 * a + 1 is composite are set
  447. */
  448. void sieve1(std::bitset<S/2> * s)
  449. {
  450. uint64_t b;
  451. for (uint64_t a = 3;; a += 2)
  452. {
  453. b = a * a / 2;
  454. if (b > s->size())
  455. break;
  456. for (; b <= s->size(); b += a)
  457. (*s)[b] = true;
  458. }
  459. }
  460.  
  461. /* Function window : A helper function for sieve2 with NMIN > 3.
  462. * note : I'm really hoping -O3 optimizes this whole function out and just
  463. * replaces it with an unrolled loop with short-circuited booleans.
  464. * Particularly for NMIN = 3, the loop doesn't acually run so the relevant
  465. * machine code is far simpler if optimized.
  466. */
  467.  
  468. inline void window(std::bitset<S/2> * s, uint64_t index)
  469. {
  470. for (int a = 1; a < NMIN / 2; ++a)
  471. (*s)[index] = (*s)[index] && (*s)[index - a];
  472. }
  473.  
  474. /* Function sieve2 : Performs a reverse scan to generate sum lookup table
  475. * argument: pointer to bitset with composite numbers set
  476. * post : bitset contains all possible sums for which no pair of summands
  477. * may be uniquely determined from their product.
  478. */
  479. void sieve2(std::bitset<S/2> * s)
  480. {
  481. // We have two goals: window every index and take care of the 3p case.
  482. // So that our loop may have definite structure, we enforce that b always
  483. // holds an index encoding a multiple of 3 and a holds the index holding
  484. // the primality of the encoded sum divided by 3.
  485.  
  486. // The highest b in our array is one of the top 3 indices, satisfying
  487. // b + (NMIN - 1) / 2 == 3 * a.
  488. uint64_t b = S / 2 - 1;
  489. if ((b + (NMIN - 1) / 2) % 3 > 1)
  490. window(s, b--);
  491. if ((b + (NMIN - 1) / 2) % 3 > 0)
  492. window(s, b--);
  493.  
  494. uint64_t a = (b + (NMIN - 1) / 2) / 3;
  495. while (b > 2)
  496. {
  497. (*s)[b] = (*s)[b] && (*s)[a];
  498. window(s, b--);
  499. window(s, b--);
  500. window(s, b--);
  501. --a;
  502. }
  503.  
  504. /* Displays the sum list to cout: a nice sanity check when debugging.
  505. */
  506. }
  507. void writeS(std::bitset<S/2> * s)
  508. {
  509. for (int i = 0; i < S/2; ++i)
  510. if ((*s)[i])
  511. std::cout << (2 * i + ((NMIN + 1) | 1)) << ", ";
  512. std::cout << std::endl;
  513. }
  514.  
  515. // ~~~~~~~~~ Beyond the Sum list ~~~~~~~~~~~
  516.  
  517. /* Function genprod : tags products according to how many sums generate them
  518. * arguments: pointer to filled sum bitset and empty product bitset
  519. * post : product bitset is tagged appropriately
  520. */
  521. void genprod(std::bitset<S/2> * s, std::bitset<NMIN*S> * ps)
  522. {
  523. uint64_t c, d;
  524. uint64_t cmax = NMIN * (S - NMIN);
  525. for (uint64_t a = 0; a < S / 2;)
  526. if ((*s)[a++])
  527. for (d = a - NMIN / 2, c = NMIN * (2 * a - 1 + NMIN % 2);
  528. d && c <= cmax; c += (--d << 1))
  529. {
  530. (*ps)[c + 1] = (*ps)[c];
  531. (*ps)[c] = true;
  532. }
  533. }
  534.  
  535. // ~~~~~ "Now I do," "Now I do too" ~~~~~~~~
  536.  
  537. /* Function gensols : scans sum and product tables for solutions
  538. * arguments: pointers to sum and product tables, reference to solution set
  539. * post : solution set is filled to sums at most 2 * sqrt(NMIN * (S - NMIN))
  540. */
  541. void gensols(std::bitset<S/2> * s, std::bitset<NMIN*S> * ps,
  542. std::set<uint64_t> &solutions)
  543. {
  544. uint64_t c, d;
  545. const uint64_t smax = sqrt(NMIN * (S - NMIN));
  546. uint64_t uniqc;
  547. for (uint64_t a = 0; a < smax - NMIN;)
  548. if ((*s)[a++])
  549. {
  550. uniqc = 0; // The solution product, or f(s) in COMET mode.
  551. for (d = a - NMIN / 2, c = NMIN * (2 * a - 1 + NMIN % 2);
  552. d; c += (--d << 1))
  553. if (!(*ps)[c + 1])
  554. {
  555. if (COMET)
  556. ++uniqc;
  557. else
  558. {
  559. if (uniqc)
  560. break;
  561. else
  562. uniqc = c;
  563. }
  564. }
  565. if (COMET)
  566. std::cout << 2 * a + ((NMIN - 1) | 1) << ", "
  567. << uniqc << std::endl;
  568. else if (!d && uniqc)
  569. solutions.insert(nums(2 * a + ((NMIN - 1) | 1), uniqc));
  570. }
  571. }
  572.  
  573. //~~~~~~~~~ Once space runs out ~~~~~~~~~~~
  574.  
  575. /* Function isOddPrime : helper to slowScheck
  576. * argument: an odd number
  577. * return : true if prime, false otherwise
  578. */
  579. bool isOddPrime(uint64_t oddo)
  580. {
  581. uint64_t amax = sqrt(oddo);
  582. for (uint64_t a = 3; a <= amax; ++a)
  583. if (!(oddo % a))
  584. return false;
  585. return true;
  586. }
  587.  
  588. /* Function slowScheck : determines if a number is a valid sum
  589. * argument: an odd number
  590. * return : true if valid (emulates (*s)[sum / 2 - 2])
  591. */
  592. bool slowScheck(uint64_t sum)
  593. {
  594. if (NMIN == 2)
  595. return !isOddPrime(sum - 2);
  596. if (sum % 3 == 0 && isOddPrime(sum / 3))
  597. return false;
  598. for (int a = 2 * (NMIN - 1); a >= NMIN; a -= 2)
  599. if (isOddPrime(sum - a))
  600. return false;
  601. return true;
  602. }
  603.  
  604. /* Function slowPcheck: Does p factor into any valid sum besides the one given?
  605. * arguments: a product and sum pair, and pointer to the known valid sums.
  606. * return : true if another pair is found, false if p is uniquely determined
  607. * (emulates (*ps)[p + 1])
  608. * note : first considers factors with minimum sum to minimize chance of
  609. * needing to call slowScheck.
  610. */
  611. bool slowPcheck(uint64_t p, uint64_t s1, std::bitset<S/2> * s)
  612. {
  613. uint64_t sum;
  614. for (uint64_t a = sqrt(p); a >= NMIN; --a)
  615. // To proceed, p must divide a, the factor sum must be odd and not s1,
  616. // and if all those succeed, we can check whether the sum is a valid
  617. // s number. If it is, then we know p isn't uniquely determined by
  618. // the provided info.
  619. if (!(p % a) && ((sum = a + p / a) % 2) && sum != s1
  620. && (sum < S && (*s)[(sum - NMIN - 1) / 2]
  621. || sum >= S && slowScheck(sum)))
  622. return true;
  623. return false;
  624. }
  625.  
  626. /* Function gensolsslow : continues where gensols left off but slower
  627. * arguments: same as gensols
  628. * note : as it runs long, we permit it to be interrupted by SIGINT
  629. * return : the largest sum that was checked to completion
  630. *
  631. * aside : In solution mode, it quits after finding a single point of
  632. * inconsistency with Sam's 3rd claim, which happens very quickly
  633. * for NMIN > 3. Thus, we are able to verify the absence of
  634. * solutions up to a very significant bound on the sum. In
  635. * comet mode, we must check every product, so this is far slower.
  636. */
  637. volatile sig_atomic_t flag = 1;
  638. void exitSlowFunction(int sig)
  639. {
  640. flag = 0;
  641. }
  642. uint64_t gensolsslow(std::bitset<S/2> * s, std::bitset<NMIN*S> * ps,
  643. std::set<uint64_t> &solutions)
  644. {
  645. uint64_t c, d;
  646. uint64_t cmax = NMIN * (S - NMIN);
  647. uint64_t a = sqrt(cmax) - NMIN;
  648. uint64_t uniqc;
  649. while (flag && a < S/2)
  650. if ((*s)[a++])
  651. {
  652. uniqc = 0;
  653. for (d = a - NMIN / 2, c = NMIN * (2 * a - 1 + NMIN % 2);
  654. d; c += (--d << 1))
  655. if (c <= cmax && !(*ps)[c + 1]
  656. || (c > cmax && !slowPcheck(c, 2 * a + ((NMIN - 1) | 1), s)))
  657. {
  658. if (COMET)
  659. ++uniqc;
  660. else
  661. {
  662. if (uniqc)
  663. break;
  664. else
  665. uniqc = c;
  666. }
  667. }
  668. if (COMET)
  669. std::cout << 2 * a + ((NMIN - 1) | 1) << ", "
  670. << uniqc << std::endl;
  671. else if (!d && uniqc)
  672. {
  673. solutions.insert(nums(2 * a + ((NMIN - 1) | 1), uniqc));
  674. if (NMIN > 3)
  675. std::cerr << "Counterexample found!!" << std::endl;
  676. }
  677. }
  678. return 2 * a + ((NMIN - 1) | 1);
  679. }
  680.  
  681. /* Procedure that utilizes all the above.
  682. */
  683. int main()
  684. {
  685. signal(SIGINT, exitSlowFunction);
  686.  
  687. std::bitset<S/2> * s = new std::bitset<S/2>;
  688. std::bitset<NMIN*S> * ps = new std::bitset<NMIN*S>;
  689.  
  690. if (VERBOSE)
  691. std::cerr << "prepared memory successfully" << std::endl;
  692.  
  693. sieve1(s);
  694. if (NMIN > 2)
  695. sieve2(s);
  696.  
  697. if (VERBOSE)
  698. std::cerr << "sieving complete, S list generated" << std::endl;
  699.  
  700. genprod(s, ps);
  701.  
  702. if (VERBOSE)
  703. std::cerr << "tagged products up to " << ps->size() << std::endl;
  704.  
  705. std::set<uint64_t> solutions;
  706. gensols(s, ps, solutions);
  707.  
  708. if (VERBOSE)
  709. std::cerr << solutions.size()
  710. << " solutions with sum at most "
  711. << (uint64_t) (2 * sqrt(NMIN * (S - NMIN))) << std::endl;
  712.  
  713. if (VERBOSE && !COMET)
  714. std::cerr << "Product table exhausted, further solutions will be "
  715. << "generated far more slowly" << std::endl;
  716.  
  717. uint64_t maxsum = gensolsslow(s, ps, solutions);
  718.  
  719. if (VERBOSE && !COMET)
  720. std::cerr << ((flag) ? "Maxed out sum table" : "SIGINT received")
  721. << " at a sum of " << maxsum
  722. << ", writing " << solutions.size()
  723. << " solutions" << std::endl;
  724.  
  725. if (VERBOSE && COMET)
  726. std::cerr << "Reached sum of " << maxsum;
  727.  
  728.  
  729. // Clear memory and present solution pairs.
  730. delete s;
  731. delete ps;
  732. if (!COMET)
  733. for (const auto& i : solutions)
  734. std::cout << (i >> 32) << ", " << (i & 0xffffffff) << std::endl;
  735.  
  736. return 0;
  737. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement