const int N = 10;
const double Bx = 2;

const double total_T = 3.0;  // total duration of adiabatic evolution
const int M = 5;  // number of Trotter steps

// elements of Bz and J are randomly generated from uniform(-2, 2)
const double Bz[] = {.3, -.9, -1.2, 1.1, 1.5, .6, -1.6, .3, -1.3, 1.9};
const double J[] = {-1.5, 1.3, -1.8, -1.3, -.6, 1.9, 1.1, -1.3, .4};

module CZ(qbit q1, qbit q2, double phi) {
    Rz(q2, 0.5*phi);
    CNOT(q1, q2);
    Rz(q2, -0.5*phi);
    CNOT(q1, q2);
}

module ZcrossZ(qbit q1, qbit q2, double phi) {
    // printf("performing sigma_z cross sigma_z Hamiltonian\n");
    Rz(q1, phi);
    Rz(q2, -phi);
    CZ(q1, q2, -2.0*phi);
}

module initialize(qbit reg[N]) {
  int n;
  for(n=0; n < N; n++) {
    PrepZ(reg[n], 0);
    H(reg[n]);
  }  
}

module red_hamiltonian(qbit reg[N], int m) {
  int n;
  for (n = 0; n < N-1; n += 2) {  // red pairs
    double phi = J[n] * (2.0*m - 1) / M;        
    ZcrossZ(reg[n], reg[n + 1], phi);
  }
}

module black_hamiltonian(qbit reg[N], int m) {
  int n;
  for (n = 1; n < N-1; n += 2) {  // black pairs
    double phi = J[n] * (2.0*m - 1) / M;        
    ZcrossZ(reg[n], reg[n + 1], phi);
  }
}

module Bz_hamiltonian(qbit reg[N], int m) {
  int n;
  for (n = 0; n < N; n++) {
    double theta1 = (1.0 - (2.0*m-1)/M) * -2 * Bx * total_T / M;
    double theta2 = (1.0 - (2.0*m-1)/M) * -2 * Bz[n] * total_T / M;     
    H(reg[n]);
    Rz(reg[n], theta1);
    H(reg[n]);
    Rz(reg[n], theta2);    
  }
}

module measure(qbit reg[N]) {
  int n;
  for (n = 0; n < N; n++) {
    MeasZ(reg[0]);
  }  
}

int main() {
  qbit reg[N];
  
  // Initialize all qubits to |+> state
  initialize(reg);
  
  int m;
  for (m=1; m <= M; m++) {
    // Z-Z Hamiltonian, executed in parallel across qbits
    red_hamiltonian(reg, m);    
    black_hamiltonian(reg, m);
    // Bz Hamiltonian, executed in parallel across qbits   
    Bz_hamiltonian(reg, m);
  }

  // Measure all qubits in Z basis
  measure(reg);

  return 0;
}