Erdős problem #647 — GPU search (no solution with 24 < n ≤ 10^22)

Problem (Erdős–Selfridge, £25). Is there an n > 24 with max_{m<n} (m + τ(m)) ≤ n + 2? Equivalently C(n): τ(n−k) ≤ k+2 for all 1 ≤ k < n. Known solutions: 2, 3, 4, 5, 6, 8, 10, 12, 24 (OEIS A087280).

Prior computation (erdosproblems.com comments, Sep 2026): full sieve to 10^12 (Idén), Lean certificate to 10^9 (Mian–Siddique, arXiv:2608.17880), residue-class searches to ≈ 6.2×10^17 (Hughes) and reportedly ≈ 9.2×10^18.

Files

file	purpose
search.py	the pipeline: base classes → wheel → GPU stage 1 (AP sieve) → GPU stage 2 (exact-rule scan) → CPU stage 3 (sympy). Checkpoints to ckpt_<tag>_d<delta>.json.
exact.py	exact reference: divisor sieve + running maximum; lists solutions of the relaxed condition τ(n−k) ≤ k+2+δ.
reduce.py, wheel.py	earlier standalone derivations (superseded by search.py, kept for reference).
main.log	log of the production run (--delta 0 --nmax 1e20).

Method

Only two facts about τ are used: d | x ⇒ τ(d) ≤ τ(x) and multiplicativity. Every stage is a kill rule that discards n only when τ(n−k) > k+2+δ is proven for some k; whatever survives goes to the next stage.

1. Brute force to BF_LIMIT = 4·10^6 with a divisor sieve (exact).
2. Base classes. For n ≡ r (mod M), M = ∏ p^e, the p-adic valuation of n−k is known exactly when it is below the cap e. T = ∏(v_p+1) is a lower bound for τ(n−k); if all valuations are exact then n−k = D·c with gcd(c, M) = 1, so τ(n−k) ≥ 2T unless n = D + k (recorded and checked directly). Lifting through 2^6 3^4 5^3 7^3 shows: every solution n > 40 is a multiple of 2520 (exceptions ≤ 24).
3. Killer shifts. For n ≡ 0 (mod 2520) let T0(k) be the bound from the 2520-part of n−k (caps 2^3, 3^2, 5, 7). For a prime p > 7 with v_p(n−k) = 1 we get τ(n−k) ≥ 4·T0(k) unless the cofactor is 1, which only happens for n < BF_LIMIT (bound computed and asserted in code). Shifts with 4·T0(k) > k+2 are killers: k ∈ {1,2,3,4,5,6,8,9,10,12,18,20,24}; those with 3·T0(k) > k+2 are strong (kill p² | n−k as well).
4. Wheel. Residues of n modulo W = 2520·11²·13²·17²·19² ≈ 5.4×10^12 surviving the killers. A p² | n−k sub-residue is normally kept (conservative), but the only way it can survive is the cofactor-1 family n−k = Du·∏q^w·p^v (v ≥ 2); all such powers are periodic mod 2520 once past the caps, so a finite enumeration (family_possible) decides whether the family can meet n ≡ r0 (mod 2520) at all. When it cannot, the sub-residue is dead too (strong_for). Result: 2,574,990 residues, density 4.8×10^−7 (4,140,360 without the family check).
5. GPU stage 1 (stage1 kernel, one thread per candidate residue, J consecutive multiples of W per chunk, n = R + jW as a 128-bit number): two further wheel primes (23, 29 by default, --extra) are applied on the fly by CRT inside the kernel, so the 2.57M stored residues become 2.57M·253·… candidates without materialising them. For primes 31 ≤ p ≤ 1000 and each killer k, the dead j form an arithmetic progression mod p (minus the sub-progression mod p² unless k is strong for that p, decided per prime by strong_for). Bitmask in shared memory in a conflict-free transposed layout; per-shift offsets come from a precomputed k·W⁻¹ mod p table and the sub-progression test is an add-and-compare counter. ~9×10^−5 of the candidates survive.
6. GPU stage 2 (stage2 kernel): survivors are 128-bit; for every shift k ≤ 64 and every prime the exact valuation is computed, T accumulated, and the candidate is dead if T > k+2 or if 2T > k+2 and n−k ≠ D (cofactor > 1). Two passes: primes ≤ 8192 first, the full list to 65536 only on what survives. Nothing survives.
7. CPU stage 3: sympy factorisation, shifts up to 3000; anything surviving is printed as *** STAGE-3 SURVIVOR.

Validation

--delta δ relaxes the budget to k+2+δ everywhere, which creates many solutions on small ranges, and --test compares the pipeline's final solution set with the exact sieve. With --bf 200000 the solutions above 2·10^5 must be recovered by the GPU path:

• δ=6, nmax=2·10^8: 211 solutions, pipeline = exact.
• δ=12, nmax=2·10^8: 2294 solutions, 868 found by the GPU path (up to n = 36,128,400), pipeline = exact.
• test_stage2.py: the stage-2 kernel on 3000 synthetic candidates above 2^64 (n ≈ 2–10×10^19) agrees verdict-for-verdict with a pure-Python replica of the rule (δ=60: 2258 alive; δ=6: 7 alive) — exercises the 128-bit limb arithmetic.
• exact.py 1e9 0: independent exact sieve reproduces OEIS A087280 (only 2,3,4,5,6,8,10,12,24 up to 10^9).

Results

• 2026-09-10, 10^20: --nmax 1e20 --tag main, first kernel version, 6814 s (1 h 54 min): no solution with 24 < n ≤ 10^20 (main.log, ckpt_main_d0.json; 4542 chunks, 1.92×10^9 stage-1 survivors, 0 stage-2 survivors).

• 2026-09-10, 10^21: --nmax 1e21 --tag run21, optimised kernel, 5009 s (1 h 23 min): no solution with 24 < n ≤ 10^21 (run21.log, ckpt_run21_d0.json; 465 launches, 1.19×10^10 stage-1 survivors, 0 stage-2 survivors).

• 2026-09-11 to 2026-09-13, 10^22: --nmax 1e22 --tag run22, same kernel, 71,050 s (19.7 h) of GPU time over four sessions with --resume (GPU power limit 280–350 W): no solution with 24 < n ≤ 10^22 (run22.log, ckpt_run22_d0.json; 3 chunks × 465 launches, 1.64×10^11 stage-1 survivors, 0 stage-2 survivors).

Combined with the brute force to 4·10^6 this independently verifies that the only solutions up to 10^22 are 2, 3, 4, 5, 6, 8, 10, 12, 24, extending the published bounds (10^12 full sieve; ≈ 6.2×10^17 and ≈ 9.2×10^18 by residue searches). This is an uncertified computation (no Lean certificate).

Throughput (RTX 3090)

version	candidates per 10^21	rate	10^21	10^22
kernel used for the 10^20 run	7.7×10^14	1.1×10^10 /s	19 h	8 days
conflict-free bitmask, family check	4.8×10^14	4.4×10^10 /s	3 h	30 h
+ 23, 29 on the fly, offset tables, two-pass stage 2	1.3×10^14	4.5×10^10 /s	0.8 h	8 h

Every kernel version was re-validated with the δ-tests and the stage-2 unit test. (The rewrite briefly emitted n0+j instead of n0+jW for survivors; the δ=12 test caught it as four missing solutions.)

Reproduce

python search.py --delta 0 --nmax 1e22 --tag run22          # fresh run to 10^22 (covers 10^20 again in minutes); --resume continues a run with the same geometry
python search.py --delta 12 --nmax 2e8 --wheel 11,13 --extra "" --B1 200 --bf 200000 --test   # validation
python search.py --delta 6 --nmax 2e8 --wheel 11,13 --extra 23 --B1 200 --bf 200000 --test      # validation (CRT path)
python test_stage2.py 60                                                                        # stage-2 kernel vs Python replica

Environment: Python 3.12, cupy-cuda13x 14.2, numba, sympy; CUDA 13.3 driver, no host C++ compiler needed (NVRTC).