Update 31 May 2019: added section at end on Entangled Qubits.

For background on quantum teleportation see Quantum Computation and Quantum Information, 10th Anniversary Edition, Michael A. Nielsen & Isaac L. Chuang, Cambridge University Press, 2010, pages 26-28.

Annotated code and output:

Three qubits are used. Initially Alice and Bob share qubits A1 and A0, and R is the unknown/random qubit in state α|0>+β|1> which will be transported:

  State A = Q(2);  int b = PRINT_BITS;  double pi = 4*atan(1);  SRAND();

 // 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;

In this example the randomly generated values of α and β are α = (-0.629227,0.197361) and β = (0.750464,-0.0438817).

Each possible measurement by Alice will be simulated separately. In each case the initial state is |A1,A0> = |00>:

 // 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); 

A:
  00 1        (1,0)
  01 0        (0,0)
  10 0        (0,0)
  11 0        (0,0)

The output shows the state indices in binary, the probability of each state (magnitude-squared of the state amplitude), and the complex state amplitudes (real,imag).

Now qubits A1 and A0 are entangled by performing a Hadamard transformation on qubit A1, followed by a controlled-NOT on qubit A0 controlled by qubit A1:

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

H1:
  00 0.5      (0.707107,0)
  01 0        (0,0)
  10 0.5      (0.707107,0)
  11 0        (0,0)

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

CX10:
  00 0.5      (0.707107,0)
  01 0        (0,0)
  10 0        (0,0)
  11 0.5      (0.707107,0)

Now Bob takes qubit A0 and travels to Mars.

At some point in the future Alice obtains qubit R which is in an unknown state α|0>+β|1>:

   // 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); }

R:
  0 0.4349   (-0.629227,0.197361)
  1 0.5651   (0.750464,-0.0438817)
Q:
  000 0.2174   (-0.444931,0.139555)
  001 0        (0,0)
  010 0        (0,0)
  011 0.2174   (-0.444931,0.139555)
  100 0.2826   (0.530658,-0.0310291)
  101 0        (0,0)
  110 0        (0,0)
  111 0.2826   (0.530658,-0.0310291)

Alice entangles R (q2) with A1 (q1) by performing a controlled-NOT on q1 controlled by q2, followed by a Hadamard transformation on q2:

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

CX21:
  000 0.2174   (-0.444931,0.139555)
  001 0        (0,0)
  010 0        (0,0)
  011 0.2174   (-0.444931,0.139555)
  100 0        (0,0)
  101 0.2826   (0.530658,-0.0310291)
  110 0.2826   (0.530658,-0.0310291)
  111 0        (0,0)

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

H2:
  000 0.1087   (-0.314614,0.0986805)
  001 0.1413   (0.375232,-0.0219409)
  010 0.1413   (0.375232,-0.0219409)
  011 0.1087   (-0.314614,0.0986805)
  100 0.1087   (-0.314614,0.0986805)
  101 0.1413   (-0.375232,0.0219409)
  110 0.1413   (-0.375232,0.0219409)
  111 0.1087   (-0.314614,0.0986805)

Now Alice measures q2 and q1, which collapses the state:

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

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

Alice sends her measurement (2 classical bits) to Bob.

If her measurement is 00, the state is:

00:
  000 0.4349   (-0.629227,0.197361)
  001 0.5651   (0.750464,-0.0438817)
  010 0        (0,0)
  011 0        (0,0)
  100 0        (0,0)
  101 0        (0,0)
  110 0        (0,0)
  111 0        (0,0)

In this case Bob does not have to perform any transformations. The state is |00>(α|0>+β|1>), i.e. the state of Bob's qubit q0 is equal to the state of the original random qubit. Recall that α = (-0.629227,0.197361) and β = (0.750464,-0.0438817) in this example.

In general Bob performs transformations on q0 depending on the bits received:

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

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

For 01 Bob performs a NOT transformation on q0:

01:
  000 0        (0,0)
  001 0        (0,0)
  010 0.5651   (0.750464,-0.0438817)
  011 0.4349   (-0.629227,0.197361)
  100 0        (0,0)
  101 0        (0,0)
  110 0        (0,0)
  111 0        (0,0)
X:
  000 0        (0,0)
  001 0        (0,0)
  010 0.4349   (-0.629227,0.197361)
  011 0.5651   (0.750464,-0.0438817)
  100 0        (0,0)
  101 0        (0,0)
  110 0        (0,0)
  111 0        (0,0)

The resulting state is |01>(α|0>+β|1>), i.e. again the state of Bob's qubit q0 is equal to the state of the original random qubit.

For 10 Bob performs a phase-flip on q0:

10:
  000 0        (0,0)
  001 0        (0,0)
  010 0        (0,0)
  011 0        (0,0)
  100 0.4349   (-0.629227,0.197361)
  101 0.5651   (-0.750464,0.0438817)
  110 0        (0,0)
  111 0        (0,0)
Z:
  000 0        (0,0)
  001 0        (0,0)
  010 0        (0,0)
  011 0        (0,0)
  100 0.4349   (-0.629227,0.197361)
  101 0.5651   (0.750464,-0.0438817)
  110 0        (0,0)
  111 0        (0,0)

The resulting state is |10>(α|0>+β|1>).

For 11 Bob performs a NOT transformation on q0, followed by a phase-flip on q0:

11:
  000 0        (0,0)
  001 0        (0,0)
  010 0        (0,0)
  011 0        (0,0)
  100 0        (0,0)
  101 0        (0,0)
  110 0.5651   (-0.750464,0.0438817)
  111 0.4349   (-0.629227,0.197361)
X:
  000 0        (0,0)
  001 0        (0,0)
  010 0        (0,0)
  011 0        (0,0)
  100 0        (0,0)
  101 0        (0,0)
  110 0.4349   (-0.629227,0.197361)
  111 0.5651   (-0.750464,0.0438817)
Z:
  000 0        (0,0)
  001 0        (0,0)
  010 0        (0,0)
  011 0        (0,0)
  100 0        (0,0)
  101 0        (0,0)
  110 0.4349   (-0.629227,0.197361)
  111 0.5651   (0.750464,-0.0438817)

The resulting state is |11>(α|0>+β|1>).

In each case, the state of Bob's qubit q0 is equal to the state of the original random qubit. Recall that α = (-0.629227,0.197361) and β = (0.750464,-0.0438817) in this example.

---

Entangled Qubits

The transported qubit inherits any entanglement it had with other qubits.

Example transporting random R0 from two entangled qubits R1,R0 with measurement 00 (QT2.out):

R:
  00 0.3834   (-0.113858,0.608675)
  01 0.1166   (-0.280054,-0.195247)
  10 0.1166   (-0.280054,-0.195247)
  11 0.3834   (-0.113858,0.608675)
00:
  0000 0.3834   (-0.113858,0.608675)
  0001 0.1166   (-0.280054,-0.195247)
  0010 0        (0,0)
  0011 0        (0,0)
  0100 0        (0,0)
  0101 0        (0,0)
  0110 0        (0,0)
  0111 0        (0,0)
  1000 0.1166   (-0.280054,-0.195247)
  1001 0.3834   (-0.113858,0.608675)
  1010 0        (0,0)
  1011 0        (0,0)
  1100 0        (0,0)
  1101 0        (0,0)
  1110 0        (0,0)
  1111 0        (0,0)

The state is |R1>|00>|R0>, which reduces to T=|R1>|R0>, which equals R:

T:
  00 0.3834   (-0.113858,0.608675)
  01 0.1166   (-0.280054,-0.195247)
  10 0.1166   (-0.280054,-0.195247)
  11 0.3834   (-0.113858,0.608675)

Example transporting random R0 from three entangled qubits R2,R1,R0 with measurement 11 (QT3.out):

R:
  000 0.1361   (0.125447,0.346873)
  001 0.1139   (-0.27759,0.192058)
  010 0.1139   (-0.27759,0.192058)
  011 0.1361   (0.125447,0.346873)
  100 0.1139   (-0.27759,0.192058)
  101 0.1361   (0.125447,0.346873)
  110 0.1361   (0.125447,0.346873)
  111 0.1139   (-0.27759,0.192058)
11: X: Z:
  00000 0        (0,0)
  00001 0        (0,0)
  00010 0        (0,0)
  00011 0        (0,0)
  00100 0        (0,0)
  00101 0        (0,0)
  00110 0.1361   (0.125447,0.346873)
  00111 0.1139   (-0.27759,0.192058)
  01000 0        (0,0)
  01001 0        (0,0)
  01010 0        (0,0)
  01011 0        (0,0)
  01100 0        (0,0)
  01101 0        (0,0)
  01110 0.1139   (-0.27759,0.192058)
  01111 0.1361   (0.125447,0.346873)
  10000 0        (0,0)
  10001 0        (0,0)
  10010 0        (0,0)
  10011 0        (0,0)
  10100 0        (0,0)
  10101 0        (0,0)
  10110 0.1139   (-0.27759,0.192058)
  10111 0.1361   (0.125447,0.346873)
  11000 0        (0,0)
  11001 0        (0,0)
  11010 0        (0,0)
  11011 0        (0,0)
  11100 0        (0,0)
  11101 0        (0,0)
  11110 0.1361   (0.125447,0.346873)
  11111 0.1139   (-0.27759,0.192058)

The state is |R2>|R1>|11>|R0>, which reduces to T=|R2>|R1>|R0>, which equals R:

T:
  000 0.1361   (0.125447,0.346873)
  001 0.1139   (-0.27759,0.192058)
  010 0.1139   (-0.27759,0.192058)
  011 0.1361   (0.125447,0.346873)
  100 0.1139   (-0.27759,0.192058)
  101 0.1361   (0.125447,0.346873)
  110 0.1361   (0.125447,0.346873)
  111 0.1139   (-0.27759,0.192058)