#include "stdafx.h"
#include <iostream>
#include <time.h>
#include <complex>

void matmulpnoif(std::complex<float> arr[], std::complex<float> out[], int numqbits, std::complex<float> a,
std::complex<float> b, std::complex<float> c, std::complex<float> d, int target)
{
    long length = 1 << (numqbits);
    long offset = 1 << (target - 1);
    long state = 0;
    while (state < length)
    {
        out[state] = arr[state] * a + arr[state + offset] * b;
        out[state + offset] = arr[state] * c + arr[state + offset] * d;
        state += 1 + offset * (((state%offset) + 1) / offset);
    }
}

void matmulpsingle(std::complex<float> arr[], std::complex<float> out[], int numqbits, std::complex<float> a,
std::complex<float> b, std::complex<float> c, std::complex<float> d, int target)
{
    long length = 1 << (numqbits);
    int shift = target - 1;
    long offset = 1 << shift;
    for (long state = 0; state < length; ++state)
    {
        if ((state >> shift) & 1)
        {
            out[state] = arr[state - offset] * c + arr[state] * d;
        }
        else
        {
            out[state] = arr[state] * a + arr[state + offset] * b;
        }
    }
}

int main()
{
    using namespace std;
    int numqbits = 25;
    long arraylength = 1 << numqbits;
    complex<float>* amplitudes = new complex<float>[arraylength];
    for (long i = 0; i < arraylength; ++i)
    {
        amplitudes[i] = complex<float>(0., 0.);
    }
    amplitudes[0] = complex<float>(1., 0.);
    complex<float> a(0., 0.);
    complex<float> b(1., 0.);
    complex<float> c(0., 0.);
    complex<float> d(1., 0.);
    int target = 1;
    int repititions = 10;
    clock_t startTime;
    //while (target <= numqbits) {
        startTime = clock();
        for (int j = 0; j < repititions; ++j) {
            complex<float>* outputs = new complex<float>[arraylength];
            matmulpsingle(amplitudes, outputs, numqbits, a, b, c, d, target);
            delete[] outputs;
        }
        cout << float(clock() - startTime) / (float)(CLOCKS_PER_SEC*repititions) << " seconds." << endl;
        startTime = clock();
        for (int k = 0; k < repititions; ++k) {
            complex<float>* outputs = new complex<float>[arraylength];
            matmulpnoif(amplitudes, outputs, numqbits, a, b, c, d, target);
            delete[] outputs;
        }
        cout << float(clock() - startTime) / (float)(CLOCKS_PER_SEC*repititions) << " seconds." << endl;
        target+=1;
    //}
    delete[] amplitudes;
    return 0;
}