# Untitled

a guest
May 12th, 2018
633
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
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
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. *
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. }