from math import sqrt
H = ((1/sqrt(2), 1/sqrt(2)), (1/sqrt(2), -1/sqrt(2)))
state = [1+0j, 0j, 0j, 0j]

def one_qubit(vector, gate, target, n=2):
    out = vector[:]
    mask = 1 << (n - 1 - target)
    for i in range(len(vector)):
        if i & mask: continue
        j = i | mask; a, b = vector[i], vector[j]
        out[i] = gate[0][0]*a + gate[0][1]*b
        out[j] = gate[1][0]*a + gate[1][1]*b
    return out

def cnot(vector, control, target, n=2):
    out = [0j] * len(vector)
    cm = 1 << (n - 1 - control); tm = 1 << (n - 1 - target)
    for i, amplitude in enumerate(vector): out[i ^ tm if i & cm else i] += amplitude
    return out

product_state = one_qubit(state, H, 0)
state = cnot(product_state, 0, 1)
missing_entangler = cnot(product_state, 1, 0)
expected = [1/sqrt(2), 0, 0, 1/sqrt(2)]
assert all(abs(a-b) <= 1e-12 for a,b in zip(state, expected))
assert abs(sum(abs(a)**2 for a in state)-1) <= 1e-12
assert missing_entangler == product_state and missing_entangler != state
print(f"PASS: 37 simulator Bell_support={[i for i,a in enumerate(state) if abs(a) > 1e-12]} missing_entangler_support={[i for i,a in enumerate(missing_entangler) if abs(a) > 1e-12]}")