// quantum computing emulation routines using real coefficients
//
// Deutsch-Jozsa problem: determine if a hidden Boolean function is constant or balanced
//
// R. Perry, Jan. 2019
//
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include "qc.h"

//------------------------------ private ------------------------------
//
//  s2	    sqrt(2)/2
//  choice  specify which hidden function to use: 0, 1, 2, 3, 4, 5, 6
//  parity  XOR of the bits
//  swap    swap two states
//  Hk	    perform Hadamard transformation on qubit k
//
static double s2;
//
static int choice( void)
{
  static int f = -1;  if( f < 0) { srand(time(0)); f = rand()%7; }

  return f;
}
//
static int parity( int x)
{
  int r = 0;  while( x) { r ^= (x & 1); x /= 2; }  return r;
}
//
static void swap( double q[], int i, int j)
{
  double t = q[i];  q[i] = q[j];  q[j] = t;
}
//
static void Hk( double q[], int M, int k)
{
  if( s2 == 0) s2 = sqrt(2)/2;

  int i = 0, j, k1 = 1 << k; // k1 = 0...010...0, k'th bit = 1

  while( i < M) // i will have k'th bit = 0, j will have k'th bit = 1
  {
    j = i | k1;  double t0 = q[i], t1 = q[j];

    q[i] = s2*(t0+t1);  q[j] = s2*(t0-t1);

    ++i;  i += (i & k1); // skip i values with k'th bit = 1
  }
}

//------------------------------ public ------------------------------
//
//  fx	  the classical hidden Boolean function
//  fq	  the quantum hidden Boolean function
//  H	  perform Hadamard transformation on all qubits
//
int fx( int x, int N) // 0 <= x <= N-1
{
  int mid = N/2;

  switch( choice() )
  {
    case 0: return 0;		// constant 0
    case 1: return 1;		// constant 1
    case 2: return (x < mid);	// x < N/2
    case 3: return !(x < mid);	// x >= N/2
    case 4: return x & 1;	// odd
    case 5: return !(x & 1);	// even
    case 6: return parity(x);	// parity
  }

  return 0;
}
//
void fq( double q[], int M) // state (x,y) -> (x,y+f(x))
{
  int mid = M/2;

  switch( choice() )
  {
    case 0: // f(x) = 0, do nothing
	break;
    case 1: // f(x) = 1, XOR last bit, i.e. swap adjacent
	for( int i = 0; i < M; i += 2) swap( q, i, i|1);
	break;
    case 2: // f(x) = (x < N/2)
	for( int i = 0; i < mid; i += 2) swap( q, i, i|1);
	break;
    case 3: // f(x) = (x >= N/2)
	for( int i = mid; i < M; i += 2) swap( q, i, i|1);
	break;
    case 4: // f(x) = (x & 1) // odd
	for( int i = 2; i < M; i += 4) swap( q, i, i|1);
	break;
    case 5: // f(x) = !(x & 1) // even
	for( int i = 0; i < M; i += 4) swap( q, i, i|1);
	break;
    case 6: // f(x) = parity(x)
	for( int i = 0; i < M; i += 2) if( parity(i/2)) swap( q, i, i|1);
	break;
  }
}
//
void H( double q[], int M)
{
  int mask = M/2;

  for( int k = 0; mask > 0; ++k) { Hk( q, M, k); mask /= 2; }
}