// simulation using one quantum bit
//
// constraint: sum |a[i]|**2 = 1 ==>> unit circle
//
// |a[0]|**2 = Prob(0)
// |a[1]|**2 = Prob(1)
//
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <tgmath.h>

typedef double complex Complex;

// swap a[i] and a[j]
void swap( Complex a[], int i, int j) { Complex t = a[i]; a[i] = a[j]; a[j] = t; }

void X( Complex a[]) { swap( a, 0, 1); } // NOT

void Z( Complex a[]) { a[1] = -a[1]; } // phase flip

void S( Complex a[]) { a[1] = CMPLX(0,1)*a[1]; } // sqrt(Z)

void H( Complex a[]) { double s2 = 1/sqrt(2); // Hadamard
  Complex t0 = a[0], t1 = a[1]; a[0] = s2*(t0+t1); a[1] = s2*(t0-t1);
}
// OpenQASM3 U with global phase gamma // untested
void U( Complex a[], double theta, double phi, double lambda, double gamma) {
  double t = theta/2; Complex t0 = a[0], t1 = a[1], p = exp(CMPLX(0,gamma)),
  c = p*cos(t), s = p*sin(t), u00 = c, u01 = -exp(CMPLX(0,lambda))*s,
  u10 = exp(CMPLX(0,phi))*s, u11 = exp(CMPLX(0,phi+lambda))*c;
  a[0] = u00*t0+u01*t1; a[1] = u10*t0+u11*t1;
}
// magnitude squared
// Equivalent to abs(x)**2 but without the square-root and square operations
// In C++ this is available as norm(x) for complex
double abs2( Complex c) { double x = creal(c), y = cimag(c); return x*x + y*y; }
//
//  0.0  ... r ...       1.0
//   |<------->|<-------->|
//       p0         p1
// 
void M( Complex a[]) { // measure
  double r = rand() / (RAND_MAX + 1.0); // 0.0 <= r < 1.0
  double p = abs2( a[0]);
  if( r < p) { printf( "M = 0\n"); a[0] = 1; a[1] = 0; }
  else { printf( "M = 1\n"); a[0] = 0; a[1] = 1; }
}
void print( const char msg[], const Complex a[]) {
  printf( "%s:", msg);
  for( int i = 0; i < 2; ++i) printf( " (%g,%g)", creal(a[i]), cimag(a[i]));
  printf( "\n");
}
int main( void) {
  srand(time(0)); Complex a[2] = { 3.0/5.0, CMPLX(0,4.0/5.0) }; // amplitudes
  print( "a", a);
  X(a); print( "X", a);
  X(a); print( "X", a);
  H(a); print( "H", a);
  Z(a); print( "Z", a);
  S(a); print( "S", a);
  S(a); print( "S", a);
  H(a); print( "H", a);
  M(a); print( "M", a);
}

/* output:

a: (0.6,0) (0,0.8)
X: (0,0.8) (0.6,0)
X: (0.6,0) (0,0.8)
H: (0.424264,0.565685) (0.424264,-0.565685)
Z: (0.424264,0.565685) (-0.424264,0.565685)
S: (0.424264,0.565685) (-0.565685,-0.424264)
S: (0.424264,0.565685) (0.424264,-0.565685)
H: (0.6,0) (0,0.8)

*/