"""[Ancillary file rsc-credits.py.]  Rigorous rational check of the two inequalities that are evaluated in floating
point by rsc-verify.py: the credit inequalities (Lemma "Credits for RSC") and the light-pair classification
(Lemma "Light pairs for RSC") of the RSC section, for p2 = 7/100, p3 = 77/1000, y = 31/100,
omega2 = 272657/250000, omega3 = 288677/250000, d3min = (3, 2, 2, 1, ...), d2min = (4, 3, 2, 1, ...).
All logarithms are replaced by truncated series with the right direction of the error:
  ln(1+x) <= sum_{k=1}^{2K+1} (-1)^(k+1) x^k / k   (alternating series, cut after a positive term; 0 < x < 1),
  ln 2    >= 2 sum_{k=0}^{K-1} 3^-(2k+1) / (2k+1)   (artanh(1/3) series, every truncation is a lower bound),
  -ln(1-p) <= sum_{k=1}^{K-1} p^k / k + p^K / (K (1-p))   (tail bounded by a geometric series),
  ln((1+z)/(1-z)) >= 2 (z + z^3/3)                   (two terms of the artanh series).
Usage: python3 rsc-credits.py"""
from fractions import Fraction as Fr

p2, p3, y = Fr(7, 100), Fr(77, 1000), Fr(31, 100)
q2, q3 = 1 - p2, 1 - p3
om2, om3 = Fr(272657, 250000), Fr(288677, 250000)

def ln1p_upper(x, K):      # 0 < x < 1
    return sum(Fr((-1) ** (k + 1)) * x ** k / k for k in range(1, 2 * K + 2))
def ln2_lower(K):
    return 2 * sum(Fr(1, 3) ** (2 * k + 1) / (2 * k + 1) for k in range(K))
def mlog1m_upper(p, K):    # -ln(1-p) <= ...
    return sum(p ** k / k for k in range(1, K)) + p ** K / (K * (1 - p))

ok = True
# (W) credits: ln om2 <= (1 - rho1) ln 2,  ln om3 <= (1 - rho1^2) (ln 2 - eps3),  1 <= om <= 2q
rho1 = 1 / (1 + y / 2 * q3)
assert rho1 == Fr(200000, 228613)
K = 7
L2 = ln2_lower(K)
lhs2, rhs2 = ln1p_upper(om2 - 1, K), (1 - rho1) * L2
lhs3, rhs3 = ln1p_upper(om3 - 1, K), (1 - rho1 ** 2) * (L2 - mlog1m_upper(p3, K))
for name, l, r in (("omega2", lhs2, rhs2), ("omega3", lhs3, rhs3)):
    print(f"credit {name}: ln omega <= {float(l):.12f}  <=  bound {float(r):.12f}  margin {float(r - l):.3e}  {'OK' if l <= r else 'FAIL'}")
    ok &= l <= r
ok &= 1 <= om2 <= 2 * q2 and 1 <= om3 <= 2 * q3
print("1 <= omega2 <= 2 q2:", 1 <= om2 <= 2 * q2, " 1 <= omega3 <= 2 q3:", 1 <= om3 <= 2 * q3)

# (L) light pairs: 3-light iff 2^-m b^-d >= (1 + p3/y) eps3 y / delta3;  2-light iff 2^-m b^-d >= eps2 y / ln 2
b = 2 + 2 * max(p2, p3 - p2) / y
eps3, eps2 = mlog1m_upper(p3, 3), mlog1m_upper(p2, 3)
z = (2 * q3 - 1) / (2 * q3 + 1)            # 2 q3 = (1+z)/(1-z)
delta3 = 2 * (z + z ** 3 / 3)
R3 = (1 + p3 / y) * eps3 * y / delta3
R2 = eps2 * y / ln2_lower(4)
print(f"b = {b}, upper bounds: R3 <= {float(R3):.6f}, R2 <= {float(R2):.6f}")
dmin3 = {1: 3, 2: 2, 3: 2}; dmin2 = {1: 4, 2: 3, 3: 2}
for m in range(1, 37):
    for d in range(1, dmin3.get(m, 1)):
        L = Fr(1, 2 ** m) / b ** d; print(f"3-light ({m},{d}): {float(L):.5f} >= {float(R3):.5f}", L >= R3); ok &= L >= R3
    for d in range(1, dmin2.get(m, 1)):
        L = Fr(1, 2 ** m) / b ** d; print(f"2-light ({m},{d}): {float(L):.5f} >= {float(R2):.5f}", L >= R2); ok &= L >= R2
print("ALL OK" if ok else "FAILED")