abc quality record

"""Verifier for "abc quality record" (statement in abc_quality_record.json).

Background: an abc triple is a + b = c with a and b positive and sharing no
common factor. Its radical rad(abc) is the product of the distinct primes
dividing a*b*c, and its quality is

    q = log(c) / log(rad(abc))

Most triples have q < 1. The abc conjecture says that q > 1 + epsilon happens
only finitely often. The highest quality known is Reyssat's
2 + 3**10 * 109 = 23**5, with q = 1.6299...

Problem: find a triple with q > 1.63, beating the record.

The comparison is made without logarithms or floating point. For positive
integers, q > Q_NUM / Q_DEN is equivalent to

    c**Q_DEN > rad**Q_NUM

which is evaluated with exact integers (c**100 > rad**163 by default).

To compute rad the verifier needs the prime factorisations of a, b and c.
It cannot find them itself, so the solver supplies them and the verifier
checks them: every listed factor must be a proven prime, and the product must
be exactly the number.

Expected submission: a JSON object with the three numbers and their
factorisations as [prime, exponent] pairs (the empty list stands for 1):

    {"a": "1", "b": "8", "c": "9",
     "factors_a": [], "factors_b": [["2", "3"]], "factors_c": [["3", "2"]]}

What the verifier checks, in order:
  1. The payload is at most MAX_BYTES characters and is a JSON object with
     exactly the six keys above.
  2. a, b, c are well-formed decimal strings of at most MAX_DIGITS digits
     (see _parse_int).
  3. a >= 1, b >= 1, a + b == c, and gcd(a, b) == 1 (then a, b and c are
     automatically pairwise coprime).
  4. Each factorisation is valid (see _parse_factorisation) and multiplies
     out to exactly its number.
  5. rad = product of all distinct primes appearing in the three lists.
  6. Size guard: c**Q_DEN and rad**Q_NUM must each have at most
     MAX_POWER_BITS bits (checked before computing them).
  7. c**Q_DEN > rad**Q_NUM.
"""
import json
import math

from sympy import isprime

# The quality threshold is Q_NUM / Q_DEN = 1.63, just above Reyssat's 1.6299.
# The tests set Q_NUM = 162 (threshold 1.62) to check that Reyssat's triple
# is then accepted.
Q_NUM = 163
Q_DEN = 100
# Maximum number of digits of a, b and c.
MAX_DIGITS = 200
# Maximum number of digits of each prime in a factorisation. Below about
# 3.3 * 10**24, sympy's isprime gives a proven answer rather than a probable one.
MAX_PRIME_DIGITS = 24
# Maximum size, in bits, of c**Q_DEN and of rad**Q_NUM. It only protects
# against modified constants: at the defaults these are at most about
# 330,000 bits, so this limit is never reached.
MAX_POWER_BITS = 1_000_000
# Maximum size of the whole submission, in characters.
MAX_BYTES = 20_000

# The keys the submission must contain.
KEYS = ("a", "b", "c", "factors_a", "factors_b", "factors_c")


# ---------------------------------------------------------------------------
# Input-parsing helpers. The same code is copied into every verifier that
# needs it, so that each file works on its own.
# ---------------------------------------------------------------------------

def _parse_int(s, max_digits, allow_neg=False):
    """Convert a strictly formatted decimal string to an int, or raise ValueError.

    Accepted: "0", "42", and "-42" when allow_neg is True.
    Rejected: anything that is not a string (for example a JSON number), "",
    "+42", " 42", "007", "-0", non-ASCII digits such as "٣", and strings with
    more than max_digits digits.
    """
    # Cheap length check before any real work: at most max_digits digits,
    # plus one character for an optional minus sign.
    if not isinstance(s, str) or len(s) > max_digits + 1:
        raise ValueError
    # Split off the sign. A "-" is only allowed where negative values make sense.
    neg = s.startswith("-")
    if neg and not allow_neg:
        raise ValueError
    body = s[1:] if neg else s
    # What remains must be 1 to max_digits ASCII digits with no leading zero
    # ("0" on its own is fine). isascii() is needed because isdigit() also
    # accepts digits from other scripts.
    if not (body.isascii() and body.isdigit()) or len(body) > max_digits or (len(body) > 1 and body[0] == "0"):
        raise ValueError
    # "-0" is not a canonical way to write zero.
    if neg and body == "0":
        raise ValueError
    # Python (3.11 and later) refuses int() on strings longer than 4300 digits
    # by default, so the number is built from chunks of 4000 digits.
    n = 0
    for i in range(0, len(body), 4000):
        chunk = body[i:i + 4000]
        n = n * 10 ** len(chunk) + int(chunk)
    return -n if neg else n


def _load(solution, keys, max_bytes):
    """Parse the submission as a JSON object with exactly the given keys, or raise."""
    # Refuse oversized payloads before handing them to the JSON parser.
    if not isinstance(solution, str) or len(solution) > max_bytes:
        raise ValueError
    d = json.loads(solution)
    # It must be an object (not a list, string or number) whose keys are
    # exactly the expected ones: nothing missing, nothing extra.
    if not isinstance(d, dict) or set(d) != set(keys):
        raise ValueError
    return d


def _parse_factorisation(factors, max_prime_digits, max_bits):
    """Validate a claimed prime factorisation and return it as a list of (p, e) ints.

    `factors` is a JSON list of [p, e] pairs of decimal strings, meaning
    n = p1**e1 * p2**e2 * ... (the empty list stands for n = 1). Each p must
    be a prime of at most max_prime_digits digits and may appear only once;
    each exponent e must be between 1 and 999999. Raises ValueError otherwise.

    max_bits guards against absurd inputs such as [["2", "999999"]], whose
    product would have a million bits. Every prime p is at least
    2**(bit length of p - 1), so the product is at least 2**bits, where bits
    is the running total of e * (bit length of p - 1). Once that total exceeds
    max_bits (the size of the number the factorisation is meant to reproduce),
    the factorisation cannot be right and is rejected. Callers can therefore
    multiply the returned pairs out without risk.
    """
    if not isinstance(factors, list):
        raise ValueError
    pairs, seen, bits = [], set(), 0
    for item in factors:
        # Each item must be a two-element list [p, e].
        if not (isinstance(item, list) and len(item) == 2):
            raise ValueError
        p = _parse_int(item[0], max_prime_digits)
        e = _parse_int(item[1], 6)
        # p must be at least 2 and not repeated; e must be at least 1.
        if p < 2 or p in seen or e < 1:
            raise ValueError
        # Size guard first (cheap), then the proven primality test.
        bits += e * (p.bit_length() - 1)
        if bits > max_bits or not isprime(p):
            raise ValueError
        seen.add(p)
        pairs.append((p, e))
    return pairs


# ---------------------------------------------------------------------------
# The actual check.
# ---------------------------------------------------------------------------

def _check(solution):
    # Steps 1-2: read the JSON object and parse a, b, c.
    sol = _load(solution, KEYS, MAX_BYTES)
    a, b, c = (_parse_int(sol[k], MAX_DIGITS) for k in ("a", "b", "c"))

    # Step 3: a genuine abc triple: positive, a + b = c, no common factor.
    if min(a, b) < 1 or a + b != c or math.gcd(a, b) != 1:
        return False

    # Steps 4-5: each factorisation must multiply out to its number; collect
    # the distinct primes of all three.
    primes = set()
    for v, key in ((a, "factors_a"), (b, "factors_b"), (c, "factors_c")):
        pairs = _parse_factorisation(sol[key], MAX_PRIME_DIGITS, v.bit_length())
        if math.prod(p ** e for p, e in pairs) != v:
            return False
        primes.update(p for p, _ in pairs)
    rad = math.prod(primes)  # the radical of a*b*c

    # Step 6: size guard, before any large power is computed.
    if Q_DEN * c.bit_length() > MAX_POWER_BITS or Q_NUM * rad.bit_length() > MAX_POWER_BITS:
        return False

    # Step 7: quality above Q_NUM / Q_DEN, i.e. c**Q_DEN > rad**Q_NUM, exactly.
    return c ** Q_DEN > rad ** Q_NUM


def verify(solution: str, rng) -> bool:
    """Platform entry point: return True only for a valid, correct solution.

    `rng` (a numpy.random.Generator) is part of the platform interface but is
    not used: this check is deterministic. Any exception (malformed JSON, a bad
    number, ...) means the submission is rejected, so this function never raises.
    """
    try:
        return _check(solution) is True
    except Exception:
        return False