// quantum computing emulation: Deutsch-Jozsa coin analogy
//
// Determine P(heads) for an NS-sided coin
//
// Usage: DJ-coin [-drE] [NS [NH]]
//
//  NS = number of sides
//  NH = number of heads
//  -d  extra debug output
//  -r  use random global phase factor
//  -E  show entanglement
//
// R. Perry, July 2018
//
#include <stdio.h>
#include <stdlib.h> // atol(), rand()
#include <tgmath.h>
#include <string.h> // memset()
#include "qce.h"

double pi;

// print X probabilities, without auxiliary bit
//
void print_X( const char *msg, State q)
{
  printf( "%s:\n", msg);

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

#if 0
// set state = |x>|1> where |x> is a superposition of 0...(NS-1)
//
void init_q_NS( State q, int rphase, unsigned long NS)
{
  double complex v = 1;  if( rphase) { double a = pi*(2*DRAND-1); v = exp(CMPLX(0,a)); }

  v /= sqrt(NS);  zero(q); 

  for( unsigned long i = 1; NS > 0; i += 2) { q.a[i] = v; --NS; }
}
#endif

// set state = H|0...0>|1>
//
void init_q( State q, int rphase)
{
  init( q, 1, rphase);  for( unsigned int k = 0; k < q.n; ++k) H( q, k);
}

// put array in random order
//
void shuffle( unsigned char R[], unsigned long n)
{
    unsigned char t;  unsigned long j;

    while( n > 1) { j = IRAND(n); --n; t = R[n]; R[n] = R[j]; R[j] = t; }
}

unsigned char *R; // for random and general (non-constant/non-balanced) functions
//
void print_R( unsigned long M)
{
  printf( "R: ");
  for( unsigned long i = 0; i < M; ++i) putchar( R[i] ? '1' : '0');
  putchar('\n');
}
//
// set R with NH 1's in first NS positions randomly
//
void init_R( unsigned long M, unsigned long NS, unsigned long NH)
{
  if( !R) R = emalloc(M);

  memset( R, 0, M);

  for( unsigned long i = 0; i < NH; ++i) R[i] = 1;

  shuffle( R, NS);
}

// f: the hidden Boolean function
//
// state (x,y) -> (x,y+f(x))
//
void f( State q)
{
  for( unsigned long i = 0; i < q.N; i += 2) if( R[i>>1]) swap( q, i, i|1);
}

// debug: display state
//
void disp( const char *msg, State q, int Eflag)
{
  print( msg, q, PRINT_BITS); if( Eflag) print( "E", q, PRINT_TANGLE);
}

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

int main( int argc, char *argv[])
{
  int debug = 0, rphase = 0, Eflag = 0;  pi = 4*atan(1);  unsigned int seed = SRAND();

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

 // number of sides and number of Heads
 //
  unsigned long NS = 2;  if( argc > 1) NS = atol(argv[1]); 

  if( NS < 1) error( "DJ-coin: bad NS arg");

  unsigned long NH = 1;  if( argc > 2) NH = atol(argv[2]);

  if( NH > NS) error( "DJ-coin: bad NH arg");

 // number of qubits is smallest q.n such that M = 2**(q.n-1) >= NS
 //
  unsigned int m = 1;  unsigned long M = 1LU << m;  while( M < NS) { ++m; M <<= 1; }

  State q = Q(m+1);

  printf( "DJ-coin: seed = %u, q.n = %u, q.N = %lu, M = %lu, NS = %lu, NH = %lu\n",
    seed, q.n, q.N, M, NS, NH);

  init_R( M, NS, NH);  print_R( M);

#if 0

 // init_q_NS effectively does H on the input qubits
 // 
  init_q_NS( q, rphase, NS);  if( debug) disp( "init_NS", q, Eflag);
 //
 // apply H to auxiliary qubit
 //
  H( q, 0);  if( debug) disp( "H", q, Eflag);

 // above not used; if NS < q.n then final H below produces mess

#else

  init_q( q, rphase);  if( debug) disp( "init", q, Eflag);

#endif

 // apply f
 //
  f( q);  if( debug) disp( "f", q, Eflag);

 // apply H
 //
 // can do all n bits (k=0,...) or skip the lsb (k=1,...)
 //
  for( unsigned int k = 0; k < q.n; ++k) { H( q, k); /* if( debug) disp( "H", q, Eflag); */ }

  if( debug) disp( "H", q, Eflag);

  print_X( "X", q);

  return 0;
}