// quantum computing emulation: Quantum Walk in 1 dimension, on a ring
//
// Usage: QW1 [-rdE] [c [p [T]]]
//
//   c = coin initialization: 0 (|0>, default), 1 (|1>), 2 (|0>+i|1>),
//        3 (sqrt(0.85)|0>-sqrt(0.15)|1>)
//   p = number of position qubits, default = 3, P = 2**p = 8 positions
//   T = number of time steps, default = min(P,10)
//  -r  use random global phase factor
//  -d  extra debug output, -dd for even more output
//  -E  show entanglement
//
// References:
//
//  J Kempe, Quantum random walks: An introductory overview,
//  Contemporary Physics, Volume 44, 2003, Issue 4 
//  https://doi.org/10.1080/00107151031000110776 
//
//  Salvador Elías Venegas-Andraca, Quantum walks: a comprehensive review
//  Quantum Information Processing, October 2012, Volume 11, Issue 5, pp 1015–1106 
//  https://doi.org/10.1007/s11128-012-0432-5
//
// R. Perry, August 2019
//
#include <stdio.h>
#include <stdlib.h> // atoi()
#include <tgmath.h>
#include "qce.h"

// display probabilities with position offset, 0 is center of the ring
//
// state = |position,coin>,  prob(position) = prob(|position,0>) + prob(|position,1>)
//
void disp( State q, unsigned long t, int debug, int E)
{
  long offset = q.N>>2;  if( debug) printf( "t=%lu:\n", t);

  for( unsigned long i = 0; i < q.N; i += 2)
    printf( "%6li %g\n", (long)(i>>1)-offset, abs2(q.a[i]) + abs2(q.a[i+1]) );

  if( debug > 1) { print( "q", q, PRINT_BITS); if( E) print( "E", q, PRINT_TANGLE); }
}

// position update, mod P
//
// move right: S |i,0> = |i+1,0>
//
// move left:  S |i,1> = |i-1,1>
//
void step( State q)
{
  double complex R = q.a[q.N-2], L = q.a[1], t;

  #define SWAP(a,b) t=a; a=b; b=t;

  for( unsigned long i = 0, k = q.N-1; i < q.N; i += 2, k -= 2)
  {
    SWAP( q.a[i], R);  SWAP( q.a[k], L);
  }
}

int main( int argc, char *argv[])
{
  int rphase = 0, debug = 0, E = 0;  unsigned int seed = SRAND();

 // command-line options
 //
  if( argc > 1 && argv[1][0] == '-') 
  {
    for( int i = 1; argv[1][i]; ++i) switch( argv[1][i])
    {
      case 'r': ++rphase; break;
      case 'd': ++debug; break;
      case 'E': ++E; break;
      default: error( "QW1: bad option");
    }
    --argc, ++argv;
  }

 // coin initialization: 0 (|0>, default), 1 (|1>), or 2 (|0>+i|1>)
 //
  int coin = 0; if( argc > 1) { coin = atoi(argv[1]);

    if( coin < 0 || coin > 3) error( "QW1: bad coin arg"); }

 // number of position qubits
 //
  unsigned int p = 3; if( argc > 2) { p = atoi(argv[2]);

    if( p < 1) error( "QW1: bad p arg"); }

 // number of time steps
 //
  unsigned long T = 1 << p;  if( T > 10) T = 10;

  if( argc > 3) { T = atoi(argv[3]); if( T < 1) error( "QW1: bad T arg"); }

  printf( "QW1: seed = %u, coin = %i, p = %u, T = %lu\n", seed, coin, p, T);

 // state = |position,coin>, initial position = P/2
 //
  State q = Q(p+1);  init( q, 0, rphase);  unsigned long mid = q.N>>1;

  switch( coin) // move initial values to middle of the ring
  {
    case 0: q.a[mid] = q.a[0]; break;
    case 1: q.a[mid+1] = q.a[0]; break;
    case 2: q.a[mid] = s2*q.a[0]; q.a[mid+1] = I*q.a[mid]; break;
    case 3: q.a[mid] = sqrt(0.85)*q.a[0]; q.a[mid+1] = -sqrt(0.15)*q.a[0]; break;
  }
  q.a[0] = 0;  if( debug) disp( q, 0, debug, E);

 // run for T time steps
 //
  for( unsigned long t = 1; t <= T; ++t)
  {
    H( q, 0); step( q); // Hadamard coin, position update

    if( debug) disp( q, t, debug, E);
  }

 // final state
 //
  if( !debug) disp( q, T, 0, 0);

  return 0;
}