import math

def h(state, wire, qubits):
    mask = 1 << (qubits - 1 - wire)
    output = state[:]
    scale = math.sqrt(0.5)
    for basis in range(len(state)):
        if basis & mask == 0:
            paired = basis | mask
            output[basis] = (state[basis] + state[paired]) * scale
            output[paired] = (state[basis] - state[paired]) * scale
    return output

def pauli(state, wire, qubits, gate):
    mask = 1 << (qubits - 1 - wire)
    output = [0j] * len(state)
    for basis, amplitude in enumerate(state):
        if gate == "X":
            output[basis ^ mask] += amplitude
        elif gate == "Z":
            output[basis] += -amplitude if basis & mask else amplitude
        else:
            raise ValueError("unsupported Pauli gate")
    return output

def cnot(state, control, target, qubits):
    if control == target:
        raise ValueError("control and target must differ")
    control_mask = 1 << (qubits - 1 - control)
    target_mask = 1 << (qubits - 1 - target)
    output = [0j] * len(state)
    for basis, amplitude in enumerate(state):
        output[basis ^ target_mask if basis & control_mask else basis] += amplitude
    return output

def distribution(state):
    qubits = len(state).bit_length() - 1
    return {format(index, f"0{qubits}b"): round(abs(amplitude) ** 2, 12)
            for index, amplitude in enumerate(state) if abs(amplitude) > 1e-12}

def in_x_basis(state, qubits):
    for wire in range(qubits):
        state = h(state, wire, qubits)
    return state

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

bell = [1.0, 0j, 0j, 0j]
bell = cnot(h(bell, 0, 2), 0, 1, 2)
messages = ((0, 0), (0, 1), (1, 0), (1, 1))
bell_states = []
for phase_bit, flip_bit in messages:
    state = bell[:]
    if flip_bit:
        state = pauli(state, 0, 2, "X")
    if phase_bit:
        state = pauli(state, 0, 2, "Z")
    bell_states.append(state)
max_overlap = 0.0
for left_index, left in enumerate(bell_states):
    for right_index, right in enumerate(bell_states):
        overlap = abs(inner(left, right))
        if left_index != right_index:
            max_overlap = max(max_overlap, overlap)
        assert abs(overlap - (1.0 if left_index == right_index else 0.0)) < 1e-12
assert distribution(bell_states[0]) == {"00": 0.5, "11": 0.5}
assert distribution(bell_states[1]) == {"01": 0.5, "10": 0.5}

ghz = [1.0] + [0j] * 7
ghz = cnot(cnot(h(ghz, 0, 3), 0, 1, 3), 0, 2, 3)
assert distribution(ghz) == {"000": 0.5, "111": 0.5}
ghz_x = distribution(in_x_basis(ghz, 3))
expected_even = {bits for bits in (format(index, "03b") for index in range(8)) if bits.count("1") % 2 == 0}
assert set(ghz_x) == expected_even and all(abs(value - 0.25) < 1e-12 for value in ghz_x.values())

zero_zero = [1.0, 0j, 0j, 0j]
one_one = [0j, 0j, 0j, 1.0]
mixture_x = {}
for component in (zero_zero, one_one):
    for outcome, probability in distribution(in_x_basis(component, 2)).items():
        mixture_x[outcome] = mixture_x.get(outcome, 0.0) + 0.5 * probability
bell_x = distribution(in_x_basis(bell, 2))
mixture_anticorrelation = mixture_x["01"] + mixture_x["10"]
assert bell_x == {"00": 0.5, "11": 0.5}
assert abs(mixture_anticorrelation - 0.5) < 1e-12
phase_mutation_x = distribution(in_x_basis(pauli(bell, 0, 2, "Z"), 2))
assert phase_mutation_x == {"01": 0.5, "10": 0.5}

invalid_control_rejected = False
try:
    cnot(bell, 0, 0, 2)
except ValueError:
    invalid_control_rejected = True
assert invalid_control_rejected
print(f"PASS: 20 Bell/GHZ harness verifies {len(bell_states)} orthogonal Bell states (max cross-overlap={max_overlap:.1e}), GHZ X-even support={len(ghz_x)}, classical-mixture X anticorrelation={mixture_anticorrelation:.3f}")