synthesize_two_block_clocks.py

#!/usr/bin/env python3
"""Exact synthesis from a supplied semilinear (n,m,h) relation.

This implements the arithmetic stage, not PDA-to-semilinear extraction.
Fourier-Motzkin uses Fractions and returns checkable infeasibility multipliers.
No external dependencies. Run this file to reproduce the accompanying JSON.
"""
from fractions import Fraction as F
from functools import reduce
from math import gcd, lcm
import json


def linear(base, *periods):
    return {"base": tuple(base), "periods": [tuple(v) for v in periods]}


def fibers(components, n):
    result = []
    for item in components:
        base = item["base"]
        if base[0] > n:
            continue
        positive = [v for v in item["periods"] if v[0] > 0]
        zero = [v for v in item["periods"] if v[0] == 0]

        def visit(i, current):
            if i == len(positive):
                if current[0] == n:
                    result.append(linear(current, *zero))
                return
            v = positive[i]
            for k in range((n - current[0]) // v[0] + 1):
                visit(i + 1, tuple(current[j] + k * v[j] for j in range(3)))
        visit(0, base)
    return result


def solve_exact(rows, dimensions):
    """rows are (coefficients, rhs, provenance): coefficients*x <= rhs."""
    def eliminate(current, d):
        # A constant contradiction carries nonnegative original-row weights.
        for a, b, proof in current:
            if not any(a) and b < 0:
                return None, proof
        if d == 0:
            return [], None
        positive = [r for r in current if r[0][-1] > 0]
        negative = [r for r in current if r[0][-1] < 0]
        reduced = [(a[:-1], b, proof) for a, b, proof in current if a[-1] == 0]
        for ap, bp, pp in positive:
            for an, bn, pn in negative:
                u, v = 1 / ap[-1], -1 / an[-1]
                proof = {i: u * pp.get(i, F(0)) + v * pn.get(i, F(0))
                         for i in pp.keys() | pn.keys()}
                reduced.append((tuple(u * x + v * y for x, y in zip(ap[:-1], an[:-1])),
                                u * bp + v * bn, proof))
        # Deduplicate exact rows to limit redundant elimination work.
        unique = {}
        for a, b, proof in reduced:
            unique.setdefault((a, b), (a, b, proof))
        prefix, failure = eliminate(list(unique.values()), d - 1)
        if failure is not None:
            return None, failure
        lower, upper = [], []
        for a, b, _ in current:
            remainder = b - sum(x * y for x, y in zip(a[:-1], prefix))
            if a[-1] < 0:
                lower.append(remainder / a[-1])
            elif a[-1] > 0:
                upper.append(remainder / a[-1])
        value = max([F(0)] + lower)
        if upper:
            value = min(value, min(upper))
        assert not lower or value >= max(lower)
        return prefix + [value], None

    solution, failure = eliminate(rows, dimensions)
    if solution is not None:
        assert all(sum(x * y for x, y in zip(a, solution)) <= b for a, b, _ in rows)
    else:
        assert all(w >= 0 for w in failure.values())
        assert all(sum(failure.get(i, 0) * row[0][j] for i, row in enumerate(rows)) == 0
                   for j in range(dimensions))
        assert sum(failure.get(i, 0) * row[1] for i, row in enumerate(rows)) < 0
    return solution, failure


def synthesize(components):
    assert all(len(v) == 3 and all(isinstance(x, int) and x >= 0 for x in v)
               for item in components for v in [item["base"]] + item["periods"])
    unbounded = [i for i, item in enumerate(components) if any(v[0] for v in item["periods"])]
    bounded = [item for i, item in enumerate(components) if i not in unbounded]
    N = 1 + max((item["base"][0] for item in bounded), default=-1)
    gs = {i: reduce(gcd, (v[0] for v in components[i]["periods"]), 0) for i in unbounded}
    M = lcm(*gs.values()) if gs else 1
    names = ["alpha"] + [f"beta_{r}" for r in range(M)] + [f"gamma_{n}" for n in range(N)]
    size = len(names)
    rows, labels = [], []

    def leq(coeff, rhs, label):
        i = len(rows)
        rows.append((tuple(map(F, coeff)), F(rhs), {i: F(1)}))
        labels.append(label)

    def equation(coeff, rhs, label):
        leq(coeff, rhs, label + " <=")
        leq([-x for x in coeff], -rhs, label + " >=")

    for j, name in enumerate(names):
        v = [0] * size
        v[j] = 1
        leq(v, 1, name + " upper")
        v[j] = -1
        leq(v, 1 if 1 <= j <= M else 0, name + " lower")
    for r in range(M):
        for i in unbounded:
            item = components[i]
            if (r - item["base"][0]) % gs[i]:
                continue
            growing = any(v[2] > 0 for v in item["periods"])
            for k, (a, b, h) in enumerate(item["periods"]):
                coeff = [0] * size
                coeff[0], coeff[1 + r] = a, b
                label = f"residue {r} component {i} period {k}"
                if growing:
                    equation(coeff, h, label)
                else:
                    leq(coeff, 0, label)
    small_fibers = [fibers(components, n) for n in range(N)]
    for n, items in enumerate(small_fibers):
        for i, item in enumerate(items):
            for k, (_, m, h) in enumerate(item["periods"]):
                coeff = [0] * size
                coeff[1 + M + n] = m
                equation(coeff, h, f"fiber {n} component {i} period {k}")
    solution, failure = solve_exact(rows, size)
    result = {"components": components, "N": N, "M": M, "variables": names,
              "constraints": [{"coeff": list(a), "rhs": b, "label": labels[i]}
                              for i, (a, b, _) in enumerate(rows)]}
    if solution is None:
        result.update({"synchronizable": False, "infeasibility_multipliers": failure})
        return result
    alpha = solution[0]
    betas = solution[1:1 + M]
    gammas = solution[1 + M:]
    bases = [components[i]["base"] for i in unbounded]
    bases += [item["base"] for items in small_fibers for item in items]
    bound = max((sum(v) for v in bases), default=0) + 2
    state_bound = N + lcm(M, alpha.denominator) + sum(v.denominator for v in betas + gammas) + 1
    result.update({"synchronizable": True, "solution": dict(zip(names, solution)),
                   "skew_construction_bound": bound, "observer_state_construction_bound": state_bound})
    return result


def floor(x):
    return x.numerator // x.denominator


def clock_height(solution, N, M, n, m):
    alpha = solution["alpha"]
    beta = solution[f"gamma_{n}"] if n < N else solution[f"beta_{n % M}"]
    H = 0
    for j in range(1, n + 1):
        c = floor(alpha * j) - floor(alpha * (j - 1))
        assert c in (0, 1)
        H = max(0, H + c)
    for j in range(1, m + 1):
        c = floor(beta * j) - floor(beta * (j - 1))
        assert c in (-1, 0, 1)
        H = max(0, H + c)
    assert H == max(0, floor(alpha * n) + floor(beta * m))
    return H


def sample_points(components, limit):
    points = set()
    for item in components:
        frontier = [item["base"]]
        seen = set()
        while frontier:
            point = frontier.pop()
            if point in seen or point[0] > limit or point[1] > limit or point[2] > 2 * limit:
                continue
            seen.add(point)
            points.add(point)
            for v in item["periods"]:
                nxt = tuple(a + b for a, b in zip(point, v))
                if nxt != point:
                    frontier.append(nxt)
    return points


def main():
    examples = {
        "saturated_descent": (True, [
            linear((0, 0, 0), (1, 0, 1), (1, 1, 0)),
            linear((0, 0, 0), (1, 1, 0), (0, 1, 0))]),
        "valley": (False, [
            linear((0, 0, 0), (1, 0, 1), (1, 1, 0)),
            linear((0, 0, 0), (1, 1, 0), (0, 1, 1))]),
        "parity_dependent_tail": (True, [
            linear((0, 0, 0), (2, 0, 2), (2, 2, 0)),
            linear((2, 1, 1), (2, 0, 2), (2, 2, 0)),
            linear((0, 0, 0), (2, 2, 0), (0, 1, 0)),
            linear((1, 0, 1), (2, 0, 2), (0, 1, 1))]),
        "finite_exception": (True, [
            linear((1, 0, 1), (1, 0, 1), (0, 1, 0)),
            linear((0, 0, 0), (0, 1, 1))]),
        "half_rate": (True, [
            linear((0, 0, 0), (2, 0, 1), (0, 1, 0)),
            linear((1, 0, 0), (2, 0, 1), (0, 1, 0))]),
        "unbounded_same_word_spread": (False, [
            linear((0, 0, 0), (1, 0, 1), (1, 0, 0), (0, 1, 0))]),
        "bounded_ambiguity": (True, [
            linear((0, 0, 0), (1, 0, 1), (0, 1, 0)),
            linear((1, 0, 0), (1, 0, 1), (0, 1, 0))])}
    formulas = {
        "saturated_descent": lambda n, m: {max(0, n - m)},
        "valley": lambda n, m: {abs(n - m)},
        "parity_dependent_tail": lambda n, m: {n + m if n % 2 else max(0, n - m)},
        "finite_exception": lambda n, m: {n if n else m},
        "half_rate": lambda n, m: {n // 2},
        "unbounded_same_word_spread": lambda n, m: set(range(n + 1)),
        "bounded_ambiguity": lambda n, m: {n} if n == 0 else {n, n - 1}}
    results = {}
    for name, (expected, components) in examples.items():
        result = synthesize(components)
        assert result["synchronizable"] == expected
        points = sample_points(components, 24)
        for n in range(25):
            for m in range(25):
                assert {h for nn, mm, h in points if (nn, mm) == (n, m)} == formulas[name](n, m)
        if expected:
            errors = [abs(h - clock_height(result["solution"], result["N"], result["M"], n, m))
                      for n, m, h in points]
            result["finite_audit_maximum_error"] = max(errors)
            assert max(errors) <= result["skew_construction_bound"]
        result["finite_audit_points"] = len(points)
        results[name] = result
    print(json.dumps({"status": "passed", "arithmetic": "exact rational Fourier-Motzkin",
                      "scope": "Supplied semilinear relations only; finite source-formula checks n,m<=24.",
                      "positive_witness_check": "Every original rational inequality rechecked.",
                      "negative_witness_check": "Nonnegative original-row combination gives 0 <= a negative rational.",
                      "not_implemented": ["PDA-to-semilinear conversion", "scaffold compiler"],
                      "formal_verification": "none", "examples": results},
                     default=lambda x: str(x) if isinstance(x, F) else x, indent=2))


if __name__ == "__main__":
    main()