// quantum computing emulation: Grover's search
//
// Usage: Grover [-erpdE] [qn [k [p [m]]]]
//
//  qn = number of qubits
//   k = number of iterations
//   p = probability of depolarizing noise, or non-resonant pulse error level
//   m = value to match
//  -e  non-resonant pulse error
//  -r  use random global phase factor
//  -p  show state probabilities
//  -d  extra debug output
//  -E  show entanglement
//
//  if m is negative or not specified then it is generated randomly
//
// R. Perry, April 2018
//
#include <stdio.h>
#include <stdlib.h> // atoi(), rand()
#include <tgmath.h>
#include "qce.h"

// print probabilities
//
void print_p( const char *msg, State q)
{
  printf( "%s:\n", msg);

  for( unsigned long i = 0; i < q.N; ++i) printf( " %lu %g\n", i, abs2(q.a[i]));
}

//---------------------------------------------------------------------

int main( int argc, char *argv[])
{
 // random phase and debug options
 //
  int rphase = 0, debug = 0, eflag = 0, pflag = 0, Eflag = 0;

  unsigned int seed = SRAND();

  if( argc > 1 && argv[1][0] == '-') 
  {
    for( int i = 1; argv[1][i]; ++i) switch( argv[1][i])
    {
      case 'e': ++eflag; break;
      case 'r': ++rphase; break;
      case 'p': ++pflag; break;
      case 'd': ++debug; break;
      case 'E': ++Eflag; break;
      default: error( "Grover: bad option");
    }
    --argc, ++argv;
  }

 // number of qubits
 //
  unsigned int qn = 2;  if( argc > 1) qn = atoi(argv[1]);

  if( qn == 0) error( "Grover: qn must be non-zero");

 // number of iterations
 //
  int k = 2;  if( argc > 2) k = atoi(argv[2]);

 // probability of depolarizing noise, or non-resonant pulse error level
 //
  double p = 0, ra = 0;  if( argc > 3) ra = atof(argv[3]); // ra can be negative

  if( !eflag) { p = ra; ra = 0; }

  if( p < 0 || p > 1) error( "Grover: bad p arg");

 // value to match
 //
  State q = Q(qn); unsigned long m; // f(m)=1, otherwise 0

  long arg = -1;  if( argc > 4) arg = atol(argv[4]);

  if( arg < 0) m = IRAND(q.N); else m = arg;

  if( m >= q.N) error( "Grover: bad m arg");

  printf( "G: seed = %u, qn = %u, k = %i, p = %g, ra = %g, m = %lu\n", seed, qn, k, p, ra, m);

  init( q, 0, rphase); // initial state |0...0>

  if( debug) { print( "q", q, 1); if( Eflag) print( "E", q, -1); }

// note that the auxiliary qubit is not included here, it is not needed
// since we will just directly flip the sign of state m in the iterations

  for( unsigned int i = 0; i < q.n; ++i)
  {
    Hra( q, i, ra);  if( debug) { print( "H", q, 1); if( Eflag) print( "E", q, -1); }
  }

  double h2 = 1 / (double) q.N; // for depolarizing noise

  if( pflag) print_p( "P", q);

  for( int iter = 1; iter <= k; ++iter)
  {
    if( debug) printf( "iter=%i\n", iter);

   // apply f, i.e. flip sign of state m
   //
    q.a[m] = -q.a[m];  if( debug) { print( "f", q, 1); if( Eflag) print( "E", q, -1); }

   // inversion about the average
   //
    double complex avg = 0; for( unsigned long i = 0; i < q.N; ++i) avg += q.a[i];

    avg *= 2.0/q.N; for( unsigned long i = 0; i < q.N; ++i) q.a[i] = avg - q.a[i];

    // printf( "avg = (%g,%g)\n", creal(avg), cimag(avg));

    if( debug) { print( "inv", q, 1); if( Eflag) print( "E", q, -1); }

    if( debug) // check max probability
    {
      unsigned long s = 0;  double p = abs2(q.a[0]), t;
      for( unsigned long i = 1; i < q.N; ++i) { t = abs2(q.a[i]); if( t > p) { p = t; s = i; } }
      if( s != m) printf( "G: s != m\n");
      printf( "P(%lu) = %g\n", s, p);
    }
    else
      if( pflag) print_p( "P", q);
    else
      printf( "%d %g\n", iter, (1-p)*abs2(q.a[m]) + p*h2);
  }

  return 0;
}