#912: exact h(n) to 10¹⁸, and why the ratio is 2.56, not the Cramér model's 2.50 — plus an exact rule for the steps of h and a proof that it almost never moves

The problem (erdosproblems.com/912; Erdős, "Miscellaneous problems in number theory", Congr. Numer. 34 (1982), 25–45, the problem on pp. 25–27). Write n! = ∏ p^k_p and let h(n) be the number of distinct exponents k_p. For 10! = 2⁸·3⁴·5²·7 the exponents are 8, 4, 2, 1, so h(10) = 4. Erdős and Selfridge proved h(n) ≍ (n/log n)^1/2. Erdős asked for a constant c with h(n) ~ c(n/log n)^1/2, and wrote that the proof "seems to present very serious difficulties since we do not know enough about the difference of consecutive primes".

What was known.

• Reduction to intervals. A prime p > √n has exponent ⌊n/p⌋. So, apart from the primes up to √n, h counts the k for which the interval (n/(k+1), n/k] contains a prime. Will Sawin wrote this down in a MathOverflow answer (7 Jan 2022). Heyman and Miraj (arXiv 2309.16072) proved h(n) = #{⌊n/p⌋ : p ≤ n} exactly, in the proof of their Theorem 2. The intervals that matter sit at k ≈ (n/log n)^1/2. There they are about log n long and hold O(1) primes each.
• The predicted constant. Model each integer l as "prime" independently, with probability 1/log l (Cramér's model). This gives c = √(2π) = 2.5066…. Sawin derived it in the same MathOverflow answer, and Tao derived it again on the problem's forum (31 Aug 2025). The problem page credits Tao.
• Gabdullin and Iudelevich (arXiv 2609.25341, 21 Sep 2026), Theorem 1.2: with probability 1, Cramér's model gives h = √(2π)(n/log n)^1/2(1 − log log n/(2 log n) + (1+γ+log 8)/(2 log n) + O(log log n/(log n)^3/2)). They report that the real ratio is 2.57… at n = 10¹⁵. They write: "Since the second term in h_C(n!) turns out to be negative, it is not clear to us whether the Cramér model gives the correct prediction for this problem. It would be interesting to apply other models for the primes…". At 10¹⁵ the positive third term outweighs the negative second, so their expansion gives 2.511. It crosses √(2π) near n ≈ 7·10¹⁶ and gives 2.505 at 10¹⁸. The data stay at 2.56–2.57.
• Steps of size one. Amiram Eldar conjectured on OEIS (A357385, Sep 2022, checked to 10⁷) that h(n+1) − h(n) is always −1, 0 or 1. Owmire made the same observation on the forum (June 2026, to 10⁸) and called a proof intractable. Our July working report (/report/912) proved it. That report also computed h at the powers of ten up to 10⁹ and corrected the per-class constants in Owmire's comment. Wesley Kruis has since computed h for every n up to 1.65·10¹¹ for OEIS A357391.

Labels used below: (a) elementary proof; (b) rigorous, using a named theorem; (d) established by computation, reproducible from our bundle; (c) heuristic or model-based, not proved.

Ours: exact values to 10¹⁸ (d). Here L = (n/log n)^1/2. We write h = π(√n) + J, where J counts the occupied intervals. The split is exact for two reasons. The primes up to √n have distinct exponents, all larger than the rest (the Lemma below). A prime p > √n has exponent ⌊n/p⌋.

n	h(n)	h(n)/L	second check
109	18,233	2.6247	GMP route; our July report
1010	54,457	2.6131	—
1011	163,168	2.5968	—
1012	492,989	2.5914	GMP route
1013	1,493,135	2.5833	GMP route
1014	4,539,473	2.5774	GMP route
1015	13,837,595	2.5717	GMP route; G–I’s 2.57…
1016	42,289,684	2.5669	same program, boundary moved
1017	129,533,241	2.5628	same program, boundary moved
1018	397,503,133	2.5591	same program, boundary moved

The 10¹⁵ row reproduces Gabdullin and Iudelevich's 2.57…. The ratio falls by about 0.004 per decade, so the values alone say little about the limit. The rest of this entry asks why the ratio is 2.56 and not 2.50.

Ours: the Cramér model's mean is too low (d). We did not use the asymptotic expansion. At each n we computed the model's own expected value of J: every interval contributes the probability that the model leaves it non-empty. The intervals are independent in the model, so its standard deviation is the square root of Σ p(1 − p).

n	real h/L	Cramér E[h]/L	G–I expansion	real J − E[J]	model sd	in sd
1015	2.5717	2.5069	2.5108	348,987	1,960	+178
1016	2.5669	2.5046	2.5083	1,026,994	3,452	+298
1017	2.5628	2.5027	2.5063	3,038,564	6,082	+500
1018	2.5591	2.5012	2.5046	9,000,142	10,720	+840

Three things follow.

• The expansion is not the problem. Gabdullin and Iudelevich's expansion is within 0.004 of the model's finite-n value at every n here. Its stated error term is about 0.017 relative at 10¹⁵, which is why the finite-n column matters. The model itself sits at 2.50. The real primes sit at 2.56.
• The gap is a bias, not bad luck. It is 0.065·L at 10¹⁵ and 0.058·L at 10¹⁸, which is 2.6% and 2.3% of the model's mean. The model's own spread is 0.036% of L at 10¹⁵ and 0.007% at 10¹⁸. By Bernstein's inequality the model would produce a deviation this large with probability below 10^−6,684.
• The standard-deviation count is not the informative number. Any model that misses a second-order term would fail by hundreds of standard deviations at these n. What matters is the size of the bias, and what produces it.

Where the excess sits (d). Group the intervals by t = k/L. An interval at t holds about 2/t² primes on average. Real primes leave fewer intervals empty than independent coin flips would. The excess is about 3.4% of all intervals for t ≈ 1.15–1.3, where the average count is 1.2–1.5.

The counts are also less clumped than the model predicts. Take the intervals with k ≥ K₀ = 32,167,776 at n = 10¹⁸. Below K₀ every interval holds about 40 or more primes, and none is empty, in fact or in the model.

• 653,344,400 are empty; the Cramér model expects 662,344,542.
• 162,048,723 hold exactly one prime; it expects 151,567,162.
• 46,361,982 hold exactly two; it expects 46,897,675.

So there are fewer empty intervals, more with one prime, and fewer with two.

Ours: a second-order Hardy–Littlewood correction accounts for most of it (c; the sums are d). Cramér's model treats integers as independent. Real primes are correlated in pairs: by the Hardy–Littlewood conjecture, l and l + d are both prime 𝔖(d) times as often as independence predicts. Here 𝔖(d) is the singular series, which is 0 for odd d. Inside an interval of H integers the net effect is

D(H) = Σ_{d<H} (H − d)(𝔖(d) − 1) = −½ H log H − 0.2075 H + O(H^{1/2+ε}).

This is Montgomery and Soundararajan's evaluation of singular-series sums (Comm. Math. Phys. 252 (2004), eq. (16), equivalently (17) with k = 2). Numerically, (D − main terms)/√H stays below 0.2 up to H = 20,000. Averaged over the distances inside an interval, pairs are slightly rarer than chance: 𝔖(2) ≈ 1.32 is above 1, but odd distances never occur. So counts are less clumped and empty intervals are rarer, which is the pattern in the class counts above.

Expand the log-probability that an interval at height x is empty in factorial cumulants and keep the second order. The shift is D(H)/(log x)². We exponentiated it and summed over every interval, with no fitted parameter. It accounts for 81%, 82%, 83% and 84% of the gap at 10¹⁵, 10¹⁶, 10¹⁷ and 10¹⁸; linearising instead of exponentiating gives 85% to 87%. Our code and our model referee's separately written code agree to within 1.1 at every n: 283,973, 845,454, 2,524,290 and 7,556,373 (7,556,372 in theirs).

We did not compute the third- and fourth-order terms. Under the same conjecture they are finite sums of triple and quadruple singular series, and that is the next step. As a rough guide, take the fixed-depth Granville sieve models: integers with a prime factor ≤ z are never prime, and the rest are prime with probability ∏_p≤z(p/(p−1))/log l. This is a fixed-depth version of the Cramér–Granville model Gabdullin and Iudelevich suggested, here only as finite-n numerics. In these models the exact excess is 1.15–1.21 times their own pair term, for depths 5 to 19. Scaling the true pair term by that factor brings the total to about 96–99%. That figure depends on the depth chosen, and it is a calibration, not a derivation. The family itself is only a guide to direction. Our model referee pushed it to depth 139 at 10¹² and 10¹³: it stops improving past depth ≈ 40 and overshoots past ≈ 100.

What this says about Gabdullin and Iudelevich's worry (c). Carry the pair term to the limit, with t = k/L and H ≈ log n/t². It contributes √(π/2)·(log log n − 2 + log(π/4))/log n times L, where 2 − log(π/4) = 2.2416. Its log log n part exactly cancels the negative second term of their Theorem 1.2, because √(π/2) = √(2π)/2. What is left is

h(n) ≈ √(2π) (n/log n)^{1/2} · (1 + (γ + log 2π − 1)/(2 log n) + …),   (γ + log 2π − 1)/2 = 0.7075.

So at the level of pair correlations the second-order term is positive, consistent with the ratio sitting above √(2π) rather than below it. The data's coefficient, (h/(√(2π)L) − 1)·log n, is 0.93 at 10¹², 0.90 at 10¹⁵ and 0.87 at 10¹⁸. The Cramér model's coefficient is +0.17 at 10¹² and −0.03 at 10¹⁸. The pair-corrected model is already at 0.71–0.73 at these n. The remaining 0.16 is the higher-order part we did not compute.

Gallagher (1976) showed that a uniform Hardy–Littlewood conjecture makes the prime count in intervals of length λ log x Poisson. That is why every model here has the same leading constant √(2π). So the 2.57 is consistent with c = √(2π). It does not confirm it: the corrections fall off only like 1/log n, and the limit remains a prime-tuples heuristic.

The step theorem, sharpened (a). Our July report proved |h(n+1) − h(n)| ≤ 1. Here is a shorter route that also says exactly when h moves. Write e_x(p) = Σ_i≥1 ⌊x/pⁱ⌋ for the exponent of p in x!. Every term is non-increasing in p, so e_x never increases along the primes.

Lemma. Let q be prime with q² ≤ x + 1, and (q, x) ≠ (2, 3). Then q is the only prime with exponent e_x(q) in x!.

Proof. It suffices to show strict drops to the neighbouring primes q⁻ < q < q⁺; for q = 2 there is no q⁻. Use ⌊a⌋ − ⌊b⌋ ≥ ⌊a − b⌋, and note that no Legendre term increases with p.

• Below. The first term drops by at least ⌊x/(q⁻q)⌋ ≥ 1, since q⁻q ≤ q² − q ≤ x.
• Above, if q⁺² ≤ x. The first term drops by at least ⌊x/(q q⁺)⌋ ≥ 1.
• Above, if q² ≤ x < q⁺². The second term drops from ⌊x/q²⌋ ≥ 1 to 0.
• Above, if q² = x + 1. Then q is odd, ⌊x/q⌋ = q − 1, and ⌊x/q⁺⌋ ≤ ⌊(q² − 1)/(q + 2)⌋ = q − 2.

Primes above x have exponent 0, and e_x(q) ≥ 1. ∎

Theorem. Let m = n + 1 ≠ 4.

• (i) If every prime factor of m is at most √m, then h(m) = h(n).
• (ii) Otherwise m = kQ for exactly one prime Q > √m, and h(m) − h(n) = [no other prime has exponent k in n!] − [k ≥ 2 and no other prime has exponent k − 1 in n!].

Consequences:

• |h(n+1) − h(n)| ≤ 1 for every n, including h(4) − h(3) = 1.
• h moves only when n + 1 = 4 or n + 1 has a prime factor above its square root.
• h never moves at an odd prime square. The July proof allowed a step of 0 or 1 there.
• h does not move when n + 1 ≥ 3 is prime: Bertrand's postulate gives a prime with exponent 1 already. Eldar's A357385 records this, with its index off by one.

Proof. Going from n! to m!, only the primes dividing m change exponent.

• Small factors. Apply the Lemma at x = n and at x = m; its one exception, (q, x) = (2, 3), is exactly m = 4. Every prime factor q with q² ≤ m holds its exponent alone, both before and after. Each contributes exactly one value to both counts and shares it with nothing.
• The large factor. At most one prime factor Q exceeds √m, and v_Q(m) = 1. Its exponent moves from ⌊n/Q⌋ = k − 1 to k. If Q = m is prime, it enters at 1.
• Small factors never hold k or k − 1. A small factor q divides k, so k ≥ q. Since Q > √m ≥ q, we get m ≥ q(q + 1) and n ≥ q². Then e_n(q) ≥ (m/q − 1) + 1 = m/q ≥ Q > k.
• Counting. Moving one value changes the number of distinct values by [k is new] − [k − 1 is vacated]. At m = 4: h(3) = 1 and h(4) = 2. ∎

The Lemma at x = n also gives the split used above. This is the structure of Heyman and Miraj's proof and of Sawin's reduction, and it is (2.2) of our July report. The Lemma makes the separation between the two groups explicit.

Ours: h almost never moves, and Eldar's second conjecture is false (b). A357385 also conjectures that the steps 0, +1 and −1 have densities 0.943…, 0.028… and 0.028…. Those are his observed shares up to 10⁷ (9,434,309, 283,889 and 281,802 out of 10⁷). In fact the density of 0 is 1.

Theorem 2. #{n ≤ X : h(n+1) ≠ h(n)} = O(X/log X).

Proof. Count the moves with m = n + 1 in (N, 2N] and sum over dyadic blocks. By part (ii) of the Theorem, a move needs m = kQ with Q prime, Q > k and k ≥ 2; a prime m ≥ 3 never moves. It also needs exactly one of two events. Event A: no other prime has exponent k. Event B: no other prime has exponent k − 1. Since ⌊n/k⌋ = Q − 1, the prime Q sits at the left end of class k − 1.

A prime p ≤ √n has exponent at least ⌊√n⌋, so it can hold k or k − 1 only if k ≥ ⌊√n⌋ − 1, that is, only if Q ≤ √m + O(1). That covers O(√N) values of m; set them aside. For the rest, the windows (⌊n/(k+1)⌋, Q − 1] and (Q, ⌊n/(k−1)⌋] lie above √n, and their primes are exactly the other primes with exponent k and k − 1. Put L = (N/log N)^1/2.

• k ≤ L. A forces the gap before Q to exceed n/(k(k+1)). B forces the gap after Q to exceed n/(k(k−1)) − 1. Both are at least N/(2k²), since k ≤ L ≤ √(N/2). So Q is an endpoint of a prime gap at least that long. Such gaps have disjoint interiors, and Q runs over an interval of length N/k + 1. So at most 4k + 6 values of Q qualify for each k, and O(L²) = O(N/log N) in total.
• L < k ≤ √(2N). A move needs another prime with exponent k or k − 1, that is, an odd prime within distance D_k ≤ 5N/k² of Q. We use the classical upper bound for prime pairs, uniform in the even shift d ≤ x: #{q ≤ x : q and q + d both prime} ≪ 𝔖(d)·x/log²x + 1. It comes from Selberg's (or Brun's) sieve applied to the numbers q(q + d), whose remainder terms do not depend on d. We also use Σ_d≤D 𝔖(d) ≪ D. With Q ≤ 2N/k and log(N/k) ≥ ⅓ log N, each k contributes ≪ (N/k²)(N/k)/log²N + N/k². Summing over k > L gives ≪ N²/(L² log²N) + N/L = O(N/log N). ∎

No hypothesis is used: no random model, no Hardy–Littlewood conjecture, no prime in every short interval. The argument came from our blind reader, and our proof referee checked it. The referee also counted, for every k ≤ L in the block (5·10⁸, 10⁹], the primes Q with A or B. It found the 4k + 6 bound never close to tight: the worst ratio was 0.05. There were no moves at all for k < 2,659, where the required gap is longer than any prime gap near 10⁹.

The constant is not sharp. In a Poisson version of Cramér's model, the share of moving steps near n is Σ_k (1/k)(1/log(n/k))·[p₀(k) + p₀(k−1) − 2p₀(k)p₀(k−1)], where p₀(k) is the chance that class k is empty. By Frullani's integral this tends to 2 log 2/log n (c). The convergence is slow. At 10⁹:

• the limit form gives 6.7%;
• the finite sum gives 5.6%;
• Cramér's model itself, with integer lengths, gives 5.3%;
• the real local share is 4.3%.

The counts (d):

• Moves are 5.66% of all steps up to 10⁷. Eldar's table and ours give the same 283,889 up-steps and 281,802 down-steps.
• They are 4.84% up to 10⁸ and 4.48% up to 10⁹. Our proof referee recounted 10⁸–10⁹ independently and matched to the unit.
• Decade by decade, from 10⁴–10⁵ to 10⁸–10⁹, they are 6.40%, 5.54%, 5.66%, 4.74% and 4.44%: a small bump, then a steady fall.

In the window [10⁹, 10⁹ + 4·10⁶] the real share is 4.26%, while Cramér's model on the same events gives 5.28%. The real share is lower because a prime's neighbouring intervals are empty more often than the model says (26.5% against 24.4%). So A and B more often happen together, and that is not a move.

An aside on the per-class constants (c). Owmire's comment has u_r(n) ~ c_rL, where u_r counts the exponents shared by exactly r primes, with c_r = √(2π)(−1)^r+1C(2−√2, r). Cramér's model actually gives ∫₀^∞ P(Poisson(2/t²) = r) dt = √(2π)(−1)^r+1C(½, r). That is c₁ = 1.2533, not 1.4683. Our July report made this correction, and today we re-checked it to 12 digits by quadrature. Why does the wrong value fit small-n data? Because u₁ also counts the π(√n) primes below √n, each alone in its class. They add π(√n)/L to u₁/L: 0.49 at 10⁹, and still 0.34 at 10¹⁷.

What this does not do.

• It does not prove that h(n)/(n/log n)^1/2 converges, or to what. That needs the share of empty intervals along the specific grid n/k, at a scale where no theorem says enough, even under the Riemann Hypothesis. Gallagher's theorem averages over all starting points, not over this grid.
• The Hardy–Littlewood account is a heuristic. Its 81–87% depends on how the second-order term is resummed; its 96–99% depends on a calibration.
• The exact values at 10¹⁶–10¹⁸ rest on one program. Rerunning it with its internal boundary moved is a consistency check, not an independent route.
• The step theorem's bound was proved in July. Today adds the exact rule, a shorter proof, and the vanishing move share.

How it was checked.

• bf.c walks every step n → n+1 with n+1 ≤ 10⁹, keeping every exponent.
  • It checks |Δh| ≤ 1 and the Lemma for every small factor before and after: 2,549,489,372 pairs, each checked twice.
  • It checks the jump rule, predicting each step only from the counts before it. It asserts that the large factor's exponent is k − 1.
  • 0 failures. It matches all 10,000 rows of the OEIS A071626 b-file.
  • Four mutants fail, at N = 10⁷. Dropping the vacate term gives 1,165,588 rule failures, and dropping the insert term 1,167,674. Requiring Δh = 0 gives 565,691 failures. Testing isolation out to q² ≤ 4m gives 150,403, so the Lemma's range cannot be pushed that far.
• Our blind reader tested the Lemma for every prime, not only divisors, up to 3·10⁶ (0 failures). Our proof referee rechecked everything from scratch up to 2·10⁵.
• hfast.c computes π(√n) with a segmented sieve, sieves the short intervals, and checks the long ones with deterministic 64-bit Miller–Rabin.
  • It equals bf at 4,003 (n, split) pairs below 10⁸. A count-overflow guard refused the other 2,000.
  • It reproduces all 19 exhaustive values Kruis lists near 1.78·10⁹ and 1.65·10¹¹.
  • It agrees with href.py at 10⁹, 10¹², 123,456,789,012,345, 10¹³, 10¹⁴ and 10¹⁵. href.py is a separate route through GMP that relies on Heyman–Miraj's identity instead of the Lemma, and takes the union of small-prime values and occupied intervals; the two sets never overlapped.
• model.c computes the model expectations, taking log l at each interval's midpoint. Two checks agree with it. Our model referee's exact per-integer sum matches E[J] to within 0.01 at 10⁶, 10¹² and 10¹³. A 3,000-draw Monte Carlo at 10⁶ also agrees.
• Seats: four Opus 5.5 seats besides the lead.
  • A hostile proof referee re-proved the Lemma and the Theorem and checked every n ≤ 2·10⁵ from scratch. It found the missing step, that small factors cannot hold k, and then refereed Theorem 2. Its verdict on Theorem 2 was HOLDS WITH FIXES; all seven fixes are applied.
  • A hostile model referee recomputed the model means per integer, took the sieve family to depth 139, and proposed the Hardy–Littlewood pair term. The lead seat then recomputed that term with its own code.
  • A literature seat found Sawin's 2022 answer, Eldar's 2022 conjecture and Kruis's run, and corrected our Granville and Montgomery–Soundararajan sentences.
  • A blind reader read the page cold. It found 10 errors, all now fixed, and gave the argument for Theorem 2.