"""Reproduce bounded checks for the Hawking-at-50 research companion.

Run: python hawking-phase-and-time-audit.py
Python 3.10+ standard library only. Emits JSON; assertions fail on disagreement.

1. QMM 2026, Eq. (8), doi:10.1007/s11128-026-05292-0.
   Q0 is the arbitrary input/control, Q1 and Q2 start in |0>. Use precisely
   CRY(0->1), Rz(1,phi), CSWAP(control=0,targets=1,2), Rz(2,-phi).
   Trace out Q0,Q1 and inspect Q2 without an extra decoder or side information.
   This is an ideal mathematical channel audit, not IBM hardware reproduction,
   an audit of the authors' QASM, or a refutation of all QMM constructions.

2. Bianchi et al., arXiv:2605.03922v2, Eqs. (3),(5),(21).
   Check the stated energy/entropy/time inequality using smooth example flux
   profiles within its assumptions. These are not astrophysical predictions.

3. QMM 2024, Sec. 2.1.1 and Eq. (3), doi:10.3390/e26121039.
   Illustrate the trace obstruction to exact canonical commutation on a finite
   Hilbert space. The general proof is trace(AB)=trace(BA); the example does not
   replace that proof. Approximate or collective infinite-space operators have
   to be specified separately. This is not a blanket verdict on QMM.

Independent AI-assisted appraisal, u/NeuronsToNirvana, 6 November 2026.
"""

import cmath
import json
import math
import random


def cycle(v, theta, phi):
    """State-vector gate application; basis order Q0,Q1,Q2, Q0 most significant."""
    w = list(v)
    c, s = math.cos(theta / 2), math.sin(theta / 2)
    for q2 in (0, 1):
        a, b = 4 + q2, 6 + q2
        w[a], w[b] = c * w[a] - s * w[b], s * w[a] + c * w[b]
    w = [z * cmath.exp((1j if (i >> 1) & 1 else -1j) * phi / 2)
         for i, z in enumerate(w)]
    w[5], w[6] = w[6], w[5]
    return [z * cmath.exp((-1j if i & 1 else 1j) * phi / 2)
            for i, z in enumerate(w)]


def output_density(v):
    return [[sum(v[2 * k + i] * v[2 * k + j].conjugate() for k in range(4))
             for j in range(2)] for i in range(2)]


def initial(psi):
    return [psi[0], 0j, 0j, 0j, psi[1], 0j, 0j, 0j]


def analytic_output(psi, theta):
    p = abs(psi[1]) ** 2
    c, s = math.cos(theta / 2), math.sin(theta / 2)
    return [[1 - p + p * c * c, p * c * s], [p * c * s, p * s * s]]


def difference(a, b):
    return max(abs(x - y) for row_a, row_b in zip(a, b)
               for x, y in zip(row_a, row_b))


def swap_input_output(v):
    w = [0j] * 8
    for i, z in enumerate(v):
        q0, q1, q2 = i >> 2, (i >> 1) & 1, i & 1
        w[4 * q2 + 2 * q1 + q0] = z
    return w


def phase_checks():
    rng = random.Random(1729)  # Reproducible seed, with no physical significance.
    errors = {"analytic_channel": 0.0, "unitarity": 0.0,
              "global_inverse": 0.0, "swap_baseline": 0.0}
    for _ in range(128):
        theta, phi = rng.uniform(-4 * math.pi, 4 * math.pi), rng.uniform(-4 * math.pi, 4 * math.pi)
        psi = [complex(rng.gauss(0, 1), rng.gauss(0, 1)) for _ in range(2)]
        norm = math.sqrt(sum(abs(z) ** 2 for z in psi))
        psi = [z / norm for z in psi]
        v = initial(psi)
        w = cycle(v, theta, phi)
        errors["analytic_channel"] = max(errors["analytic_channel"],
            difference(output_density(w), analytic_output(psi, theta)))
        columns = [cycle([complex(i == j) for i in range(8)], theta, phi)
                   for j in range(8)]
        for i in range(8):
            for j in range(8):
                inner = sum(columns[i][k].conjugate() * columns[j][k] for k in range(8))
                errors["unitarity"] = max(errors["unitarity"], abs(inner - (i == j)))
        recovered = [sum(columns[j][k].conjugate() * w[k] for k in range(8)) for j in range(8)]
        errors["global_inverse"] = max(errors["global_inverse"], max(abs(x - y) for x, y in zip(v, recovered)))
        target = [[psi[i] * psi[j].conjugate() for j in range(2)] for i in range(2)]
        errors["swap_baseline"] = max(errors["swap_baseline"],
            difference(output_density(swap_input_output(v)), target))
    theta, phi = math.pi / 4, 0.6
    plus = [1 / math.sqrt(2), 1 / math.sqrt(2)]
    minus = [1 / math.sqrt(2), -1 / math.sqrt(2)]
    a = output_density(cycle(initial(plus), theta, phi))
    b = output_density(cycle(initial(minus), theta, phi))
    assert all(e < 1e-12 for e in errors.values()), errors
    assert difference(a, b) < 1e-12
    # Numerically propagate both branches of an entangled reference + input,
    # then trace control and memory. The reference is untouched by the cycle.
    def bell_overlap_after(operation):
        branches = [operation(initial([1 / math.sqrt(2), 0j])),
                    operation(initial([0j, 1 / math.sqrt(2)]))]
        rho = [[sum(branches[i // 2][2 * k + i % 2] *
                    branches[j // 2][2 * k + j % 2].conjugate()
                    for k in range(4)) for j in range(4)] for i in range(4)]
        return ((rho[0][0] + rho[0][3] + rho[3][0] + rho[3][3]) / 2).real
    bell_overlap = bell_overlap_after(lambda v: cycle(v, math.pi, 0.6))
    swap_bell_overlap = bell_overlap_after(swap_input_output)
    assert abs(bell_overlap - 0.5) < 1e-12
    assert abs(swap_bell_overlap - 1.0) < 1e-12
    return {"random_cases": 128, "maximum_errors": errors,
            "theta": theta, "phi": phi,
            "plus_minus_output_difference": difference(a, b),
            "example_output_real": [[z.real for z in row] for row in a],
            "basis_perfect_case_theta_pi_bell_overlap": bell_overlap,
            "swap_baseline_bell_overlap": swap_bell_overlap,
            "scope": "zero memory/output ancillas, output Q2 only, no decoder or side information"}


def time_checks():
    hbar, entropy, energy, alpha = 0.7, 5.1, 2.3, 3 / (48 * math.pi)
    lower = hbar * entropy ** 2 / (16 * math.pi ** 2 * alpha * energy)
    profiles = []
    for shape in (0.0, 0.2, 0.8, -0.6):
        duration = lower * (1 + shape ** 2 / 2)
        amplitude = entropy / (4 * math.pi * alpha * duration)
        steps = 4096
        k = [-amplitude * (1 + shape * math.sin(2 * math.pi * i / steps))
             for i in range(steps + 1)]
        def integral(values):
            return duration / steps * (sum(values[1:-1]) + (values[0] + values[-1]) / 2)
        measured_entropy = -4 * math.pi * alpha * integral(k)
        measured_energy = alpha * hbar * integral([x * x for x in k])
        assert abs(measured_entropy - entropy) < 1e-10
        assert abs(measured_energy - energy) < 1e-10
        assert duration >= lower - 1e-12
        profiles.append({"shape": shape, "duration_over_bound": duration / lower,
                         "entropy_error": abs(measured_entropy - entropy),
                         "energy_error": abs(measured_energy - energy)})
    return {"equation": "tau >= hbar*S_C^2/(16*pi^2*alpha*M_C)",
            "profiles": profiles,
            "scope": "conformal-vacuum s-wave model; positive illustrative inputs; not a lifetime estimate"}


def finite_commutator_check():
    scale = 1 / math.sqrt(2)
    x = [[0j, scale], [scale, 0j]]
    p = [[0j, -1j * scale], [1j * scale, 0j]]
    def product(a, b):
        return [[sum(a[i][k] * b[k][j] for k in range(2))
                 for j in range(2)] for i in range(2)]
    xp, px = product(x, p), product(p, x)
    comm = [[xp[i][j] - px[i][j] for j in range(2)] for i in range(2)]
    trace = comm[0][0] + comm[1][1]
    assert abs(trace) < 1e-12
    return {"dimension": 2, "hbar": 1,
            "trace_commutator": [trace.real, trace.imag],
            "trace_required_by_exact_CCR": [0, 2],
            "scope": "illustration of exact finite-dimensional trace obstruction; no approximate/collective representation specified"}


if __name__ == "__main__":
    print(json.dumps({"qmm_channel": phase_checks(), "purification_bound": time_checks(),
                      "finite_commutator": finite_commutator_check()}, indent=2))