#!/usr/bin/env python3
"""Exact arithmetic checks for the fixed-x convergence constants.

This file does not replace the symbolic denominator-contiguity certificate.
It verifies the numerical inequalities used after that exact identity is
known:

  c_N/c_(N-1) >= 29*N^2,
  sum k^3*(1/3)^k < 5,
  beta(x_official) < 4*10^-19 < 1.
"""

from fractions import Fraction as F


R = 151931373056001
X = F(1, R)
TRANSFER_BOUND = 400_000_000


# Clearing the positive denominator (6N+1)(6N+5), the difference between
#
# 576*N^2*(2N+1)^2 / ((6N+1)(6N+5))
#
# and 29*N^2 has numerator
#
# N^2*(431 + 1260*N + 1260*N^2).
assert all(coefficient > 0 for coefficient in (431, 1260, 1260))


# Exact closed form for sum_{k>=1} k^3 t^k at t=1/3.
theta3_sum = F(1, 3) * (1 + F(4, 3) + F(1, 9)) / (1 - F(1, 3)) ** 4
assert theta3_sum == F(33, 8)
assert theta3_sum < 5


# Entrywise constants in
#
#     |(M_n)_(i,j)| <= K_(i,j) * u^(i+2-j),  u >= 3,
#
# on |x|=1/4.  Each left side below is the exact sum-of-absolute-
# coefficients estimate described in the report.
raw_estimates = [
    [
        F(4 * (144 + 288 + 144) + (99 + 333 + 229 + 114 + 40 + 64), 8),
        F(4 * (432 + 864 + 432) + (243 + 909 + 868 + 80 + 272), 8),
        F(4 * (432 + 864 + 432) + (153 + 648 + 860 + 360), 8),
        F(4 * 144, 8),
    ],
    [F(1), F(3), F(3), F(1)],
    [
        F(1) + F(9 + 63 + 158 + 168 + 64, 4 * 144),
        F(1) + F(36 + 189 + 316 + 168, 4 * 72),
        F(1) + F(54 + 189 + 158, 4 * 36),
        F(1),
    ],
    [
        F(1) + F(54 + 378 + 948 + 1008 + 384, 4 * 288)
        + F(18 + 45 + 251 + 1086 + 1384 + 576, 16 * 288),
        F(1) + F(153 + 657 + 1292 + 2064 + 1072, 4 * 144)
        + F(72 + 702 + 1069 + 2508 + 1512, 16 * 144),
        F(1) + F(180 + 891 + 1450 + 1116, 4 * 72)
        + F(108 + 864 + 1385 + 1422, 16 * 72),
        F(1) + F(6 + 33 + 58 + F(14, 9), 16)
        + F(4 + 32 + 63, 64),
    ],
]
raw_caps = [
    [400, 1200, 1200, 72],
    [1, 3, 3, 1],
    [2, 4, 4, 1],
    [5, 13, 17, 9],
]
for row, caps in zip(raw_estimates, raw_caps):
    for estimate, cap in zip(row, caps):
        assert estimate <= cap
        assert cap < 10_000


beta = (
    F(TRANSFER_BOUND, 29)
    * X**2
    / (F(1, 4) * (1 - X))
)
assert beta == F(3125, 1307443596565949700399927)
assert beta < F(4, 10**19)
assert beta < 1


print("PASS: exact fixed-x convergence constants")
print("sum k^3/3^k =", theta3_sum)
print("entrywise transfer constants < 10000")
print("beta =", beta)
print("beta < 4e-19 < 1")