#define PI 3.141592653589793238462643383279502884197
const int N = 4;

// NxN and NxNxNxN rotation angles for single and double excitation operators
// generated in Python as [random.uniform(0, 2*math.pi) for _ in range(N*...*N)]
const double Theta_p_q[N][N] = {3.2182142487054763, 3.7594365739635673, 0.7679973723407487, 5.005325008037422, 1.995481713596621, 0.12747368999445358, 3.346551736496727, 1.5377881080327662, 4.33258219753088, 0.461922013898618, 5.919638122624508, 1.853058773043441, 1.08697553231529, 2.2583944071585247, 1.2285307593407533, 5.190407212320072};
const double Theta_p_q_r_s[N][N][N][N] = {1.3679496193514813, 5.3987021473931485, 3.7662849087964174, 2.191305104063104, 3.215552711333309, 1.486382716317953, 0.2572307363006586, 4.78569957275046, 4.787031144801062, 1.3444389251718336, 0.42688896624930817, 1.8080220007155932, 0.8505376652589783, 5.860116935203846, 2.5943617870182973, 2.1883936951975946, 5.700141765748832, 2.1970559292456135, 3.8775361457772837, 2.229341106472796, 5.271361790475311, 1.0653919596797994, 2.8011268953022563, 2.9190648832098653, 2.4939581042612446, 0.6152344327331015, 3.456377109494123, 0.28204090367162754, 0.774220878411388, 5.407549605682287, 2.5929873060649995, 3.704762570757387, 4.230628804334122, 4.531773693038241, 0.061560287232719335, 0.03533994899234157, 3.645292140442946, 0.26093530847612345, 1.6620362191786842, 4.989004894164902, 2.8525757263708296, 5.678955234001865, 5.298583063509366, 1.5808331714915942, 5.179721443729388, 0.37420467536999713, 2.646598084578072, 5.197633070894743, 4.197227936979659, 1.8612864228449677, 1.5204252304758463, 4.617259521457752, 0.4450404292194258, 5.922360509132177, 1.4705760743296976, 3.148016150025384, 2.4012497414914966, 0.7784759728212509, 3.397665162250324, 5.702197802866682, 5.267081686666654, 5.418615389423272, 0.16303072005059183, 5.889582747385034, 3.9002567920600315, 2.938537761256362, 3.794484658818093, 1.1861105452823373, 3.910890894199179, 5.941276387811751, 2.3183394678014717, 1.2206942962042053, 0.6930426381884766, 3.6655083109747353, 5.51889970382442, 2.67865577854184, 1.527186957340491, 4.866415818215523, 4.250712941931898, 0.12326389592931845, 0.9017092080420118, 3.8827922895013143, 0.7746198586100097, 6.183705019122961, 4.687234603481983, 3.2041512394266185, 6.058124722365269, 1.4120578937655803, 1.752748264782318, 1.1471867473168882, 0.596938775700471, 0.1720446712682702, 5.664830165172676, 1.5836415778105082, 0.6480181950370385, 0.5810089916659814, 1.8719873727895782, 0.492672094311885, 3.228909104112608, 4.194341037782101, 0.3066754200950564, 5.549517743242185, 3.4383093080591425, 1.131006601711999, 1.750860743043907, 1.7737980527245725, 3.470865351091518, 0.4015339001736433, 0.1328788892641507, 3.311076433008854, 3.7029543230587243, 4.075762053550507, 5.353945986667215, 3.1473488845527315, 4.763720270502149, 2.7237039309590902, 3.290358874993329, 1.0282441649725276, 4.3378270125156675, 1.676265070528847, 6.217970805793189, 5.2899160982192095, 3.5974612236524925, 1.3042709286390664, 4.417537717629104, 4.781749341223548, 3.54738768196523, 4.285896329730382, 2.0043186882638753, 5.969716862384182, 0.5845174346070312, 0.25763852990918723, 1.5145245142426182, 4.894083134267615, 2.957543700616756, 6.013650820008401, 1.8994059951567266, 3.0163547932843158, 5.726142673083595, 3.0832655185004203, 3.7944007436345646, 3.7478378158358856, 0.6774444695014695, 3.6746544344459116, 6.069315930151495, 5.1568851190528076, 1.4708555157457441, 2.8082444294049136, 6.047958023987658, 6.166010238312782, 0.5833579691799105, 3.2471790371005524, 4.841950374442266, 4.9625898103288355, 1.7360918528505656, 1.273558641943351, 6.236016881985402, 4.867138293052107, 3.4085033940312432, 5.537924547540386, 5.482420606468638, 3.073471398022417, 3.89847443381525, 0.21071263585971908, 3.4204788002293585, 1.2862371171892573, 1.6680990653613195, 5.576193563539006, 1.4458059644726329, 4.96637989818178, 4.151430191258576, 6.119272262234779, 0.3768389211894691, 3.9197588402550685, 1.3530727526178512, 3.6033549063186565, 5.783847478144294, 3.2128957607612403, 1.8723912594186767, 1.945154965789314, 5.985269903720078, 2.216880944763532, 0.7641833618397538, 2.6544410136344307, 2.7394907973021168, 2.4000687986670197, 1.6341253770942281, 0.6937065188174155, 2.3682699831940512, 0.9211373429746709, 4.11991859258512, 1.2988469800350073, 1.7199685561276337, 0.45984674955422994, 5.789531821585711, 4.281666467797316, 4.734666729386387, 5.831602291735034, 5.2751677932654335, 5.964638996971862, 2.959372906369318, 4.866393326552442, 0.5158107121323209, 5.752075569564825, 4.505837144868799, 0.7529446060146769, 3.5843322473911567, 0.7050794490492214, 0.5778163780739158, 3.8421226289509334, 4.973310497645709, 3.3901400892868585, 4.867794053516748, 4.305914249140959, 2.8696533179420634, 1.77783134745889, 3.879027125517393, 5.723132091219646, 4.460405217879545, 5.79143981135287, 3.6130118744438855, 2.457554912810849, 0.31459583668211744, 0.7184519159766944, 3.218389103278122, 1.599832915960389, 2.504199941780175, 1.4544114866060491, 2.1517464720381634, 2.733828663529408, 3.5986027927980646, 5.226575868382455, 0.9879239990096947, 1.1670274845452207, 3.783803447171394, 5.806177526684089, 4.8548397377780415, 4.7513555274988235, 6.068892664156556, 0.6256289838472665, 1.4795901287467799, 4.471253347039421, 2.2285942887214567, 0.844591345361827, 0.5114549948183631, 1.5989550892182625, 0.9333583354413291, 0.6018205076052627, 4.301476538208704, 2.820005746606613, 4.014584238343068, 3.0093584140015572, 1.8451466370892362, 4.490647193354604, 1.2088121152848075, 3.119917030457536};


scaff_module Y_Rx(qbit q) {  // Y in Table A1 of Whitfield et al 2010 really means Rx(-pi/2)
  Rx(q, -PI/2);
}


scaff_module Y_Rxdagger(qbit q) {
  Rx(q, PI/2);
}


scaff_module CNOTLadder(qbit reg[N], int controlStartIndex, int controlStopIndex) {
  // Applies a ladder of CNOTs, as in the dashed-CNOT notation at bottom of
  // Table A1 of Whitfield et al 2010
  int delta = controlStopIndex > controlStartIndex ? 1 : -1;
  for (int index = controlStartIndex;; index += delta) {
    CNOT(reg[index], reg[index - 1]);
    if (index == controlStopIndex)
      break;
  }
}


scaff_module DoubleExcitationOperator(qbit reg[N], int p, int q, int r, int s) {
  // Prerequisite: p > q > r > s
  // See Double Excitation Operator circuit in Table A1 of Whitfield et al 2010
  for (int i = 1; i <= 8; i++) {
    // NB(pranav): for some reason, ScaffCC wouldn't let me move the M operator
    // into a separate scaff_module. Seems like a ScaffCC bug.
    switch (i) {  // M operator (as defined in A1)
      case 1: H(reg[p]); H(reg[q]); H(reg[r]); H(reg[s]); break;
      case 2: Y_Rx(reg[p]); Y_Rx(reg[q]); Y_Rx(reg[r]); Y_Rx(reg[s]); break;
      case 3: H(reg[p]); Y_Rx(reg[q]); H(reg[r]); Y_Rx(reg[s]); break;
      case 4: Y_Rx(reg[p]); H(reg[q]); Y_Rx(reg[r]); H(reg[s]); break;
      case 5: Y_Rx(reg[p]); Y_Rx(reg[q]); H(reg[r]); H(reg[s]); break;
      case 6: H(reg[p]); H(reg[q]); Y_Rx(reg[r]); Y_Rx(reg[s]); break;
      case 7: Y_Rx(reg[p]); H(reg[q]); H(reg[r]); Y_Rx(reg[s]); break;
      case 8: H(reg[p]); Y_Rx(reg[q]); Y_Rx(reg[r]); H(reg[s]); break;
    }
    CNOTLadder(reg, p, q+1); CNOT(reg[q], reg[r]); CNOTLadder(reg, r, s+1);
    Rz(reg[s], Theta_p_q_r_s[p][q][r][s]);
    CNOTLadder(reg, s+1, r); CNOT(reg[q], reg[r]); CNOTLadder(reg, q+1, p);
    switch (i) {  // M dagger
      case 1: H(reg[p]); H(reg[q]); H(reg[r]); H(reg[s]); break;
      case 2: Y_Rxdagger(reg[p]); Y_Rxdagger(reg[q]); Y_Rxdagger(reg[r]); Y_Rxdagger(reg[s]); break;
      case 3: H(reg[p]); Y_Rxdagger(reg[q]); H(reg[r]); Y_Rxdagger(reg[s]); break;
      case 4: Y_Rxdagger(reg[p]); H(reg[q]); Y_Rxdagger(reg[r]); H(reg[s]); break;
      case 5: Y_Rxdagger(reg[p]); Y_Rxdagger(reg[q]); H(reg[r]); H(reg[s]); break;
      case 6: H(reg[p]); H(reg[q]); Y_Rxdagger(reg[r]); Y_Rxdagger(reg[s]); break;
      case 7: Y_Rxdagger(reg[p]); H(reg[q]); H(reg[r]); Y_Rxdagger(reg[s]); break;
      case 8: H(reg[p]); Y_Rxdagger(reg[q]); Y_Rxdagger(reg[r]); H(reg[s]); break;
    }
  }
}


scaff_module SingleExcitationOperator(qbit reg[N], int p, int q) {
  // Prerequisite: p > q
  // See Single Excitation Operator circuit in Table A1 of Whitfield et al 2010
  H(reg[p]); H(reg[q]);
  CNOTLadder(reg, p, q+1);
  Rz(reg[q], Theta_p_q[p][q]);
  CNOTLadder(reg, q+1, p);
  H(reg[p]); H(reg[q]);

  Y_Rx(reg[p]); Y_Rx(reg[q]);
  CNOTLadder(reg, p, q+1);
  Rz(reg[q], Theta_p_q[p][q]);
  CNOTLadder(reg, q+1, p);
  Y_Rx(reg[p]); Y_Rx(reg[q]);
}


int main() {
  qbit reg[N];

  // enumerate all N > p > q > r > s >= 0 and apply Double Excitation Operator
  for (int p = 0; p < N; p++) {
    for (int q = 0; q < p; q++) {
      for (int r = 0; r < q; r++) {
        for (int s = 0; s < r; s++) {
          DoubleExcitationOperator(reg, p, q, r, s);
        }
      }
    }
  }

  // enumerate all N > p > q >= 0 and apply Single Excitation Operator
  for (int p = 0; p < N; p++) {
    for (int q = 0; q < p; q++) {
      SingleExcitationOperator(reg, p, q);
    }
  }

  return 0;
}