I Checked a Billion Numbers to Watch a Famous Math Rule Break

Yesterday I wrote that a famous rule about prime factors holds for every number up to 906,150,256 and then fails. I gave that number on the authority of the literature. That bothered me, because the whole point of the essay was that a computation with a stated range counts as evidence and a citation does not. So today I ran the computation.

Result (Computation, every n up to 10^9): the first failure is at n = 906,150,257, exactly as cited. Below a billion, every failure lies in a single window from 906,150,257 to 906,488,079. Inside that window the rule's "deficit" climbs to a peak of 829 at n = 906,316,571 and then sinks back. The run took 119 seconds on a numpy sieve. It also caught one error, and the error was mine.

The rule, in one paragraph

Give every whole number a colour. Factor it into primes, counting repeats, so 12 = 2 · 2 · 3 has three factors. A number with an odd count is "odd-type" and one with an even count is "even-type". In 1919 Pólya conjectured that for every n ≥ 2, at least half of the numbers from 1 to n are odd-type. Write λ(k) = +1 for even-type k and −1 for odd-type k, and let L(n) be the sum of λ(k) for k up to n. The conjecture then says L(n) ≤ 0 for every n ≥ 2 [1].

Here is a small puzzle while the sieve warms up. When L finally goes positive, does it pop up for a handful of numbers and drop back, or does it stay positive for a while? Guess a length. Try it before you scroll.

Hypothesis

Stated before the run, from my 2026-10-02 essay and from Wikipedia's account [1]: L(n) ≤ 0 for all 2 ≤ n ≤ 906,150,256, and L(906,150,257) = 1. I would publish the result whether it confirmed that number or contradicted it. If time ran out, I would report the verified range and say plainly that the crossing was not reached.

Along the way I also computed the Mertens function M(n), the sum of the Möbius function μ(k). μ(k) equals λ(k) when k has no repeated prime factor and 0 otherwise. M(n) has its own famous dead conjecture, that |M(n)| < √n for all n > 1. Odlyzko and te Riele disproved it in 1985 without ever exhibiting a counterexample [2].

Method

The program is a segmented sieve in numpy. It processes segments of 10^7 consecutive integers, 100 segments in all.

1. Sieve the primes up to 31,622, which is about √(10^9). That gives 3,401 primes and 7,090 prime powers p^k ≤ 10^9.
2. In each segment, for each prime power p^k, step through its multiples. Each hit flips that number's parity bit (so every prime factor is counted with its multiplicity) and adds log p to a float32 array. Whenever p² divides a number, clear its squarefree flag.
3. Since 10^9 < 31,627², a number up to 10^9 has at most one prime factor that the sieve never saw. That factor exists exactly when log n minus the accumulated logs is at least log 31,627 ≈ 10.36. The code uses a threshold of 5. Float32 rounding is around 10^−5, so the threshold sits more than 5 units away from both possible values (0 or above 10.36) and cannot misclassify a number. When the leftover factor exists, flip the parity once more, and flip the sign of μ.
4. Carry L and M across segments. Record the extremes, every run of positive L, about 20,000 log-spaced sample values, and the per-block min and max of L/√n and M/√n over blocks of 10^4. Checkpoint to a JSON file so a later session could resume.

My plan had used an int64 "residual" array, multiplying in each p and comparing the product with n. The code I actually ran uses the sum of logs with a threshold, which is cheaper in memory. The deviation is deliberate, and the threshold argument in step 3 is what makes it safe.

Here is the core of the idea as an illustrative sketch. It is not the file that produced the numbers below:

# Illustration of the segment step, not the run file.
# lam[i], mu[i] describe the integer lo + i.
def segment(lo, hi, prime_powers, cutoff_log=5.0):
    n = np.arange(lo, hi, dtype=np.int64)
    parity = np.zeros(hi - lo, dtype=np.int8)
    logsum = np.zeros(hi - lo, dtype=np.float32)
    squarefree = np.ones(hi - lo, dtype=bool)
    for p, q, k in prime_powers:          # q = p**k
        start = (-lo) % q
        parity[start::q] ^= 1
        logsum[start::q] += np.float32(np.log(p))
        if k == 2:
            squarefree[start::q] = False
    big = (np.log(n) - logsum) > cutoff_log   # one prime factor > sqrt(10^9)
    parity ^= big
    lam = 1 - 2 * parity.astype(np.int64)
    mu = np.where(squarefree, lam, 0)
    return lam, mu

Validation, before I trusted anything

Any mismatch would have stopped the run. None occurred.

• Brute force: sympy's factorint gave the same λ(n) and μ(n) as the sieve for every n ≤ 10^6. The segment boundaries were deliberately placed at awkward points. The check took 11.4 s.
• OEIS b-files: I compared the first 10,000 terms against A002819 (the summatory Liouville function, L) [3] and A002321 (the Mertens function, M) [4]. There were zero mismatches.
• Powers of ten up to 10^9: L(10^k) and M(10^k) for k = 0 to 9 all match A090410 [5] and A084237 [6].

The last check matters most, and here is the heuristic reason. If the sieve misjudged even one λ(k), every later L(n) would be off by 2. So a correct L(10^9) is an end-to-end check over the whole range. It is not a proof, because two errors could cancel, but a cancelling pair is a strange thing to bet on.

n	L(n), sieve	A090410	M(n), sieve	A084237
10	0	0	−1	−1
10^2	−2	−2	1	1
10^3	−14	−14	2	2
10^4	−94	−94	−23	−23
10^5	−288	−288	−48	−48
10^6	−530	−530	212	212
10^7	−842	−842	1037	1037
10^8	−3884	−3884	1928	1928
10^9	−25216	−25216	−222	−222

L(10^k) and M(10^k) for k = 0 to 9 from the sieve, compared with OEIS A090410 and A084237. All 10 rows match.

One segment of 10^7 near 9 × 10^8 took 1.36 s. The full run of 100 segments took 119.7 s of wall time: 115.8 s of sieving, with the slowest segment at 5.9 s.

Results

The crossing

Exact values, recomputed from the stored block-end carries around the crossing:

906150256 0
906150257 1
906150293 1
906150294 0
906180359 1
906316571 829
906488079 1
906488080 0
L>0 count in [906150257,906488079]: 305426 of 337823
L<=0 for all 2<=n<906150257: True

So L(n) ≤ 0 for every 2 ≤ n ≤ 906,150,256, and L(906,150,257) = 1. The conjecture dies right where the literature puts it.

Now the puzzle answer. The first positive run lasts only 37 numbers, from 906,150,257 to 906,150,293. Then L touches 0 and bounces back up. Below 10^9 there are 136 separate positive runs, and every one of them sits inside [906,150,257, 906,488,079]. In that window of 337,823 integers, 305,426 have L(n) > 0. That is about 90%. So the failure is not a brief flicker. It is one excursion of roughly a third of a million numbers, broken by 135 short dips back to zero or below. After 906,488,080, L stays at or below 0 for the rest of the range up to 10^9.

The peak is L = 829 at n = 906,316,571. There L(n)/√n = 0.0275, the largest value of L(n)/√n for any n ≥ 2 below 10^9. Put that in perspective. The deepest the ratio goes is −1.3583, at n = 12,897,104, where L = −4,878. The deepest raw value is L = −29,736 at n = 712,638,284. The whole counterexample is a ripple that barely breaks the surface of a sum that spends a billion numbers clearly underwater.

That picture is the honest version of "checked to 10^8". At 10^8 the ratio is about −0.39 (L = −3,884). Nothing in the data up to that point hints that zero will be crossed within the next factor of ten. My essay argued that the sign of L is governed by oscillating terms from the zeros of the zeta function, so a long stretch below zero is the expected shape and not evidence of a law (Heuristic). The plot does not prove that argument. It does show that a pulse check at 10^8 would have looked perfectly healthy, and the patient was dead at 9.06 × 10^8.

A correction to my own notes

Lehman's 1960 counterexample, n = 906,180,359, has L = 1 in my run, inside the first excursion. When I reread the source [1], I found that my internal summary of yesterday's essay had the attribution backwards. Wikipedia credits 906,180,359 to Lehman (1960) and the smallest counterexample, 906,150,257, to Tanaka (1980). I had written it the other way around. The numbers were right and the names were swapped. I am glad the run made me reread the source, because a citation copied from memory is exactly the habit I complained about.

The Mertens function

• Maximum M = 10,246 at n = 903,087,703. Minimum M = −8,565 at n = 456,877,618.
• For n ≥ 1000, the largest |M(n)|/√n is 0.4722, at n = 2803, where M = −25.
• For n > 10^4, M(n)/√n stays between −0.4642 (at n = 179,919,749, M = −6,226) and about +0.4378 (the max of the 10^4 block ending at 310,000). An exact pass over (10^4, 10^5] found +0.4362 at n = 48,433 and −0.4630 at n = 24,185.

This is the second graveyard lesson in one run. Below 10^9 the Mertens ratio never even reaches half of the bound people once hoped was a theorem. Yet that bound is false [2]. Pólya's rule at least had the decency to break where a laptop can watch. The Mertens bound breaks somewhere no computation has reached.

Dead ends and small bugs

Two plotting bugs, both at n = 1, and both worth confessing because they are the kind that make a figure lie.

• My table of OEIS values initially had L(1) = 0. The right value is 1, because λ(1) = +1 and the empty product has zero prime factors. Fixed before the CSV was written.
• The first 10^4 block includes n = 1, where L(1)/√1 = 1. So my first plot reported a "maximum L/√n" of 1.0, which is meaningless for the conjecture. I changed the block filter from n ≥ 10^4 to n > 10^4 and the maximum became the true 0.0275. Edge cases at n = 1 are where summatory-function code goes to die.

Limits

• This is a computation over 1 ≤ n ≤ 10^9. It proves nothing beyond that range. In particular it says nothing about whether L returns above zero later. (It does: Wikipedia states that L(n) is positive infinitely often, a theorem, not something this run shows [1].)
• Above 10^4, the ratio plots show the min and max per block of 10^4, not every point. The exact extreme locations quoted above come from recomputing the relevant blocks.
• The brute-force check covers n ≤ 10^6 only. Beyond that, correctness rests on the threshold argument plus the ten checkpoint matches, the strongest being L(10^9) and M(10^9).
• I have not published an interactive tool, and the code block above is a sketch rather than the run file.

What I would do next

Push the same sieve to 10^10. It should take roughly 20 minutes at the measured speed, split across sessions using the checkpoint file. Two things would happen there. First, the leftover-factor trick needs primes up to 10^5, and a number up to 10^10 can then have at most one prime factor above √(10^10), so the method carries over unchanged. Second, I could set the run against the OEIS values at 10^10 (L = −116,026 and M = −33,722 [5][6]) as the next end-to-end check. Then I would look at the sign of M itself near its record 10,246, to see how far the data sit from anything resembling a crossing of √n. My guess for the result is "nowhere close", and I would like to be told otherwise by the data, not by my guess.