// quantum computing emulation library
//
// R. Perry, Jan. 2018
//
#define _QCE_C_

#include <stdio.h>
#include <tgmath.h>
#include <stdlib.h> // rand(), srand()
#include <string.h> // memcpy()
#include <limits.h> // CHAR_BIT
#include <float.h>  // DBL_EPSILON
#include <time.h>   // time()
#include "qce.h"

//------------------------------ private ------------------------------

//  bitrev	bit reverse, used by qfft()
//
static unsigned long bitrev( unsigned long x, unsigned int nbits)
{
  unsigned long r = 0;

  for( unsigned int i = 0; i < nbits; ++i) { r <<= 1; r |= x & 1; x >>= 1; }

  return r;
}

// if C11 CMPLX macro is not available
// https://gcc.gnu.org/onlinedocs/gcc-8.2.0/gcc/Other-Builtins.html#index-_005f_005fbuiltin_005fcomplex
//
#ifndef CMPLX
#define CMPLX(x,y) __builtin_complex((x),(y))
#endif

//------------------------------ public ------------------------------

// pi and s2 values created by main/constants.c
//
const double pi = 0x1.921fb54442d18p+1; // 3.14159265358979312 = 4*atan(1)
const double s2 = 0x1.6a09e667f3bcdp-1; // 0.707106781186547573 = sqrt(2)/2

// set srand() using time(0) or SRAND environment variable
//
unsigned int SRAND( void)
{
  char *env = getenv( "SRAND");  unsigned int seed = 0;  if( env) seed = atoi(env);

  if( seed == 0) seed = time(0);

  srand( seed);  return seed;
}

//-------------------- functions for complex only --------------------

// Hadamard transformation on qubit k with non-resonant pulse error
// see http://fog.misty.com/perry/qc/pulse/notes.html
// ra = a/d would be 0 for no error
//
void Hra( State_Complex q, unsigned int k, double ra)
{
  if( k >= q.n) error( "H: bad index");

  double rb = sqrt(1-ra*ra), A = pi/(4*rb), C = cos(A), S = sin(A);

  double raS = ra*S, rbS = rb*S;

  double complex c00 = CMPLX(C,-raS), c11 = CMPLX(C,raS);

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

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

    q.a[i] = c00*t0 + rbS*t1;  q.a[j] = rbS*t0 - c11*t1;

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

// Hadamard transformation on qubit k with phase shift error
// see http://fog.misty.com/perry/qc/pulse/notes.html
// Note also that Hrp*Hrp = I
// rp = p/(pi/2) would be 1 for no error
//
void Hrp( State_Complex q, unsigned int k, double rp)
{
  if( k >= q.n) error( "H: bad index");

  double p = rp*(pi/2), C = cos(p), S = sin(p);

 // printf( "C = %g, 1-C = %g, S = %g, 1-S = %g\n", C, 1-C, S, 1-S);
 // with rp=1.0: C = 6.12323e-17, 1-C = 1, S = 1, 1-S = 0
 // use more accurate values for this case:
 //
  if( rp == 1.0) { C = 0; S = 1; }

  double complex c01 = CMPLX(S,-C), c10 = CMPLX(S,C);

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

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

    q.a[i] = s2*(t0 + c01*t1);  q.a[j] = s2*(c10*t0 - t1);

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

// quantum FFT
//
void qfft( State_Complex q, int inv)
{
  fft( q.a, q.N, inv);  double s = 1/sqrt(q.N);

 // scale values and bit-reverse the indices
 // see https://stackoverflow.com/questions/932079/in-place-bit-reversed-shuffle-on-an-array
 //
  for( unsigned long i = 0, j = 0; i < q.N; ++i)
  {
    j = bitrev( i, q.n);  

    if( i < j) { swap_Complex( q, i, j);  q.a[i] *= s;  q.a[j] *= s; }

    else if( i == j) q.a[i] *= s;
  }
}
  
// print in binary, n bits
//
void printb( unsigned long x, unsigned int n)
{
  unsigned long mask = 1LU << (n-1);

  while( mask) { putchar( (x & mask) ? '1' : '0'); mask >>= 1; }
}

//------------------------- general functions -------------------------

// print error message and exit
//
void error( const char *msg)
{
  fprintf( stderr, "%s\n", msg);  exit(1);
}

// malloc with exit on error
//
void *emalloc( unsigned long size)
{
  void *r = malloc(size);  if( !r) { error( "qce: malloc() failed"); }  return r;
}

//--------- _Complex and _Real functions defined by templates ---------
//
#undef QTYPE_COMPLEX
#define QTYPE_COMPLEX
#undef QTYPE
#undef qtype
#undef State
#define QTYPE Complex
#define qtype qtype_Complex
#define State State_Complex
#include "template.c"
//
#undef QTYPE_COMPLEX
#undef QTYPE
#undef qtype
#undef State
#define QTYPE Real
#define qtype qtype_Real
#define State State_Real
#include "template.c"