// qgates.h  #include after #include "qfunct.h"
//           arrays built dynamically, call  gate_init()  only once
//           in C++  c_vec keth(2);                      |    |
//                   keth = matvec(hada, ket);  |keth| = |hada| * |ket|
//                                                       |    |

static c_mat hada;   // Hadamard matrix built by  hada_init();
static c_mat paulix; // equivalent to not
static c_mat pauliy; // equivalent to rotate about Y axis
static c_mat pauliz; // equivalent to rotate about Z axis
static c_mat cnot;   // controlled not
static c_mat ccnot;  // two controlled not, Toffoli
static c_mat swap2;  // swap2
static c_mat cswap;  // controlled swap
static c_mat Rk_make(double k); // rotation e^(2 Pi i / 2^k)
static c_mat cRk_make(double k); // controlled rotation e^(2 Pi i / 2^k)
static c_mat rsx;    // root sigma x
static c_mat rsxt;   // root sigma x conjugate transpose
static c_vec ccnot3(c_vec keta, c_vec ketb, c_vec ketc); 
static c_vec cnot2(c_vec keta, c_vec ketb); 
static c_vec rsx2(c_vec keta, c_vec ketb); 
static c_vec rsxt2(c_vec keta, c_vec ketb); 

static c_vec ket0(2);
static c_vec ket1(2);
static c_vec bra0(2);
static c_vec bra1(2);

static void hada_init(void) // hadamard
{
  c_vec row(2);
  row[0] = cmplx(1.0/sqrt(2.0),0);
  row[1] = cmplx(1.0/sqrt(2.0),0);
  hada.push_back(row);
  row[0] = cmplx(1.0/sqrt(2.0),0);
  row[1] = cmplx(-1.0/sqrt(2.0),0);
  hada.push_back(row);
} // end hada_init  call only once

static void paulix_init(void) // single qbit not
{
  c_vec row(2);
  row[0] = cmplx(0,0);
  row[1] = cmplx(1.0,0);
  paulix.push_back(row);
  row[0] = cmplx(1.0,0);
  row[1] = cmplx(0,0);
  paulix.push_back(row);
} // end hada_init  call only once

static void pauliy_init(void) // single qbit rotate
{
  c_vec row(2);
  row[0] = cmplx(0,0);
  row[1] = cmplx(0,-1.0);
  pauliy.push_back(row);
  row[0] = cmplx(0,1.0);
  row[1] = cmplx(0,0);
  pauliy.push_back(row);
} // end pauliy_init  call only once

static void pauliz_init(void) // single qbit rotate
{
  c_vec row(2);
  row[0] = cmplx(1.0,0);
  row[1] = cmplx(0,0);
  pauliz.push_back(row);
  row[0] = cmplx(0,0);
  row[1] = cmplx(0,-1.0);
  pauliz.push_back(row);
} // end pauliz_init  call only once

static void cnot_init(void)
{
  c_vec row(4);
  row[0] = cmplx(1,0);
  row[1] = cmplx(0,0);
  row[2] = cmplx(0,0);
  row[3] = cmplx(0,0);
  cnot.push_back(row);
  row[0] = cmplx(0,0);
  row[1] = cmplx(1,0);
  cnot.push_back(row);
  row[1] = cmplx(0,0);
  row[3] = cmplx(1,0);
  cnot.push_back(row);
  row[3] = cmplx(0,0);
  row[2] = cmplx(1,0);
  cnot.push_back(row);
} // end cnot_init  call only once

static void ccnot_init(void) // Tofffoli
{
  c_vec row(8);
  for(int i=0; i<8; i++) row[i] = cmplx(0,0);
  row[0] = cmplx(1,0);
  ccnot.push_back(row);
  row[0] = cmplx(0,0);
  row[1] = cmplx(1,0);
  ccnot.push_back(row);
  row[1] = cmplx(0,0);
  row[2] = cmplx(1,0);
  ccnot.push_back(row);
  row[2] = cmplx(0,0);
  row[3] = cmplx(1,0);
  ccnot.push_back(row);
  row[3] = cmplx(0,0);
  row[4] = cmplx(1,0);
  ccnot.push_back(row);
  row[4] = cmplx(0,0);
  row[5] = cmplx(1,0);
  ccnot.push_back(row);
  row[5] = cmplx(0,0);
  row[7] = cmplx(1,0);
  ccnot.push_back(row);
  row[7] = cmplx(0,0);
  row[6] = cmplx(1,0);
  ccnot.push_back(row);
} // end ccnot_init  call only once

static void swap2_init(void)
{
  c_vec row(4);
  for(int i=0; i<4; i++) row[i] = cmplx(0,0);
  row[0] = cmplx(1,0);
  swap2.push_back(row);
  row[0] = cmplx(0,0);
  row[2] = cmplx(1,0);
  swap2.push_back(row);
  row[2] = cmplx(0,0);
  row[1] = cmplx(1,0);
  swap2.push_back(row);
  row[1] = cmplx(0,0);
  row[3] = cmplx(1,0);
  swap2.push_back(row);
} // end swap2_init  call only once

static void cswap_init(void)
{
  c_vec row(8);
  for(int i=0; i<8; i++) row[i] = cmplx(0,0);
  row[0] = cmplx(1,0);
  cswap.push_back(row);
  row[0] = cmplx(0,0);
  row[1] = cmplx(1,0);
  cswap.push_back(row);
  row[1] = cmplx(0,0);
  row[2] = cmplx(1,0);
  cswap.push_back(row);
  row[2] = cmplx(0,0);
  row[3] = cmplx(1,0);
  cswap.push_back(row);
  row[3] = cmplx(0,0);
  row[4] = cmplx(1,0);
  cswap.push_back(row);
  row[4] = cmplx(0,0);
  row[6] = cmplx(1,0);
  cswap.push_back(row);
  row[6] = cmplx(0,0);
  row[5] = cmplx(1,0);
  cswap.push_back(row);
  row[5] = cmplx(0,0);
  row[7] = cmplx(1,0);
  cswap.push_back(row);
} // end cswap_init  call only once

static c_mat Rk_make(double k) // rotation e^(2 Pi i / 2^k)
{
  c_mat Rk;
  c_vec row(2);
  double Pi = 3.1415926535897932384626433832795028841971;
  double p;
  cmplx x;
  p = pow(2.0,k);
  x = cmplx(0, 2.0*Pi/p);
  cout << "p=" << p << ", x=" << x << ", exp(x)=" << exp(x) << endl;
  row[0] = cmplx(1,0);
  row[1] = cmplx(0,0);
  Rk.push_back(row);
  row[0] = cmplx(0,0);
  row[1] = exp(x);
  Rk.push_back(row);
  return Rk;
} // end Rk_make

static c_mat cRk_make(double k) // controlled rotation e^(2 Pi i / 2^k)
{
  c_mat cRk;
  c_vec row(4);
  double Pi = 3.1415926535897932384626433832795028841971;
  double p;
  cmplx x;
  p = pow(2.0,k);
  x = cmplx(0, 2.0*Pi/p);
  for(int i=0; i<4; i++) row[i] = cmplx(0,0);
  row[0] = cmplx(1,0);
  cRk.push_back(row);
  row[0] = cmplx(0,0);
  row[1] = cmplx(1,0);
  cRk.push_back(row);
  row[1] = cmplx(0,0);
  row[2] = cmplx(1,0);
  cRk.push_back(row);
  row[2] = cmplx(0,0);
  row[3] = exp(x);
  cRk.push_back(row);
  return cRk;
} // end cRk_make

static void rsx_init(void)   // root sigma x
{
  c_vec row(4);
  row[0] = cmplx(1,0);
  row[1] = cmplx(0,0);
  row[2] = cmplx(0,0);
  row[3] = cmplx(0,0);
  rsx.push_back(row);
  row[0] = cmplx(0,0);
  row[1] = cmplx(1,0);
  rsx.push_back(row);
  row[1] = cmplx(0,0);
  row[2] = cmplx(1.0/2.0, 1.0/2.0);  // 1/2 1+i
  row[3] = cmplx(1.0/2.0, -1.0/2.0); // 1/2 1-i
  rsx.push_back(row);
  row[2] = cmplx(1.0/2.0, -1.0/2.0); // 1/2 1-i
  row[3] = cmplx(1.0/2.0, 1.0/2.0);  // 1/2 1+i
  rsx.push_back(row);
} // end rsx_init  call only once

static void rsxt_init(void)   // root sigma x complex transpose
{
  c_vec row(4);
  row[0] = cmplx(1,0);
  row[1] = cmplx(0,0);
  row[2] = cmplx(0,0);
  row[3] = cmplx(0,0);
  rsxt.push_back(row);
  row[0] = cmplx(0,0);
  row[1] = cmplx(1,0);
  rsxt.push_back(row);
  row[1] = cmplx(0,0);
  row[2] = cmplx(1.0/2.0, -1.0/2.0); // 1/2 1-i
  row[3] = cmplx(1.0/2.0, 1.0/2.0);  // 1/2 1+i
  rsxt.push_back(row);
  row[2] = cmplx(1.0/2.0, 1.0/2.0);  // 1/2 1+i
  row[3] = cmplx(1.0/2.0, -1.0/2.0); // 1/2 1-i
  rsxt.push_back(row);
} // end rsxt_init  call only once

static c_vec cnot2(c_vec keta, c_vec ketb)
{
  c_vec ketab(4);
  c_vec gout(4);
  c_vec ketout(2);
  // cnot_init(); needed previously
  ketab = outp(keta, ketb);
  gout = matvec(cnot, ketab);
  ketout = factor(gout, keta);
  return ketout;
} // end cnot2  two input kets, one output ket, first not changed.

static c_vec ccnot3(c_vec keta, c_vec ketb, c_vec ketc)
{
  c_vec ketabc(8);
  c_vec gout(8);
  c_vec ketout(2);
  ketabc = outp(outp(keta, ketb), ketc);
  gout = matvec(ccnot, ketabc);
  ketout = factor3(gout, keta, ketb);
  return ketout;
} // end ccnot3  three input kets, just last may be changed.

static c_vec rsx2(c_vec keta, c_vec ketb)
{
  c_vec ketab(4);
  c_vec gout(4);
  c_vec ketout(2);
  ketab = outp(keta, ketb);
  gout = matvec(rsx, ketab);
  // prtvec(gout);
  ketout = factor(gout, keta);
  return ketout;
} // end rsx2
 
static c_vec rsxt2(c_vec keta, c_vec ketb)
{
  c_vec ketab(4);
  c_vec gout(4);
  c_vec ketout(2);
  ketab = outp(keta, ketb);
  gout = matvec(rsxt, ketab);
  // prtvec(gout);
  ketout = factor(gout, keta);
  return ketout;
} // end rsxt2

static void gate_init()
{

  ket0[0] = cmplx(1.0,0.0);
  ket0[1] = cmplx(0.0,0.0);
  ket1[0] = cmplx(0.0,0.0);
  ket1[1] = cmplx(1.0,0.0);
  bra0[0] = cmplx(1.0,0.0);
  bra0[1] = cmplx(0.0,0.0);
  bra1[0] = cmplx(0.0,0.0);
  bra1[1] = cmplx(1.0,0.0);

  hada_init();   // Hadamard matrix built by  hada_init();
  paulix_init(); // equivalent to not
  pauliy_init(); // equivalent to rotate about Y axis
  pauliz_init(); // equivalent to rotate about Z axis
  cnot_init();   // controlled not
  ccnot_init();  // two controlled not, Toffoli
  swap2_init();  // swap2
  cswap_init();  // controlled swap
  rsx_init();    // root sigma x
  rsxt_init();   // root sigma x conjugate transpose
} // end gate_init