import math
import random

def inner(left, right):
    return sum(a.conjugate() * b for a, b in zip(left, right))

def matvec(matrix, vector):
    return tuple(sum(matrix[row][column] * vector[column] for column in range(2)) for row in range(2))

def normalize(vector):
    norm = math.sqrt(inner(vector, vector).real)
    return tuple(value / norm for value in vector)

def eigh2(matrix):
    if abs(matrix[0][0].imag) > 1e-12 or abs(matrix[1][1].imag) > 1e-12 or abs(matrix[1][0] - matrix[0][1].conjugate()) > 1e-12:
        raise ValueError("observable must be Hermitian")
    a, d, off = matrix[0][0].real, matrix[1][1].real, matrix[0][1]
    center = (a + d) / 2
    radius = math.sqrt(((a - d) / 2) ** 2 + abs(off) ** 2)
    values = (center + radius, center - radius)
    if radius < 1e-15:
        vectors = ((1 + 0j, 0j), (0j, 1 + 0j))
    else:
        vectors = tuple(normalize((off, value - a)) if abs(off) > 1e-15
                        else ((1 + 0j, 0j) if abs(value - a) < 1e-12 else (0j, 1 + 0j))
                        for value in values)
    return values, vectors

def born_probabilities(state, eigenvectors):
    return tuple(abs(inner(vector, state)) ** 2 for vector in eigenvectors)

def seeded_estimate(values, probabilities, shots, seed):
    if shots < 1:
        raise ValueError("shots must be positive")
    rng = random.Random(seed)
    samples = [values[0] if rng.random() < probabilities[0] else values[1] for _ in range(shots)]
    mean = sum(samples) / shots
    variance = sum((sample - mean) ** 2 for sample in samples) / shots
    return mean, math.sqrt(variance / shots)

scale = math.sqrt(0.5)
observables = (
    ((1 + 0j, 0j), (0j, -1 + 0j)),
    ((0j, 1 + 0j), (1 + 0j, 0j)),
    ((1 + 0j, 1j), (-1j, -1 + 0j)),
    ((2 + 0j, 0j), (0j, 2 + 0j)),
)
states = ((1, 0), (scale, scale), (scale, 1j * scale), (0.6, 0.8j))
max_residual = 0.0
max_sampling_z = 0.0
measurement_cases = 0
for observable_index, observable in enumerate(observables):
    values, eigenvectors = eigh2(observable)
    assert all(isinstance(value, float) for value in values)
    assert abs(inner(eigenvectors[0], eigenvectors[1])) < 1e-12 or abs(values[0] - values[1]) < 1e-12
    for value, vector in zip(values, eigenvectors):
        residual = max(abs(actual - value * expected) for actual, expected in zip(matvec(observable, vector), vector))
        max_residual = max(max_residual, residual)
    for state_index, state in enumerate(states):
        probabilities = born_probabilities(state, eigenvectors)
        assert abs(sum(probabilities) - 1.0) < 1e-12
        analytic = inner(state, matvec(observable, state)).real
        spectral = sum(value * probability for value, probability in zip(values, probabilities))
        assert abs(analytic - spectral) < 1e-12
        estimate, stderr = seeded_estimate(values, probabilities, 30000, 1400 + 10 * observable_index + state_index)
        allowance = 6 * stderr + 1e-12
        assert abs(estimate - analytic) <= allowance
        max_sampling_z = max(max_sampling_z, abs(estimate - analytic) / max(stderr, 1e-12))
        measurement_cases += 1
assert max_residual < 1e-12

nonhermitian_rejected = False
try:
    eigh2(((0j, 1j), (1j, 0j)))
except ValueError:
    nonhermitian_rejected = True
assert nonhermitian_rejected
zero_shots_rejected = False
try:
    seeded_estimate((1, -1), (0.5, 0.5), 0, 1)
except ValueError:
    zero_shots_rejected = True
assert zero_shots_rejected
print(f"PASS: 14 observable notebook solves {len(observables)} Hermitian spectra and {measurement_cases} seeded estimates; last estimate={estimate:.4f}, max eigen-residual={max_residual:.2e}, max sampling z={max_sampling_z:.2f}")