// Quantum algorithms namespace implementation file

#include <stdafx.h>
#include "resource.h"

#include<math.h>
#include<iostream.h>
#include<fstream.h>
#include<iomanip.h>
#include<stdlib.h>
#include<time.h> 

#include "qalgos.h"

// use functions from the CCAlgos and CModexp namespaces
// without the need for the :: operartor
using namespace CCAlgos;
using namespace CModexp;

void CQAlgos::qft(CHilbert* hilb, THREADPARAMS* ptp)
{
	// number of qubits to operate qft on 
	int n_qbs = hilb->numbits();

	float percentage = 0;

	// num gates in qft circuit
	int numgates = (n_qbs)*(n_qbs - 1)/2;	

	// counter to number of gates applied during the qft
	int gate = 0;

	// qubits upon which the gates act
	int butterfly_index[1];
	int twiddle_bits[2];
	
	// number of qubits to operate on
	int l = n_qbs;
	int i, j;

	// always same matrix representation -> only
	// need to create once
	CButterfly* bgate = new CButterfly();

	// create and apply circuit
	for(i = l - 1; i > 0; i--)
	{
		butterfly_index[0] = i;

		bgate->applysub(hilb, butterfly_index);

		for(j = 0; j < l-i; j++)
		{
			twiddle_bits[0] = i-1;
			twiddle_bits[1] = l-j-1;

			// need to create this each time as action
			// is qubit dependant
			CTwiddle* twiddlegate = new CTwiddle(twiddle_bits[0], twiddle_bits[1]);

			twiddlegate->applysub(hilb, twiddle_bits);
			
			delete twiddlegate;
		}
	}

	butterfly_index[0] = 0;

	bgate->applysub(hilb, butterfly_index);
	
	delete bgate;
	
	hilb->reverse();	// reverse as reulsts are presneted in 
						// reverse basis state order

//	cout << "Qft complete" << endl;
}

void CQAlgos::subqft(int s_i, int e_i, CHilbert* hilb, THREADPARAMS* ptp)
{
	// (MFC) get activity progress object ptr
	CProgressCtrl *progress_ctrl_ptr = (CProgressCtrl*) ptp->pView->GetDlgItem(IDC_ACTIVITY_PROGRESS);
	progress_ctrl_ptr->SetPos(0);

	float percentage = 0;
	
	// num gates in the qft circuit for this number of qubits
	int numgates = (e_i - s_i + 1)*(e_i - s_i + 2)/2;
	
	// counter to number of gates applied during the qft
	int gate = 0;

	double progress = 0;
	
	// qubits upon which the gates act
	int butterfly_index[1];
	int twiddle_bits[2];

	// number of qubits upon which to perform the qft
	int l = e_i-s_i+1;
	int i, j;

	// always same matrix representation -> only
	// need to create once
	CButterfly* bgate = new CButterfly();

	for(i = l-1; i > 0; i--)
	{
		butterfly_index[0] = s_i + i;
		bgate->applysub(hilb, butterfly_index);
		gate++;

		// (MFC)
		progress = ((float)gate/numgates)*100.0;
		progress_ctrl_ptr->SetPos((int)progress);

		for(j = 0; j < l-i; j++)
		{

			twiddle_bits[0] = s_i+i-1;
			twiddle_bits[1] = s_i+l-j-1;

			// need to create this each time as action
			// is qubit dependant
			CTwiddle* twiddlegate = new CTwiddle(twiddle_bits[0], twiddle_bits[1]);
			
			twiddlegate->applysub(hilb, twiddle_bits);
			gate++;

			// (MFC)
			progress = ((float)gate/numgates)*100.0;
			progress_ctrl_ptr->SetPos((int)progress);
		
			delete twiddlegate;
		}
	}

	butterfly_index[0] = s_i;

	bgate->applysub(hilb, butterfly_index);
	gate++;

	// (MFC)
	progress = 100;
	progress_ctrl_ptr->SetPos((int)progress);

	delete bgate;

	hilb->reverse();	// reverse as results are presented in 
						// reverse basis state order

//	cout << "Qft complete" << endl;
}

bool CQAlgos::shor(int n, THREADPARAMS* ptp, int fileout)
{
	// (MFC) get overall sim progress object ptr
	CProgressCtrl *progress_ctrl_ptr = (CProgressCtrl*) ptp->pView->GetDlgItem(IDC_SIM_PROGRESS);
	progress_ctrl_ptr->SetPos(0);

	// register initialisation
	// form is |reg2|reg1> |F(A)|A >
	// note that in the literature reg1 is usualy 2log n*n / log 2
	// shorter register gives more spread in QFT but much faster
	// run time

	// register sizes
	int n_qbs_reg1 = (int)ceil(log(n)/log(2));		
	int n_qbs_reg2 = (int)ceil(log(n)/log(2));

	// number of qubits in the entire space
	int n_qbs_entire_space = n_qbs_reg1 + n_qbs_reg2;

	// integer lists representing the qubits in each register
	int* reg1 = new int[n_qbs_reg1];
	int* reg2 = new int[n_qbs_reg2];

	// assign qubit indices to the registers
	int k = 0;
	for(int i = 0; i < n_qbs_reg1; i++)
		reg1[i] = i;
	for(int j = n_qbs_reg1; j < n_qbs_entire_space; j++)
		reg2[k++] = j;

	// create statevector space for the entire system
	CHilbert* state_vec = new CHilbert(n_qbs_entire_space);

	// get coprime x
	int x = coprime(n);

	// (MFC)
	CString string;
	string.Format("Coprime x: %d", x);
	ptp->res->activity_string = string;
	ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	
	
	// (MFC)
	if(fileout)
		ptp->fout << "Performing modular exponentiation on reg1, result in reg2 " << endl;
	ptp->res->activity_string = "Performing modular exponentiation on reg1, result in reg2 ";
	ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	

	// peform modula exponentiation on reg1, result in reg2
	state_vec->modexp(n_qbs_reg1, n_qbs_reg2, x, n);

	// (MFC)
	progress_ctrl_ptr->SetPos(10);
	if(fileout)
	{
		ptp->fout << *state_vec;
		ptp->fout << "Measuring reg2..." << endl; 
	}	
	ptp->res->activity_string = "Measuring reg2...";
	ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	

	int reg2val = state_vec->measure(n_qbs_reg2, reg2, ptp);

	// (MFC)
	string.Format("Value : %d", reg2val);
	ptp->res->activity_string = string;
	ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	
	if(fileout)
		ptp->fout << *state_vec;
	progress_ctrl_ptr->SetPos(55);
	if(fileout)
		ptp->fout << "Getting space for reg1, ignoring reg2" << endl;
	ptp->res->activity_string = "Getting space for reg1, ignoring reg2";
	ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	
	
	// ptr to a smaller space of reg1 
	CHilbert* sub_state_vec;

	// get subspace representing reg1
	// i.e. project entire space onto reg1 subspace
	sub_state_vec = state_vec->subspaceprojection(n_qbs_reg1, reg1);
	
	// (MFC)
	progress_ctrl_ptr->SetPos(65);

	// tidy up
	delete state_vec;

	// (MFC)
	if(fileout)
		ptp->fout << *sub_state_vec;
	if(fileout)
		ptp->fout << "Performing qft on reg1" << endl;
	ptp->res->activity_string = "Performing qft on reg1";
	ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	

	// perform qft on reg1, since this is the entire space now
	// we could just use qft, but the MFC progress stuff hasn't been 
	// added to it yet
	subqft(0, n_qbs_reg1-1, sub_state_vec, ptp);

	// (MFC)
	progress_ctrl_ptr->SetPos(95);
	if(fileout)
	{
		ptp->fout << *sub_state_vec;
		ptp->fout << "Measuring reg1..." << endl; 
	}
	ptp->res->activity_string = "Measuring reg1...";
	ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	

	int reg1val = sub_state_vec->measure(n_qbs_reg1, reg1, ptp);

	// (MFC)
	string.Format("Value : %d", reg1val);
	ptp->res->activity_string = string;
	ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	
	if(fileout)
		ptp->fout << *sub_state_vec;
	progress_ctrl_ptr->SetPos(80);

	// tidy up
	delete sub_state_vec;
	delete reg1;
	delete reg2;

	// failure indicators
	bool failed = false;
	bool facfailed[2];
	facfailed[0] = facfailed[1] = false;

	// test to see if valid value measured
	// 0 -> all period information lost
	if(reg1val == 0)
	{
		cout << endl << "******** Zero measured in reg1 ********" << endl << endl;
		ptp->res->activity_string = "******** Zero measured in reg1 ********";
		ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	
		failed = true;
	}
	else
	{
		cout << "Using rational approximation to find order";
		ptp->res->activity_string = "Using rational approximation to find order";
		ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	
		cout << reg1val << " and " << (int)pow(2,n_qbs_reg1) << endl;
		cout << "Approx is: ";

		// do rational approximation to find the order
		int r = approxdenom((double)reg1val/pow(2, n_qbs_reg1), n);
		cout << "Order is: " << r << endl; 

		// calculated order must be even 
		if(r%2 != 0)
		{
			cout << endl << "******** Odd order ********" << endl << endl;
			ptp->res->activity_string = "******** Odd order ********";
			ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	
			failed = true;
		}
		else
		{
			// slight variation on the original Shor algorithm
			// we need the gcd of y+-1 and n -> euclids algorithm will
			// first divide by n -> do modexp to cut out this step
			int y = modexp(x, r/2, n);
			cout << "Possible factors are gcd(" << y << "+-1," << n << ")" << endl;
	
			int factors[2];
	
			// calculate potential factors
			factors[0] = euclid(y-1, n);
			factors[1] = euclid(y+1, n);

			cout << "Possible factors are: " << factors[0] << " and " << factors[1] << endl;
			
			// (MFC)
			CString string;
			string.Format("Possible factors are: %d and %d", factors[0], factors[1]);
			ptp->res->activity_string = string;
			ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	

			// check if potential factors are indeed factors
			for(i = 0; i < 2; i++)
			{
				if(factors[i] == n || factors[i] == 1)
				{
					string.Format("******** %d - trivial factor ********", factors[i]);
					ptp->res->activity_string = string;
					ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);	
					facfailed[i] = true;
				}
				else if(!(n % factors[i])) 
				{
					string.Format("******** %d - confirmed factor ********", factors[i]);
					ptp->res->activity_string = string;
					ptp->pWnd->SendMessage(WM_USER_THREAD_UPDATE_BOX);						
				}

			}

		}
	}

	// (MFC)
	progress_ctrl_ptr->SetPos(100);

	if(failed || (facfailed[0] && facfailed[1]))
		return false;	// both factors are invalid
	else
		return true;	// all okay
}