"""Every tolerance the corpus holds an engine to, each derived from that engine's arithmetic.

None is fitted to a result. Each is a first-order worst-case bound, and so is loose: an engine
that passes is within it, and the run also records how far inside it each engine is. Where an
assumption enters (a library's backward-error constant, which the library does not state), the
function says so and states the number it takes.

Floating point (unit roundoff u: 2^-53 in double, 2^-24 in single). One gate acting on w qubits
at once, applied to a normalised state, errs in the Euclidean norm by at most C(w)·u:

  C(w) = (2√2 + (2^w − 1) + 2 + π) · 2^(w/2)

- 2√2 u: a complex product of two rounded numbers (Higham, Accuracy and Stability of Numerical
  Algorithms, 2nd ed., §3.6);
- 2^w − 1: the additions of a row of 2^w terms, summed in any order;
- 2 + π: the matrix entries, each within 2u from its sine and cosine, and moved by at most π·u
  by the rounding of an angle of at most 2π;
- 2^(w/2): the norm of the matrix of absolute values, at most √(2^w) for a unitary.

Errors of successive gates add (each gate is unitary, so it carries the error before it unchanged
in norm). A fused gate of width w costs C(w) for each gate fused into it, since its matrix is a
product of theirs computed in the same arithmetic.
"""
import math

U_DOUBLE = 2.0**-53
U_SINGLE = 2.0**-24


def C(w):
    return (2 * math.sqrt(2) + (2**w - 1) + 2 + math.pi) * 2 ** (w / 2)


def gate_count(steps):
    n = 0
    for s in steps:
        for x in s.get("then", [s]):
            if "g" in x:
                n += 1
    return n


def _first_order(total):
    # The bound's own second-order term: (1 + x)·x for x below 1.
    return total * (1 + total) if total < 1 else float("inf")


def double_state(gates, n, width):
    return _first_order(gates * C(width) * U_DOUBLE)


def single_state(gates, n, width):
    """Single precision (complex64 amplitudes) with the same model and u = 2^-24."""
    return _first_order(gates * C(width) * U_SINGLE)


def contraction(steps, width):
    """A tensor network contracted pairwise: a contraction with an index of 2^width values is a
    row sum of that many terms, so each step costs C(width). The model assumes, as for a gate, that
    the norm of an intermediate's matrix of absolute values is at most √(2^width) times its own;
    for a network of unitaries contracted in a sensible order that holds, and the measured
    errors (recorded beside it) are the check on the assumption."""
    return _first_order(steps * C(width) * U_DOUBLE)


DD_TOLERANCE = 2.0**-52 * 1024  # mqt.core's dd::RealNumber::eps, read from its header


def decision_diagram(gates, n):
    """A decision-diagram engine merges edge weights that lie within its tolerance τ of an entry
    of its table, so each weight moves by at most τ; an amplitude is a product of n + 1 weights,
    so each gate moves it by at most (n + 1)·τ relatively, and the double arithmetic adds the
    state vector's bound."""
    return _first_order(gates * (n + 1) * DD_TOLERANCE) + double_state(gates, n, 2)


LAPACK_CONSTANT = 16  # p(m, n) in LAPACK's SVD backward error p(m, n)·u·||A||, taken as 16·m


def svd_steps(steps):
    """The two-site updates a matrix product state makes for a circuit: a gate on qubits d apart is
    routed by d − 1 swaps there and d − 1 back, each an update, and the gate itself one more."""
    total = 0
    for s in steps:
        for x in s.get("then", [s]):
            if "g" in x and len(x["q"]) >= 2:
                total += 2 * (max(x["q"]) - min(x["q"])) - 1
    return total


def mps(updates, chi, truncation):
    """A matrix product state: each two-site update is a contraction and an SVD of a 2χ by 2χ
    matrix, whose backward error LAPACK states as p(m, n)·u with p "a modestly growing function"
    of the shape; we take p = 16·2χ. `truncation` bounds what the engine's truncations moved the
    state: a truncation that drops weight w moves a normalised state by at most 2√w (√w dropped,
    at most w more by renormalising), and the unitary steps after it carry that unchanged, so
    these add. The weights are the ones the engine reports (truncation_from)."""
    svd = updates * (C(2) + LAPACK_CONSTANT * 2 * chi) * U_DOUBLE
    return _first_order(svd) + truncation


def truncation_from(weights=None, total=None, count=0):
    """2·Σ√wᵢ from each truncation's weight when the engine reports them one by one; from a total
    W over `count` truncations, 2·√(count·W), which bounds that sum (Cauchy–Schwarz)."""
    if weights is not None:
        return 2 * sum(math.sqrt(max(w, 0.0)) for w in weights)
    return 2 * math.sqrt(max(count, 0) * max(total or 0.0, 0.0))


def density(gates, n):
    """A density matrix: each gate is applied on both sides (2·C(2)), each depolarizing channel
    is a sum of up to 16 Pauli conjugations (16·C(2)); an expectation value of a Pauli string
    errs by at most √(2^n) times the error of ρ in the Frobenius norm."""
    return math.sqrt(2**n) * gates * (2 + 16) * C(2) * U_DOUBLE


def dem(noise_sites):
    """A detector error model's probabilities. A channel's per-Pauli rate, (1 − √(1 − 4p/3))/2 or
    its two-qubit form, is within 3u absolutely; merging two mechanisms, p + q − 2pq, adds at most
    3u more; a mechanism merges at most every noise channel of the circuit. Two engines differ by at
    most twice that: 2 · noise_sites · 6u, written 16u per site."""
    return noise_sites * 16 * U_DOUBLE


# Fixed point: the extension's own bound (quantum-sim.md §5). Per inexact gate, every component
# it writes is within 1 + √2/2 units of 2^-30 of the ideal gate applied to the same values: one
# unit from two products rounded half a unit each, √2/2 from table constants within half a unit.
PER_COMPONENT = 1 + math.sqrt(2) / 2
UNIT = 2.0**-30


def fixed_point_op(n, base, controls, param):
    """The Euclidean bound, in units, for one op of the bit-exact op-log: base gate codes are
    quantum-sim.md's (§3.0): I X Y Z H S S† T T† P = 0..9."""
    c = len(controls)
    if base in (0, 1, 2, 3, 5, 6):
        return 0.0  # entries 0, ±1, ±i: exact
    if base == 9 and param <= 2:
        return 0.0  # P(1) = -1 and P(2) = i: exact
    if base == 4:
        # H writes both amplitudes of each pair: four components per pair.
        return PER_COMPONENT * math.sqrt(4 * 2 ** (n - 1 - c))
    # T, T† and P(k) write one amplitude of each pair: two components.
    return PER_COMPONENT * math.sqrt(2 * 2 ** (n - 1 - c))


def fixed_point_bound(n, ops):
    """The bit-exact method's distance from the ideal state, bounded by its published per-gate
    bound summed over the op-log it ran, each gate carrying the error before it: a fixed-point
    gate's matrix is within 2·2^-30 of a unitary, so it can grow that error by at most that."""
    total = 0.0
    inexact = 0
    for base, controls, target, param in ops:
        b = fixed_point_op(n, base, controls, param)
        if b:
            inexact += 1
        total += b
    return total * UNIT * (1 + 2 * UNIT) ** inexact


def l1_distance(shots, outcomes, alpha):
    """ε with P(Σ|X_k/shots − p_k| ≥ ε) ≤ alpha for counts X of `shots` draws over `outcomes`
    categories: Weissman, Ordentlich, Seroussi, Verdú and Weinberger (2003), (2^K − 2)·exp(−shots·ε²/2).
    The L1 distance is twice the total variation. No normal approximation. It holds the whole
    distribution at once, where the per-outcome test below holds each count alone."""
    return math.sqrt(2 * (outcomes * math.log(2) + math.log(1 / alpha)) / shots)


def bernstein(shots, p, alpha):
    """t with P(|X − shots·p| ≥ t) ≤ alpha for X binomial(shots, p), from Bernstein's
    inequality 2·exp(−t² / (2(V + t/3))), V = shots·p·(1 − p). No normal approximation."""
    b = math.log(2 / alpha)
    v = shots * p * (1 - p)
    return b / 3 + math.sqrt(b * b / 9 + 2 * b * v)


# ------------------------------------------------------------------ the engine's other methods


def pauli_rounding(rotations, terms, coeff):
    """Pauli propagation's rounding, beside its truncation bound. A rotation maps each pair of
    strings (Q, Q') it splits by a 2x2 rotation of their coefficients, so the coefficients' 2-norm
    is unchanged by the exact map (it is the observable's Hilbert-Schmidt norm, |coeff|), and each
    rotation adds at most 16u of it: each new coefficient is two products and a sum (3u of
    |c| + |c'|), its sine and cosine within an ulp each and the angle within u·2π, so (5 + 2π)u per
    coefficient, √2 of that over the pair. The value sums at most `terms` coefficients, which by
    Cauchy-Schwarz is √terms times their 2-norm, with (terms − 1)u of rounding on the sum. So the
    value is within |coeff|·√terms·(16·rotations + terms)·u of the exact propagation's."""
    return abs(coeff) * math.sqrt(max(terms, 1)) * (16 * rotations + max(terms, 1)) * U_DOUBLE


def pauli_bound_rounding(bound, terms, rotations):
    """The truncation bound is a sum of every dropped magnitude, at most terms·rotations of them,
    each sum rounded: the bound as computed is within that many u of its value."""
    return bound * max(terms, 1) * max(rotations, 1) * U_DOUBLE


SEC_PI_8 = 1 / math.cos(math.pi / 8)


def stabrank_exact(t, n, coeff=1.0):
    """The stabilizer-rank sum, exactly: every term is 2^(−k/2)·e^(iπj/8), counted exactly per
    (k, j), and the counts evaluated once in double precision. Each of the at most 2n + 1 powers k
    sums eight counts times a constant rounded within u/2 (8u), the powers are summed (2n + 1)u,
    each power of 1/√2 is within u, and the scale (1/(2cos π/8))^t is within 2tu; every error is
    relative to the sum of the terms' magnitudes, at most 2^t (one for each branch, each amplitude
    of a stabilizer state at most 1) times the scale: sec(π/8)^t. So (2t + 2n + 13)·u·sec(π/8)^t."""
    return abs(coeff) * (2 * t + 2 * n + 13) * U_DOUBLE * SEC_PI_8**t


def cutrank(n):
    """The cut-rank sum of an amplitude: 2^n terms, each a product of n one-qubit weights
    W_q/√2 (each entry of magnitude at most 1/√2, computed as a product of two 2x2 matrices and a
    scaling, within 6u) and signs ±1, so each term's magnitude is at most 2^(−n/2) and its value
    within (2√2 + 6)·n·u ≤ 9nu of its products; the passes add each term into at most 2n sums (a
    pass per qubit, a merge per dimension lost), 2nu more. The sum of the terms' magnitudes is
    at most 2^n·2^(−n/2) = 2^(n/2), so an amplitude is within (11n + 2)·u·2^(n/2)."""
    return (11 * n + 2) * U_DOUBLE * 2 ** (n / 2)


def hoeffding(samples, spread, alpha):
    """t with P(|mean − E| ≥ t) ≤ alpha for the mean of `samples` independent values each within an
    interval of width `spread` (Hoeffding's inequality, 2·exp(−2·samples·t²/spread²)). No normal
    approximation."""
    return spread * math.sqrt(math.log(2 / alpha) / (2 * samples))


def binomial_upper(trials, p, alpha):
    """The least x with P(X ≥ x) ≤ alpha for X binomial(trials, p): a count of at least x
    happens with probability at most alpha when each trial succeeds with probability at most p."""
    tail = 0.0
    for x in range(trials, -1, -1):
        term = math.comb(trials, x) * p**x * (1 - p) ** (trials - x)
        if tail + term > alpha:
            return x + 1
        tail += term
    return 0


def density_steps(n, steps):
    """A density matrix run as the engine runs one (its rows and columns a state of 2n qubits):
    `steps` superoperators, each a gate on the rows or its conjugate on the columns, or a channel
    that is a sum of at most 16 conjugations, each within 16·C(4)·u at the fusion width of 4 the
    run takes; the error of ρ in the Frobenius norm is at most their sum, and an expectation value
    of a Pauli string errs by at most √(2^n) times it."""
    return math.sqrt(2**n) * steps * 16 * C(4) * U_DOUBLE


def noisy_steps(gates, noise):
    """The superoperators a noisy run applies (wai-quantum quantum_noisy::DensityMatrix::evolve):
    two for each gate, a two-qubit gate's depolarizing channel, and for each qubit a gate touches
    its one-qubit depolarizing channel (a gate on one qubit or three or more), dephasing and
    damping. `gates` is a list of qubit counts."""
    total = 0
    for k in gates:
        total += 2
        if k == 2 and noise.get("depolarizing_2q"):
            total += 1
        for _ in range(k):
            total += (1 if k != 2 and noise.get("depolarizing_1q") else 0) + (1 if noise.get("dephasing") else 0) + (1 if noise.get("amplitude_damping") else 0)
    return total


def matching_rounding(fired, detectors):
    """Exact matching weighs each edge ln((1 − p)/p) rounded to 1/16384, so a path of at most
    detectors + 1 edges is within (detectors + 1)/32768 of its exact weight; a correction pairs
    the fired detectors in at most `fired` paths, and the one it finds least on its integers is
    within twice that of the least on exact weights."""
    return max(fired, 1) * (detectors + 1) / 16384