// quantum computing emulation: Simon's problem
//
// Determine the period of a function:
//
//   f(y) == f(x)  <<==>>  y = x XOR period
//
// Usage: Simon [-rdE] [qn [k]]
//
// qn must be even; qn/2 is the number of bits in the periodic function
//
// k is the bit number to set high, or negative to use a random periodic function
//
//  -r  use random global phase factor
//  -d  extra debug output
//  -E  show entanglement
//
// R. Perry, Jan. 2018
//
#include <stdio.h>
#include <stdlib.h> // rand(), free()
#include <tgmath.h>
#include "qce.h"

int debug = 0; // debug option

// printlu: print unsigned long array
//
void printlu( const char *msg, unsigned long R[], unsigned long N)
{
  printf( "%3s:", msg);

  for( unsigned long i = 0; i < N; ++i) { printf( " %lu:%lu", i, R[i]); } putchar( '\n');
}

static unsigned long *F; // function with random period
static unsigned long nF; // number of elements

// f_FR: classical b ^ F(a)
//
unsigned long f_FR( unsigned long a, unsigned long b, unsigned long unused)
{
  return b ^ F[a];
}

// FR: (a,b) -> (a,b+F(a)), random periodic F()
//
unsigned long FR( State q)
{
  unsigned int n = q.n >> 1;
  unsigned long N = 1LU << n;

  if( F && nF != N) { free(F); F = NULL; }

  if( !F) { F = emalloc( N*sizeof( unsigned long)); nF = N; }

  unsigned long t, m, i, *R = emalloc( N*sizeof(unsigned long));

  for( i = 0; i < N; ++i) { R[i] = i; F[i] = 0; }

 // R = 0 ... (N-1) in random order
 //
  while( i > 1) { m = IRAND(i); --i; t = R[i]; R[i] = R[m]; R[m] = t; }

  if( debug) printlu( "R", R, N);

  unsigned long X, Y, j, k, period = 1 + IRAND(N-1); // [1,...,N-1]

 // create function table
 //
  for( i = j = 0, k = 1; i < N; i += 2, ++k)
  {
    X = R[j];  while( F[X]) { ++j; X = R[j]; } // next unused value

    Y = X ^ period;  F[X] = F[Y] = k;
  }

  if( debug) printlu( "F", F, N);

  operate( q, q.n >> 1, f_FR, 0);

  free( R);

  return period;
}

// f_k1: classical set bit k high
//
unsigned long f_k1( unsigned long a, unsigned long b, unsigned long K)
{
  return b ^ (a | K);
}

// k1: (a,b) -> (a,b+f(a)), set bit k high
//
unsigned long k1( State q, unsigned int k)
{
  unsigned long period = 1LU << k;  operate( q, q.n >> 1, f_k1, period);

  return period;
}

// dot: dot product
//
int dot( unsigned long x, unsigned long y)
{
  int r = 0;  unsigned long z = x & y;

  while( z > 0) { if( z & 1) { r ^= 1; } z >>= 1; }

  return r;
}

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

int main( int argc, char *argv[])
{
 // random phase and debug options
 //
  int rphase = 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 'r': ++rphase; break;
      case 'd': ++debug; break;
      case 'E': ++Eflag; break;
      default: error( "Simon: bad option");
    }
    --argc, ++argv;
  }

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

  if( qn == 0 || (qn & 1) ) error( "Simon: qn must be non-zero and even");

 // bit number to set high, or negative to use random periodic function
 //
  int k = 0;  if( argc > 2) k = atoi(argv[2]);

  if( (k >= 0) && (k >= (qn >> 1)) ) error( "Simon: k out of range");

  State q = Q(qn);  printf( "Simon: seed = %u, qn = %u, k = %d\n", seed, qn, k);

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

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

 // apply H to the top half qubits
 //
  for( unsigned int i = q.n >> 1; i < q.n; ++i)
  {
    H( q, i);  if( debug) { print( "H", q, 0);  if( Eflag) print( "E", q, -1); }
  }

 // apply the periodic function
 //
  unsigned long period;

  if( k < 0)
  {
    period = FR( q);  if( debug) { print( "FR", q, 0);  if( Eflag) print( "E", q, -1); }
  }
  else
  {
    period = k1( q, k);  if( debug) { print( "k1", q, 0);  if( Eflag) print( "E", q, -1); }
  }

  printf( "period = %lu\n", period);

 // apply H again to the top half qubits
 //
  for( unsigned int i = q.n >> 1; i < q.n; ++i)
  {
    H( q, i);  if( debug) { print( "H", q, 0);  if( Eflag) print( "E", q, -1); }
  }
 
 // now only states (a,b) with a orthogonal to the period should be non-zero
 //
  if( debug) printf( "dot products with the period\n");

  unsigned int shift = q.n >> 1;  double p = 0;

  for( unsigned long i = 0; i < q.N; ++i)
  {
    int d = dot(i>>shift,period);
    if( d == 0) p += abs2(q.a[i]);
    if( debug) printf( "%d = %lu:(%g,%g)\n", d, i, creal(q.a[i]), cimag(q.a[i]) );
  }

  printf( "P(orthogonal to period) = %g, 1-P = %g\n", p, 1-p);

  return 0;
}