mertens_meter.py

1import numpy as np234KNOWN = [5    14.1347,6    21.0220,7    25.0109,8    30.4249,9    32.9351,10    37.5862,11    40.9187,12    43.3271,13]1415N_MERTENS = 5000016N_STACK = 20017N_READOUT = 2000018N_SAMPLES = 819219BAND = (8.0, 55.0)20PEAK_FACTOR = 3.0212223def mobius(n):24    mu = np.ones(n + 1, dtype=np.int64)25    composite = np.zeros(n + 1, dtype=bool)26    for p in range(2, n + 1):27        if composite[p]:28            continue29        composite[p::p] = True30        mu[p::p] *= -131        square = p * p32        if square <= n:33            mu[square::square] = 034    mu[0] = 035    return mu363738def mertens(mu):39    running = np.cumsum(mu[1:].astype(np.int64))40    return np.concatenate((np.zeros(1, dtype=np.int64), running))414243def stack(mu, n):44    nodes = 045    brightness = np.empty(n, dtype=np.int64)46    for b in range(1, n + 1):47        brightness[b - 1] = mu[b::b][: n // b].sum()48        nodes += int((np.gcd(np.arange(b + 1), b) == 1).sum())49    return nodes, brightness505152def readout(m, n):53    breaches = 054    for x in range(1, n + 1):55        if int(m[x // np.arange(1, x + 1)].sum()) != 1:56            breaches += 157    return breaches585960def normalised(m, n):61    x = np.arange(1, n + 1, dtype=np.float64)62    return x, m[1 : n + 1] / np.sqrt(x)636465def spectrum(x, values, samples):66    log_x = np.log(x)67    grid = np.linspace(log_x[0], log_x[-1], samples)68    resampled = np.interp(grid, log_x, values)69    step = (grid[-1] - grid[0]) / (samples - 1)70    amplitude = np.fft.rfft(resampled * np.hanning(samples))71    gamma = np.fft.rfftfreq(samples, d=step) * 2.0 * np.pi72    return gamma, np.abs(amplitude) ** 2, step737475def peaks(gamma, power, band, factor):76    inside = (gamma > band[0]) & (gamma < band[1])77    threshold = np.median(power[inside]) * factor78    found = []79    for i in range(1, power.size - 1):80        if not inside[i]:81            continue82        if power[i] > power[i - 1] and power[i] > power[i + 1] and power[i] > threshold:83            found.append(i)84    return found, threshold858687def main():88    mu = mobius(N_MERTENS)89    m = mertens(mu)9091    nodes, brightness = stack(mu[: N_STACK + 1], N_STACK)92    depth = m[N_STACK // np.arange(1, N_STACK + 1)]93    agreements = int((brightness == depth).sum())94    breaches = readout(m, N_READOUT)9596    x, values = normalised(m, N_MERTENS)97    gamma, power, step = spectrum(x, values, N_SAMPLES)98    found, threshold = peaks(gamma, power, BAND, PEAK_FACTOR)99100    print(f"domain  N_mertens = {N_MERTENS}  N_stack = {N_STACK}")101    print(f"samples {N_SAMPLES}  log step {step:.8f}  bin width {gamma[1]:.6f}")102    print(f"stack   nodes {nodes}  brightness min {brightness.min()} max {brightness.max()}")103    print(f"stack   sum mu(kb) equals M(floor(N/b)) at {agreements} of {N_STACK} denominators")104    print(f"readout sum M(floor(x/n)) = 1 through x = {N_READOUT}, breaches {breaches}")105    print(f"band    {BAND[0]:g} to {BAND[1]:g}  peaks {len(found)}  threshold {threshold:.4g}")106    print(f"mertens M({N_MERTENS}) = {int(m[-1])}")107    print()108    print("known      detected  error")109    for target in KNOWN:110        best = min(found, key=lambda i: abs(gamma[i] - target))111        print(f"{target:8.4f}   {gamma[best]:7.2f}   {abs(gamma[best] - target):5.2f}")112113114if __name__ == "__main__":115    main()