"""Independent check of the m = 5 tuple: prime counting by the Lucy_Hedgehog (Legendre-type) method in O(n^{3/4})
operations, vectorised with numpy, plus a deterministic Miller-Rabin test (bases valid for n < 3.3e24).
Shares no code with zhang.py (segmented sieve).  Checks pi(k), the first prime > k, that 5274442853 is the
k-th prime > k, and the diameter."""
import numpy as np, math, time, sys
def primepi(n):
    r = math.isqrt(n)
    nsmall = n // r                      # small values v = 0..n//r
    Ssmall = np.arange(nsmall + 1, dtype=np.int64) - 1; Ssmall[0] = 0
    i = np.arange(1, r + 1, dtype=np.int64)
    Slarge = np.concatenate([[0], n // i - 1]).astype(np.int64)   # Slarge[i] = S(n//i)
    for p in range(2, r + 1):
        if Ssmall[p] == Ssmall[p - 1]:
            continue
        sp = Ssmall[p - 1]; p2 = p * p
        imax = min(r, n // p2)
        if imax >= 1:
            ii = np.arange(1, imax + 1, dtype=np.int64); d = ii * p
            vals = np.where(d <= r, Slarge[np.minimum(d, r)], Ssmall[np.minimum(n // d, nsmall)])
            Slarge[1:imax + 1] -= vals - sp
        if p2 <= nsmall:
            v = np.arange(p2, nsmall + 1, dtype=np.int64)
            Ssmall[p2:] -= Ssmall[v // p] - sp
    return int(Slarge[1])
def is_prime(n):
    if n < 2: return False
    for q in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):
        if n % q == 0: return n == q
    d, s = n - 1, 0
    while d % 2 == 0: d //= 2; s += 1
    for a in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41):
        x = pow(a, d, n)
        if x in (1, n - 1): continue
        for _ in range(s - 1):
            x = x * x % n
            if x == n - 1: break
        else: return False
    return True
if __name__ == "__main__":
    # self-test against known values
    for n, v in [(10**6, 78498), (10**7, 664579), (10**8, 5761455), (10**9, 50847534)]:
        t = time.time(); got = primepi(n); print("pi(%d) = %d (expected %d) %s %.1fs" % (n, got, v, got == v, time.time() - t), flush=True)
    k = 234360000
    t = time.time(); pk = primepi(k); print("pi(k) =", pk, "%.1fs" % (time.time() - t), flush=True)
    print("primes in (k, 234360017]:", [x for x in range(k + 1, 234360018) if is_prime(x)])
    last = 5274442853
    t = time.time(); pl = primepi(last); print("pi(5274442853) =", pl, "%.1fs" % (time.time() - t), flush=True)
    print("is_prime(last) =", is_prime(last), " pi(last) - pi(k) =", pl - pk, " (should be k =", k, ")")
    print("diameter =", last - 234360017)