#!/usr/bin/env python3
from fractions import Fraction

R = 200
C = Fraction(4029639598, 25970038185)
ACTIVE = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 43]
ALLOCATION = [
    (1, 2, 4029639598, 25970038185),
    (1, 3, 597400199, 51940076370),
    (2, 3, 3432239399, 51940076370),
    (2, 5, 30432359, 51940076370),
    (3, 5, 177892183, 11130016365),
    (3, 7, 6149035241, 311640458220),
    (4, 5, 1, 45),
    (5, 7, 5536939, 903305676),
    (5, 11, 29881539, 3312120812),
    (6, 11, 37106416381, 5713408400700),
    (6, 13, 25678291319, 5713408400700),
    (7, 13, 1, 120),
    (8, 17, 1, 153),
    (9, 13, 3520187831, 34280450404200),
    (9, 17, 1970727917, 623280916440),
    (9, 19, 34256599957, 17140225202100),
    (10, 19, 1, 231),
    (11, 19, 117890228609, 51420675606300),
    (11, 23, 34208283533, 25710337803150),
    (12, 23, 1, 325),
    (13, 23, 1, 378),
    (14, 29, 1, 435),
    (15, 29, 1, 496),
    (16, 29, 213831204731, 174324746984752),
    (16, 31, 3197967125333, 5752716650496816),
    (17, 31, 1, 630),
    (18, 31, 142776482538959, 143817916262420400),
    (18, 37, 2286634159870117, 5321262901709554800),
    (19, 31, 1, 780),
    (20, 41, 1, 861),
    (21, 43, 1, 946),
    (22, 41, 1, 1035),
    (23, 43, 1, 1128),
    (24, 37, 10658441105846658564953, 17935134484658062065193800),
    (24, 43, 796497000815658011939, 3587026896931612413038760),
    (25, 31, 1, 1326),
    (26, 43, 1, 1431),
    (27, 47, 1, 1540),
    (28, 53, 1, 1653),
    (29, 47, 1, 1770),
    (30, 43, 1, 1891),
    (31, 61, 1, 2016),
    (32, 41, 1, 2145),
    (33, 41, 1, 2278),
    (34, 61, 1, 2415),
    (35, 47, 1, 2556),
    (36, 47, 1, 2701),
    (37, 41, 1, 2850),
    (38, 73, 1, 3003),
    (39, 79, 1, 3160),
    (40, 43, 1, 3321),
    (41, 47, 1, 3486),
    (42, 61, 1, 3655),
    (43, 47, 1, 3828),
    (44, 71, 1, 4005),
    (45, 83, 1, 4186),
    (46, 73, 1, 4371),
    (47, 67, 1, 4560),
    (48, 61, 1, 4753),
    (49, 53, 1, 4950),
    (50, 53, 1, 5151),
    (51, 67, 1, 5356),
    (52, 59, 1, 5565),
    (53, 103, 1, 5778),
    (54, 107, 1, 5995),
    (55, 103, 1, 6216),
    (56, 109, 1, 6441),
    (57, 113, 1, 6670),
    (58, 79, 1, 6903),
    (59, 103, 1, 7140),
    (60, 103, 1, 7381),
    (61, 107, 1, 7626),
    (62, 83, 1, 7875),
    (63, 67, 1, 8128),
    (64, 83, 1, 8385),
    (65, 67, 1, 8646),
    (66, 113, 1, 8911),
    (67, 83, 1, 9180),
    (68, 79, 1, 9453),
    (69, 71, 1, 9730),
    (70, 73, 1, 10011),
    (71, 103, 1, 10296),
    (72, 89, 1, 10585),
    (73, 103, 1, 10878),
    (74, 101, 1, 11175),
    (75, 79, 1, 11476),
    (76, 101, 1, 11781),
    (77, 107, 1, 12090),
    (78, 97, 1, 12403),
    (79, 83, 1, 12720),
    (80, 101, 1, 13041),
    (81, 101, 1, 13366),
    (82, 103, 1, 13695),
    (83, 89, 1, 14028),
    (84, 103, 1, 14365),
    (85, 101, 1, 14706),
    (86, 89, 1, 15051),
    (87, 89, 1, 15400),
    (88, 101, 1, 15753),
    (89, 97, 1, 16110),
    (90, 97, 1, 16471),
    (91, 97, 1, 16836),
    (92, 97, 1, 17205),
    (93, 97, 1, 17578),
    (94, 97, 1, 17955),
    (95, 97, 1, 18336),
    (96, 97, 1, 18721),
    (97, 191, 1, 19110),
    (98, 197, 1, 19503),
    (99, 163, 1, 19900),
    (100, 191, 1, 20301),
    (101, 109, 1, 20706),
    (102, 131, 1, 21115),
    (103, 199, 1, 21528),
    (104, 179, 1, 21945),
    (105, 173, 1, 22366),
    (106, 199, 1, 22791),
    (107, 127, 1, 23220),
    (108, 199, 1, 23653),
    (109, 151, 1, 24090),
    (110, 211, 1, 24531),
    (111, 193, 1, 24976),
    (112, 193, 1, 25425),
    (113, 197, 1, 25878),
    (114, 181, 1, 26335),
    (115, 167, 1, 26796),
    (116, 197, 1, 27261),
    (117, 173, 1, 27730),
    (118, 193, 1, 28203),
    (119, 179, 1, 28680),
    (120, 173, 1, 29161),
    (121, 163, 1, 29646),
    (122, 137, 1, 30135),
    (123, 137, 1, 30628),
    (124, 167, 1, 31125),
    (125, 197, 1, 31626),
    (126, 197, 1, 32131),
    (127, 157, 1, 32640),
    (128, 157, 1, 33153),
    (129, 149, 1, 33670),
    (130, 131, 1, 34191),
    (131, 163, 1, 34716),
    (132, 199, 1, 35245),
    (133, 211, 1, 35778),
    (134, 151, 1, 36315),
    (135, 151, 1, 36856),
    (136, 149, 1, 37401),
    (137, 193, 1, 37950),
    (138, 199, 1, 38503),
    (139, 197, 1, 39060),
    (140, 197, 1, 39621),
    (141, 173, 1, 40186),
    (142, 149, 1, 40755),
    (143, 173, 1, 41328),
    (144, 199, 1, 41905),
    (145, 167, 1, 42486),
    (146, 163, 1, 43071),
    (147, 163, 1, 43660),
    (148, 173, 1, 44253),
    (149, 173, 1, 44850),
    (150, 157, 1, 45451),
    (151, 163, 1, 46056),
    (152, 191, 1, 46665),
    (153, 199, 1, 47278),
    (154, 193, 1, 47895),
    (155, 163, 1, 48516),
    (156, 173, 1, 49141),
    (157, 197, 1, 49770),
    (158, 179, 1, 50403),
    (159, 167, 1, 51040),
    (160, 181, 1, 51681),
    (161, 163, 1, 52326),
    (162, 199, 1, 52975),
    (163, 193, 1, 53628),
    (164, 197, 1, 54285),
    (165, 167, 1, 54946),
    (166, 181, 1, 55611),
    (167, 191, 1, 56280),
    (168, 199, 1, 56953),
    (169, 173, 1, 57630),
    (170, 193, 1, 58311),
    (171, 181, 1, 58996),
    (172, 193, 1, 59685),
    (173, 191, 1, 60378),
    (174, 197, 1, 61075),
    (175, 179, 1, 61776),
    (176, 191, 1, 62481),
    (177, 197, 1, 63190),
    (178, 191, 1, 63903),
    (179, 199, 1, 64620),
    (180, 191, 1, 65341),
    (181, 197, 1, 66066),
    (182, 199, 1, 66795),
    (183, 197, 1, 67528),
    (184, 197, 1, 68265),
    (185, 199, 1, 69006),
    (186, 193, 1, 69751),
    (187, 193, 1, 70500),
    (188, 191, 1, 71253),
    (189, 199, 1, 72010),
    (190, 197, 1, 72771),
    (191, 197, 1, 73536),
    (192, 197, 1, 74305),
    (193, 199, 1, 75078),
    (194, 197, 1, 75855),
    (195, 197, 1, 76636),
    (196, 199, 1, 77421),
    (197, 199, 1, 78210),
    (198, 199, 1, 79003),
    (199, 347, 1, 79800),
    (200, 241, 1, 80601),
]

def factor(a):
    out = {}
    p = 2
    while p * p <= a:
        while a % p == 0:
            out[p] = out.get(p, 0) + 1
            a //= p
        p += 1
    if a > 1:
        out[a] = out.get(a, 0) + 1
    return out

def primes_through(n):
    return [a for a in range(2, n + 1) if factor(a) == {a: 1}]

def previous_prime(p):
    q = p - 1
    while factor(q) != {q: 1}:
        q -= 1
    return q

def require(condition, message):
    """Fail unconditionally, including when Python is run with ``-O``."""
    if not condition:
        raise RuntimeError(message)

x = {}
for r, q, num, den in ALLOCATION:
    require((r, q) not in x and 1 <= r <= R and r + 1 <= q <= 2*r + 1,
            f"invalid or duplicate allocation entry ({r}, {q})")
    x[r, q] = Fraction(num, den)
    require(x[r, q] > 0, f"nonpositive allocation entry ({r}, {q})")

for r in range(1, R + 1):
    row_sum = sum((v for (s, q), v in x.items() if s == r), Fraction())
    require(row_sum == Fraction(1, (r + 1)*(2*r + 1)),
            f"row identity failed at r={r}")

loads = {p: Fraction() for p in primes_through(2*R + 1)}
for (r, q), value in x.items():
    for p, exponent in factor(q).items():
        loads[p] += exponent*value
for p, load in loads.items():
    require(load <= C/Fraction(p - 1),
            f"finite capacity exceeded at p={p}")
require([p for p, load in loads.items() if load == C/Fraction(p - 1)]
        == ACTIVE, "active-prime list does not match the saturated rows")

# For r>200 route the whole row to the least prime q>r.  Only primes
# 211 <= p <= 401 can overlap the finite certificate; check them exactly.
for p in [q for q in primes_through(401) if q > R]:
    p0 = previous_prime(p)
    tail = sum((Fraction(1, (r + 1)*(2*r + 1))
                for r in range(max(R + 1, p0), p)), Fraction())
    require(loads.get(p, Fraction()) + tail <= C/Fraction(p - 1),
            f"finite/tail overlap capacity exceeded at p={p}")

# For p>401, Nagura gives p < 6*p0/5, and hence
# (p-1)*sum_{r=p0}^{p-1} alpha_r < 3/25 < C.
require(Fraction(3, 25) < C, "Nagura-tail numerical inequality failed")
print("PASS: exact infinite allocation certificate")