// quantum computing emulation: Quantum Teleportation
//
// Usage: QT [r]
//
//  r = number of entangled qubits, default = 1
//
// Only one random qubit will be transported,
// but it can be entangled with r-1 other qubits.
//
// R. Perry, May 2019
//
#include <stdio.h>
#include <stdlib.h>
#include <tgmath.h>
#include <float.h>
#include "qce.h"

int main( int argc, char *argv[])
{
  State A = Q(2);  int b = PRINT_BITS;  unsigned int seed = SRAND(), r = 1;

 // number of entangled qubits
 //
  if( argc > 1) r = atoi(argv[1]);

 // random qubit to transport: R = alpha|0> + beta|1>
 //
  double complex alpha = DRAND*exp(CMPLX(0,pi*(2*DRAND-1)));

  double complex beta = sqrt(1-abs2(alpha))*exp(CMPLX(0,pi*(2*DRAND-1)));

  State R = Q(1);  R.a[0] = alpha;  R.a[1] = beta;

  printf( "QT: r = %u, seed = %u, alpha = (%g,%g), beta = (%g,%g)\n\n",
    r, seed, creal(alpha), cimag(alpha), creal(beta), cimag(beta));

  if( r > 1)
  {
    State T = Q(r-1);  init( T, 0, 0);  R = product( T, R);

   // apply H to qubits 1...(R.n-1)
   //
    for( unsigned int k = 1; k < R.n; ++k) H( R, k);

   // apply CX to qubit 0 controlled by qubits 1...(R.n-1)
   //
    for( unsigned int c = 1; c < R.n; ++c) CX( R, c, 0);
  } 

 // for each of Alice's possible measurements of q2,q1
 //
  for( unsigned int m = 0; m < 4; ++m)
  {
   // initial state |00>
   //
    init( A, 0, 0);  if( m == 0) print( "A", A, b); 

   // Alice and Bob jointly entangle A1,A0
   //
    H( A, 1);  if( m == 0) print( "H1", A, b);

    CX( A, 1, 0);  if( m == 0) print( "CX10", A, b);

   // q = combined state of R and A: |...,q2,q1,q0> = |R>|A1,A0>
   //
    State q = product( R, A);

    if( m == 0) { print( "R", R, b); print( "Q", q, b); }

   // Alice entangles q2 with q1
   //
    CX( q, 2, 1);  if( m == 0) print( "CX21", q, b);

    H( q, 2);  if( m == 0) print( "H2", q, b);

   // Alice measures, collapse state
   //
    collapse( q, 2, m >> 1);  collapse( q, 1, m & 1);

    printb( m, 2);  print( "", q, b);

   // Bob transformations
   //
    if( m & 1) { X( q, 0); print( "X", q, b); }

    if( m >> 1) { Z( q, 0); print( "Z", q, b); }

    if( r > 1) // show reduced state, should equal R
    {
      State T = Q(q.n-2);  double s = 0;

     //    ...,q3,q2,q1,q0 qubits 2,1 are eliminated by the measurement
     // ...,T2,T1,m2,m1,R0 -> ...,T2,T1,R0
     //
      for( unsigned int i = 0, j = 0; i < T.N; ++i)
      {
        j = ((i>>1)<<3) | (m <<1) | (i&1);  T.a[i] = q.a[j];  s += abs2(q.a[j]);
      } 

     // normalize and compare
     //
      s = sqrt(s);  for( unsigned int i = 0; i < T.N; ++i) T.a[i] /= s;

      print( "T", T, b);  print( "R", R, b);

      double diff = 0; for( unsigned int i = 0; i < T.N; ++i) diff += abs2(T.a[i]-R.a[i]);

     // threshold for zero norm squared: 2 * sum { |e + i*e|**2 } = 2 * (2*e*e*(T.N/2))
     //
      if( diff > 2*DBL_EPSILON*DBL_EPSILON*T.N) printf( "T != R, diff = %g\n", diff);
    }
  }

  if( r < 2) printf( "\nalpha = (%g,%g), beta = (%g,%g)\n",
    creal(alpha), cimag(alpha), creal(beta), cimag(beta));

  return 0;
}