quantum.cpp

#include "quantum.hpp"
#include <cmath>
#include <cstdint>
#include <random>
#include <stdexcept>
#include <immintrin.h> //optimization of bit scattering (_pdep_u64)


std::random_device rnd_device;
std::mt19937 rnd_engine{rnd_device()};

inline qufloat_t rand_float(qufloat_t a = 0, qufloat_t b = 1) {
	std::uniform_real_distribution<qufloat_t> dist(a, b);
	return dist(rnd_engine);
}


std::vector<uint8_t> bit_range(uint8_t start, uint8_t end) {
	std::vector<uint8_t> v;
	for (uint8_t i = start; i <= end; i++) v.push_back(i);
	return v;
}

std::ostream& operator<<(std::ostream& os, const std::vector<amp_t>& vec) {
	os << '[';
	for (size_t i = 0; i < vec.size(); i++) {
		os << vec[i];
		if (i != vec.size() - 1)
			os << ", ";
	}
	os << ']';
	return os;
}

std::string print_bits(uint64_t x, int bit_num) {
	std::string str = "";
	bool printing = (bit_num != 0);
	const int start = bit_num == 0 ? 0 : 64 - bit_num;
	
	for (int i = start; i < 64; i++) {
		const bool b = x >> (63-i) & 1;
		if (b || i >= 63) printing = true;
		if (printing)
			str += b ? '1' : '0';
	}
	return str;
}


void qureg::resize(uint_fast8_t num) {
	if (num == 0) throw std::out_of_range("Can't have 0 qubits");
	if (num > max_qubits) throw std::out_of_range("Quantum register is too big");
	size = num;
	state.resize(1 << num, amp_t());
}

qureg::qureg(uint8_t num) {
	resize(num);
	const uint64_t state_size = 1 << num;
	
	//state gets filled with equal magnitudes but random phases
	double amp = 1. / sqrt(state_size);
	for (amp_t& v : state)
		v = std::polar(amp, rand_float(0, 2*M_PI));
}

qureg::qureg(uint8_t num, uint64_t bits) {
	resize(num);
	const uint64_t state_size = 1 << num;
	
	if (bits >> num) throw std::out_of_range("Invalid collapsed state requested");
	
	state[bits] = amp_t(1, 0);
}


double qureg::probability(uint64_t bits) const {
	return norm(state[bits]);
}

uint64_t qureg::measure(void) const {
	//we'll look for state amplitude whose range this is in
	double x = rand_float();
	
	//Kaham sum to aid accuracy
	double sum = 0;
	double c = 0;
	
	const uint64_t len = state.size();
	for (uint64_t i = 0; i < len; i++) {
		double y = probability(i) - c;
		double t = sum + y;
		c = (t - sum) - y;
		
		sum = t;
		if (x <= sum) return i;
	}
	
	return state.size() - 1;
}

uint64_t qureg::collapse(void) {
	uint64_t m = measure();
	
	std::fill(state.begin(), state.end(), amp_t());
	state[m] = amp_t(1, 0);
	
	return m;
}


uint64_t scatter(uint64_t bits, uint64_t mask) {
	return _pdep_u64(bits, mask); //TODO add fallback
}

uint64_t reorder(uint64_t bits, const std::vector<uint8_t>& order) {
	uint64_t result = 0;
	for (uint8_t i = 0; i < order.size(); i++)
		result |= ((bits & (1 << i)) >> i) << order[i];
	return result;
}

void qureg::gate(qugate gate, bit_list bits, bit_list c_on, bit_list c_off) {
	const uint8_t num = bits.size();
	const uint8_t cbit_num = c_on.size() + c_off.size();
	if (num != gate.bits) throw std::invalid_argument("Given different amount of bits than gate accepts");
	if (gate.bits + cbit_num > size) throw std::invalid_argument("Gate is too big for register");
	
	//we break down the operation into the blocks where the data inside only affects itself
	//can be parallelized for systems with a lot of qubits
	
	uint64_t cmask = 0, cbits = 0;
	for (int i = 0; i < c_on.size(); i++)
		cbits |= 1 << c_on[i];
	for (int i = 0; i < c_off.size(); i++)
		cmask |= 1 << c_off[i];
	cmask |= cbits;
	
	const uint64_t blocks = 1LU << (size - num - cbit_num);
	const uint64_t block_size = 1LU << num;
	
	uint64_t mask = (~0LU >> (64 - size)) ^ cmask;
	for (int i = 0; i < num; i++) mask ^= 1LU << bits[i];
	
	//input and output of gate
	std::vector<amp_t*> in(block_size);
	std::vector<amp_t> out(block_size);
	
	for (uint64_t b = 0; b < blocks; b++) {
		const uint64_t block = scatter(b, mask) | cbits;
		
		//prepare block into inputs
		for (uint64_t i = 0; i < block_size; i++)
			in[i] = &state[block | reorder(i, bits)];
		
		gate(in, out);
		
		//write back
		for (uint64_t i = 0; i < block_size; i++)
			*in[i] = out[i];
	}
}


void gate_X(amp_t** in, amp_t* out, void* _) {
	out[0] = *in[1];
	out[1] = *in[0];
}
qugate qugate::pauli_X() {
	return qugate(gate_X);
}
void qureg::op_X(uint8_t bit) {
	gate(qugate::pauli_X(), {bit});
}

void gate_Y(amp_t** in, amp_t* out, void* _) {
	out[0] = amp_t(in[1]->imag(), -in[1]->real());
	out[1] = amp_t(-in[0]->imag(), in[0]->real());
}
qugate qugate::pauli_Y() {
	return qugate(gate_Y);
}
void qureg::op_Y(uint8_t bit) {
	gate(qugate::pauli_Y(), {bit});
}

void gate_Z(amp_t** in, amp_t* out, void* _) {
	out[0] = *in[0];
	out[1] = -*in[1];
}
qugate qugate::pauli_Z() {
	return qugate(gate_Z);
}
void qureg::op_Z(uint8_t bit) {
	gate(qugate::pauli_Z(), {bit});
}

void gate_hadamard(amp_t** in, amp_t* out, void* _) {
	constexpr const qufloat_t c = M_SQRT1_2;
	out[0] = c * (*in[0] + *in[1]);
	out[1] = c * (*in[0] - *in[1]);
}
qugate qugate::hadamard() {
	return qugate(gate_hadamard);
}
void qureg::op_hadamard(uint8_t bit) {
	gate(qugate::hadamard(), {bit});
}

void qureg::op_X_range(uint8_t a, uint8_t b) {
	qugate g = qugate::pauli_X();
	for (uint8_t i = a; i <= b; i++) gate(g, {i});
}
void qureg::op_Y_range(uint8_t a, uint8_t b) {
	qugate g = qugate::pauli_Y();
	for (uint8_t i = a; i <= b; i++) gate(g, {i});
}
void qureg::op_Z_range(uint8_t a, uint8_t b) {
	qugate g = qugate::pauli_Z();
	for (uint8_t i = a; i <= b; i++) gate(g, {i});
}
void qureg::op_hadamard_range(uint8_t a, uint8_t b) {
	qugate g = qugate::hadamard();
	for (uint8_t i = a; i <= b; i++) gate(g, {i});
}

void qureg::op_CX(uint8_t bit, uint8_t control) {
	gate(qugate::pauli_X(), {bit}, {control});
}
void qureg::op_CY(uint8_t bit, uint8_t control) {
	gate(qugate::pauli_Y(), {bit}, {control});
}
void qureg::op_CZ(uint8_t bit, uint8_t control) {
	gate(qugate::pauli_Z(), {bit}, {control});
}
void qureg::op_CCX(uint8_t bit, uint8_t c1, uint8_t c2) {
	gate(qugate::pauli_X(), {bit}, {c1, c2});
}
void qureg::op_CCY(uint8_t bit, uint8_t c1, uint8_t c2) {
	gate(qugate::pauli_Y(), {bit}, {c1, c2});
}
void qureg::op_CCZ(uint8_t bit, uint8_t c1, uint8_t c2) {
	gate(qugate::pauli_Z(), {bit}, {c1, c2});
}

void gate_swap(amp_t** in, amp_t* out, void* _) {
	out[0] = *in[0];
	out[1] = *in[2];
	out[2] = *in[1];
	out[3] = *in[3];
}
qugate qugate::swap() {
	return qugate(gate_swap, 2);
}
void qureg::op_swap(uint8_t a, uint8_t b) {
	gate(qugate::swap(), {a, b});
}


void gate_increment(amp_t** in, amp_t* out, void* data) {
	const int bits = ((int*)data)[0];
	const int add = ((int*)data)[1];
	
	for (int i = 0; i < (1 << bits); i++)
		out[i] = *in[(i - add) & ((1 << bits) - 1)];
}
qugate qugate::increment(int by, uint8_t bit_num) {
	int* dat = (int*)malloc(2 * sizeof(int));
	dat[0] = bit_num;
	dat[1] = by;
	return qugate(gate_increment, bit_num, dat, true);
}