import numpy as np

I = np.matrix([[1, 0], [0, 1]])
Z = np.matrix([[1, 0], [0, -1]])
X = np.matrix([[0, 1], [1, 0]])

def Mi(nqbits, i, M):
    result = 1
    for j in range(nqbits):
        if(j != i):
            result = np.kron(result, I)
        else:
            result = np.kron(result, M)
    return result



def H_interaction(nqbits):
    interaction_terms = [Mi(nqbits, i, Z) @ Mi(nqbits, i+1, Z) for i in range(nqbits)]
    return sum(interaction_terms)

def H_field(nqbits, g):
    field_terms = [g*Mi(nqbits, i, X) for i in range(nqbits)]
    return sum(field_terms)

def H(nqbits, g):
    return (-H_interaction(nqbits) + H_field(nqbits, g)).real