def METHOD_NAME(self, fields, rng, rational_function):
    phi = fields[-1]
    inv_sqrt_M = self.inv_sqrt_matrix(fields, rational_function)
    eta = g.lattice(phi)
    rng.cnormal(eta, sigma=2.0**-0.5)  # 1/sqrt(2)
    phi @= inv_sqrt_M * eta
    return g.norm2(eta)