Simple quantum computer emulator

main.cpp

#include "stdafx.h"

template<size_t N>
class qgate {
public:
	static const size_t c_inputs = N;
};

class qstate {
private:
	vector<complex<double>> m_ket;
	size_t m_qubits;

public:
	explicit qstate(size_t N, size_t basis = 0)
		: m_ket((size_t) 1u << N, 0.0) {
		m_ket[basis] = 1.0;
		m_qubits = N;
	}

	template<typename T>
	void apply(T const& gate, array<size_t, T::c_inputs> const& indices) {
		vector<complex<double>> newKet(m_ket.size(), 0.0);

		size_t mask = 0;
		for (size_t k = 0; k < T::c_inputs; ++k) {
			mask |= (size_t) 1u << indices[k];
		}
		mask = ~mask;

		for (size_t i = 0; i < m_ket.size(); ++i) {
			for (size_t j = 0; j < m_ket.size(); ++j) {
				if ((i ^ j) & mask)
					continue;

				size_t input = 0, output = 0;
				for (size_t k = 0; k < T::c_inputs; ++k) {
					input |= ((i >> indices[k]) & 1) << k;
					output |= ((j >> indices[k]) & 1) << k;
				}
				newKet[i] += gate.getCoefficient(input, output)*m_ket[j];
			}
		}

		m_ket = newKet;
	}

	inline size_t qubits() const {
		return m_qubits;
	}
	inline size_t size() const {
		return m_ket.size();
	}
	inline complex<double> operator[](size_t basis) const {
		return m_ket[basis];
	}
};

inline ostream& operator<<(ostream& os, qstate const& state) {
	for (size_t i = 0; i < state.size(); ++i) {
		auto v = state[i];

		os << setw(20) << v << " |";
		for (size_t j = 0; j < state.qubits(); ++j) {
			os << (((i >> j) & 1) ? '1' : '0');
		}
		double p = v.real()*v.real() + v.imag()*v.imag();
		os << ">   " << (p*100.0) << "%" << endl;
	}
	return os;
}

class qnot : public qgate<1> {
public:
	inline complex<double> getCoefficient(size_t input, size_t output) const {
		return (input ^ output) ? 1.0 : 0.0;
	}
};

class qhadamard : public qgate<1> {
public:
	inline complex<double> getCoefficient(size_t input, size_t output) const {
		return (input & output) ? -M_SQRT1_2 : M_SQRT1_2;
	}
};

class qsrn : public qgate<1> {
public:
	inline complex<double> getCoefficient(size_t input, size_t output) const {
		return (input & ~output) ? -M_SQRT1_2 : M_SQRT1_2;
	}
};

class qrotation : public qgate<1> {
private:
	double s, c;
public:
	inline qrotation(double angle) {
		s = sin(angle);
		c = cos(angle);
	}
	inline complex<double> getCoefficient(size_t input, size_t output) const {
		return (input ^ output) ? (input ? s : -s) : c;
	}
};

class qphase : public qgate<1> {
private:
	complex<double> m_phase;
public:
	inline qphase(double angle)
		: m_phase(cos(angle), sin(angle)) {
	}
	inline complex<double> getCoefficient(size_t input, size_t output) const {
		return (input == output) ? (input ? m_phase : 1.0) : 0.0;
	}
};

class cnot : public qgate<2> {
public:
	inline complex<double> getCoefficient(size_t input, size_t output) const {
		return ((input ^ output) == ((input & 1) << 1)) ? 1.0 : 0.0;
	}
};

class cphase : public qgate<2> {
private:
	complex<double> m_phase;
public:
	inline cphase(double angle)
		: m_phase(cos(angle), sin(angle)) {
	}
	inline complex<double> getCoefficient(size_t input, size_t output) const {
		return (input == output) ? ((input == 3) ? m_phase : 1.0) : 0.0;
	}
};

class qswap : public qgate<2> {
public:
	inline complex<double> getCoefficient(size_t input, size_t output) const {
		return ((((input&1) << 1)|(input >> 1)) == output) ? 1.0 : 0.0;
	}
};

class cswap : public qgate<3> {
public:
	inline complex<double> getCoefficient(size_t input, size_t output) const {
		if (input & 1)
			return ((((input & 2) << 1) | ((input & 4) >> 1) | (input & 1)) == output) ? 1.0 : 0.0;
		else
			return (input == output) ? 1.0 : 0.0;
	}
};

class qtoffoli : public qgate<3> {
public:
	inline complex<double> getCoefficient(size_t input, size_t output) const {
		return ((input ^ output) == ((input & (input >> 1) & 1) << 2)) ? 1.0 : 0.0;
	}
};

template<size_t N>
class qft : public qgate<N> {
private:
	double m_scale;
	complex<double> m_root;

public:
	inline qft() {
		size_t n = (size_t) 1u << N;
		m_scale = 1.0 / sqrt(n);
		double theta = M_PI*2.0 / n;
		m_root = complex<double>(cos(theta), sin(theta));
	}

	inline complex<double> getCoefficient(size_t input, size_t output) const {
		return m_scale*pow(m_root, input*output);
	}
};

class qhalf_adder : public qgate<2> {
public:
	inline complex<double> getCoefficient(size_t input, size_t output) const {
		return (output == ((input+1)>>1)) ? 1.0 : 0.0;
	}
};

int _tmain(int argc, _TCHAR* argv[])
{
	qstate state(3);
	cout << state << endl;

	state.apply(qnot(), { 1 });
	cout << state << endl;

	state.apply(qft<3>(), { 0, 1, 2 });
	cout << state << endl;

	state.apply(qft<3>(), { 0, 1, 2 });
	cout << state << endl;

	cin.get();

	return 0;
}

stdafx.h

// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//

#pragma once

#include <SDKDDKVer.h>

#include <stdio.h>
#include <tchar.h>

#define NOMINMAX
#define _USE_MATH_DEFINES

#include <Windows.h>

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <map>
#include <unordered_map>
#include <set>
#include <unordered_set>
#include <functional>
#include <algorithm>
#include <stdint.h>
#include <complex>
#include <bitset>
#include <array>
#include <iomanip>
#include <math.h>

using namespace std;