diff --git a/exercises/ex05_solution/qsimulator.cpp b/exercises/ex05_solution/qsimulator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..09845d15d7d1e8e153fc528cb2ac0c07a0dbe390 --- /dev/null +++ b/exercises/ex05_solution/qsimulator.cpp @@ -0,0 +1,293 @@ +#include <vector> +#include <iostream> +#include <complex> +#include <string> +#include <random> + + +class QSimulator +{ +public: + typedef unsigned state_type; + typedef double scalar_type; + typedef std::complex<scalar_type> complex_type; + typedef std::vector<complex_type> wf_type; + typedef std::pair<scalar_type,scalar_type> dmatr_type; + QSimulator(unsigned n, bool v); + + void H(const wf_type& iwf, wf_type& owf, unsigned qubitId) const; + void CNOT(const wf_type& iwf, wf_type& owf, unsigned cqid, unsigned qid) const; + void Z(wf_type& wf, unsigned qid) const; + void X(const wf_type& iwf, wf_type& owf, unsigned qid) const; + void Uf(const wf_type& iwf, wf_type& owf, state_type(*f)(state_type st, unsigned hqid)) const; + + bool measureQubit(wf_type& wf, unsigned qid) const; // physical measurement + dmatr_type tomography(const wf_type& iwf, unsigned qid) const; // state of the qubit qid + void printWF(wf_type wf) const; + state_type maxState() const { return maxState_; } + +private: + const unsigned n_; + const state_type maxState_; + const bool verbose_; + + mutable std::mt19937 rndEngine_; + mutable std::uniform_real_distribution<scalar_type> probDist_; + mutable std::uniform_int_distribution<unsigned> qidDist_; +}; + +QSimulator::QSimulator(unsigned n, bool v) +: verbose_(v), n_(n), maxState_( (1<<n)-1 ), probDist_(0.0,1.0), qidDist_(1,n-1), rndEngine_(time(NULL)) +{ +} + +void QSimulator::H(const wf_type& iwf, wf_type& owf, unsigned qubitId) const{ + owf = iwf; + for(state_type st=0; st<=maxState(); st++){ + if(std::abs(iwf[st])!=0.0){ + state_type flippedSt = st ^ (1<<qubitId); + int sign = (st & (1<<qubitId)) >> qubitId; // 1 if qubitId is set 0 otherwise + sign = 2*sign-1; // 1 if qubitId is set -1 otherwise + + owf[flippedSt] = ( static_cast<scalar_type>(sign)*iwf[flippedSt] + iwf[st] )/std::sqrt(2); + owf[st] = ( -static_cast<scalar_type>(sign)*iwf[st] + iwf[flippedSt] )/std::sqrt(2); + } + } + if(verbose_){ + std::cout << "wavefunction after the Hadamard gate:" << std::endl; + printWF(owf); + } +} + +void QSimulator::CNOT(const wf_type& iwf, wf_type& owf, unsigned cqid, unsigned qid) const{ + owf = iwf; + for(state_type st=0; st<=maxState(); st++){ + if(std::abs(iwf[st])!=0.0){ + if (st & (1<<cqid)){ // control qubit set + state_type flippedSt = st ^ (1<<qid); // flip the other qubit + owf[flippedSt] = iwf[st]; + owf[st] = iwf[flippedSt]; + } + } + } + if(verbose_){ + std::cout << "wavefunction after the CNOT gate:" << std::endl; + printWF(owf); + } +} + +// ====== convention: output qubit id=0 +void QSimulator::Uf(const wf_type& iwf, wf_type& owf, state_type(*f)(state_type st, unsigned hqid) ) const{ + owf = iwf; + unsigned hiddenControlBitId = qidDist_(rndEngine_); + for(state_type st=0; st<=maxState(); st++){ + if(std::abs(iwf[st])!=0.0){ + if( f(st,hiddenControlBitId) ){ + state_type newState = st ^ 1; // y = y XOR f(x) + owf[st] = iwf[newState]; + owf[newState] = iwf[st]; + } + } + } + if(verbose_){ + std::cout << "owavefunction after the Uf gate:" << std::endl; + printWF(owf); + } +} + +bool QSimulator::measureQubit(wf_type& wf, unsigned qid) const{ + dmatr_type dmatr = tomography(wf, qid); + + scalar_type randNum = probDist_(rndEngine_); + bool measured1 = (randNum<dmatr.second)? 1: 0; + if(verbose_){ + std::cout << "probability to measure |0>: " << dmatr.first << ", for |1>: " << dmatr.second << std::endl; + std::cout << "randNum = " << randNum << ",measured1 = " << measured1 << std::endl; + } + +// ====== drop projected out components, renormalize the rest + for(state_type st=0; st<=maxState(); st++){ + if(std::abs(wf[st])!=0.0){ + if(st & (1<<qid)){ // bit set + if(measured1){ + wf[st] /= std::sqrt(dmatr.second); // renormalize + }else{ + wf[st] = 0; + } + }else{ // bit not set + if(measured1){ + wf[st] = 0; + }else{ + wf[st] /= std::sqrt(dmatr.first); // renormalize + } + } + } + } +// ====== drop projected out components, renormalize the rest + + if(verbose_){ + std::cout << "wavefunction after the measurement:" << std::endl; + printWF(wf); + } + + if(measured1){ + return true; + }else{ + return false; + } +} + +QSimulator::dmatr_type QSimulator::tomography(const wf_type& wf, unsigned qid) const{ + scalar_type beta2 = 0; + scalar_type alpha2 = 0; + for(state_type st=0;st<=maxState(); st++){ + if(st & (1<<qid)){ // bit set + beta2 += std::real( std::conj(wf[st]) * wf[st] ); + }else{ + alpha2 += std::real( std::conj(wf[st]) * wf[st] ); + } + } + return dmatr_type(alpha2,beta2); +} + +void QSimulator::X(const wf_type& iwf, wf_type& owf, unsigned qid) const{ // NOT gate + owf = iwf; + for(state_type st=0; st<=maxState(); st++){ + state_type newState = st ^ (1<<qid); + owf[newState] = iwf[st]; + owf[st] = iwf[newState]; + } + if(verbose_){ + std::cout << "wavefunction after the X gate:" << std::endl; + printWF(owf); + } +} + +void QSimulator::Z(wf_type& wf, unsigned qid) const{ + for(state_type st=0; st<=maxState(); st++){ + if(std::abs(wf[st])!=0.0){ + if(st & (1<<qid)){ + wf[st] *= -1; + } + } + } + if(verbose_){ + std::cout << "wavefunction after the Z gate:" << std::endl; + printWF(wf); + } +} + +void QSimulator::printWF(wf_type wf) const{ + for(state_type st=0; st<=maxState(); st++){ + std::cout << wf[st] << "\t"; + } + std::cout << std::endl; +} + + +QSimulator::state_type f_const0(QSimulator::state_type st, unsigned hiddenControlBitId){ + return 0; +} +QSimulator::state_type f_const1(QSimulator::state_type st, unsigned hiddenControlBitId){ + return 1; +} +QSimulator::state_type f_balanced(QSimulator::state_type st, unsigned hiddenControlBitId){ + return (st & 1<<hiddenControlBitId)>>hiddenControlBitId; +} + +std::string usage(const std::string& prog) +{ + return "usage: " + prog + " nQubitsForDeutschJosza"; +} +int main(int argc, const char** argv) +{ + if( argc != 2 ){ + std::cerr << usage(argv[0]) << std::endl; + return -1; + } + unsigned DJN = stoi(argv[1],nullptr,10); + if(DJN<2){ + std::cerr << "At least 2 qubits will be needed for the Deutsch-Josza algorithm" << std::endl; + return -1; + + } +// ============== Teleportation begin + // qubit index convention: unknown id=0, Alice Bell id=1, Bob Bell id=2 + bool verbose = false; + int n = 3; + std::cout << "Teleportation:" << std::endl; + QSimulator teleport(n, verbose); + QSimulator::wf_type iwf(teleport.maxState()+1,0); + QSimulator::wf_type owf(teleport.maxState()+1,0); + + iwf[1] = 1/3.0; + iwf[0] = 2*std::sqrt(2)/3; // 1/3 |001> + 2sqrt(2)/3 |000> + + std::cout << "initial wavefunction:" << std::endl; + teleport.printWF(iwf); + + QSimulator::dmatr_type dmatr = teleport.tomography(iwf, 0); + std::cout << "State of the Alice's |psi> qubit (id=0) is " << std::sqrt(dmatr.first) << "|0> + " << std::sqrt(dmatr.second) << "|1>" << std::endl; + +// ====== generate Bell00 state + teleport.H(iwf,owf,1); + iwf.swap(owf); + teleport.CNOT(iwf,owf,1,2); + iwf.swap(owf); +// ====== generate Bell00 state + + teleport.CNOT(iwf,owf,0,1); + iwf.swap(owf); + teleport.H(iwf,owf,0); + iwf.swap(owf); + bool M1 = teleport.measureQubit(iwf,0); + bool M2 = teleport.measureQubit(iwf,1); + if(M2){ + teleport.X(iwf,owf,2); + iwf.swap(owf); + } + if(M1){ + teleport.Z(iwf,2); + } + + std::cout << "final wavefunction:" << std::endl; + teleport.printWF(iwf); + + dmatr = teleport.tomography(iwf, 2); + std::cout << "State of the Bob's (id=2) qubit is " << std::sqrt(dmatr.first) << "|0> + " << std::sqrt(dmatr.second) << "|1>" << std::endl; +// ============== Teleportation end + +// ============== Deutsch-Jozsa begin + // convention: read off bit id=0. + std::cout << "\n\n\nDeutsch-Jozsa with " << DJN << " qubits:" << std::endl; + QSimulator DJozsa(DJN, false); + QSimulator::wf_type dj_iwf(DJozsa.maxState()+1,0); + QSimulator::wf_type dj_owf(DJozsa.maxState()+1,0); + QSimulator::state_type initial_state = 1; // 1 in the read off bit + dj_iwf[initial_state] = QSimulator::complex_type(1.0,0.0); + + for(unsigned qid = 0; qid<DJN; qid++){ + DJozsa.H(dj_iwf, dj_owf, qid); + dj_iwf.swap(dj_owf); + } + DJozsa.Uf(dj_iwf, dj_owf, &f_balanced); // for balanced function + //DJozsa.Uf(dj_iwf, dj_owf, &f_const1); // for constant function + dj_iwf.swap(dj_owf); + for(unsigned qid = 1; qid<DJN; qid++){ + DJozsa.H(dj_iwf, dj_owf, qid); + dj_iwf.swap(dj_owf); + } + bool allZeros = true; + for(unsigned qid = 1; qid<DJN; qid++){ + allZeros = allZeros && (!DJozsa.measureQubit(dj_iwf,qid)); + if(!allZeros){ + break; + } + } + if(allZeros){ + std::cout << "the function f is constant" << std::endl; + }else{ + std::cout << "the function f is balanced" << std::endl; + } +// ============== Deutsch-Jozsa end +} diff --git a/exercises/ex05_solution/qsimulator.ipynb b/exercises/ex05_solution/qsimulator.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..d157965b4e872c0ebcbfc3b2957a270e0379d200 --- /dev/null +++ b/exercises/ex05_solution/qsimulator.ipynb @@ -0,0 +1,350 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 48, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + } + ], + "source": [ + "%pylab inline --no-import-all\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.cm as cm\n", + "from numpy import sqrt, exp, sinh\n", + "from scipy import constants, special\n", + "import random \n", + "import copy" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# random number generators:\n", + "# random.randint(a, b)" + ] + }, + { + "cell_type": "code", + "execution_count": 103, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "class QSimulator:\n", + " uniform_float_dist = random.uniform(0.0, 1.0)\n", + " \n", + " def __init__(self,n, v):\n", + " self.n_ = n\n", + " self.max_State = (1<<n)-1\n", + " self.verbose = v\n", + " #uniform_int_dist = rnd.randint(1, n-1)\n", + " \n", + " def CNOT(self, iwf, cqid, qid):\n", + " owf = copy.deepcopy(iwf)\n", + " for st in range(self.max_State+1):\n", + " if np.absolute(iwf[st]) > 0.0 :\n", + " if st & (1<<cqid): # control qubit set\n", + " flippedSt = st ^ (1<<qid) # flips qubit \n", + " owf[flippedSt] = iwf[st]\n", + " owf[st] = iwf[flippedSt]\n", + " \n", + " if self.verbose:\n", + " print(\" wave function after CNOT gate:\\n\", \"\\n\")\n", + " print(owf, \"\\n\")\n", + " return owf\n", + " \n", + " def tomography(self, wf, qid):\n", + " alpha = 0\n", + " beta = 0\n", + " for st in range(self.max_State+1):\n", + " if st & (1<<qid):\n", + " beta += (np.conj(wf[st]) * wf[st]).real\n", + " else:\n", + " alpha += (np.conj(wf[st]) * wf[st]).real\n", + " return (alpha, beta)\n", + " \n", + " def H(self, iwf, qid):\n", + " owf = copy.deepcopy(iwf)\n", + " \n", + " for st in range(self.max_State+1):\n", + " if np.absolute(iwf[st]) > 0.0 :\n", + " \n", + " flippedSt = st ^ (1<<qid) # flips qubit\n", + " sign = (st & (1<<qid)) >> qid # 1 if qubitId is set (=1), 0 otherwise\n", + " sign = 2*sign-1 # 1 if qubitId is set (=1), -1 otherwise\n", + " \n", + " owf[flippedSt] = ( float(sign)*iwf[flippedSt] + iwf[st]) /sqrt(2)\n", + " owf[st] = ( -1*float(sign)*iwf[st] + iwf[flippedSt] ) /sqrt(2)\n", + " \n", + " if self.verbose:\n", + " print(\" wf after H gate:\\n\", owf, \"\\n\")\n", + " \n", + " return owf\n", + " \n", + " def measureQubit(self, iwf, qid):\n", + " probs_pair = self.tomography(iwf, qid)\n", + " rnd_num = random.uniform(0.0, 1.0)\n", + " measure1 = (rnd_num < probs_pair[1])\n", + " \n", + " if self.verbose:\n", + " print(\" The probability to measure |0>: \", probs_pair[0], \", for |1>: \", probs_pair[1], \"\\n\")\n", + " print(\"Random number:\", rnd_num, \", we measured 1: \", measure1, \"\\n\")\n", + " \n", + " for st in range(self.max_State+1):\n", + " if np.abs(iwf[st])!= 0.0:\n", + " if st & (1<<qid):\n", + " if measure1:\n", + " iwf[st] /= sqrt(probs_pair[1])\n", + " else:\n", + " iwf[st] = 0\n", + " else:\n", + " if measure1:\n", + " iwf[st] = 0\n", + " else:\n", + " iwf[st] /= sqrt(probs_pair[0])\n", + " \n", + " if self.verbose:\n", + " print(\"wavefunction after measurement:\\n\", iwf, \"\\n\")\n", + " \n", + " return (iwf,measure1)\n", + " \n", + " def X(self, iwf, qid):\n", + " owf = copy.deepcopy(iwf)\n", + " for st in range(self.max_State+1):\n", + " flipped = st ^(1<<qid)\n", + " owf[st] = iwf[flipped]\n", + " owf[flipped] = iwf[st]\n", + " return owf\n", + " \n", + " def Z(self, iwf, qid):\n", + " owf = copy.deepcopy(iwf)\n", + " for st in range(self.max_State+1):\n", + " if np.absolute(iwf[st]) > 0.0:\n", + " if st & (1<<qid):\n", + " owf *= -1\n", + " \n", + " return owf\n", + " \n", + " def Uf(self, iwf, func):\n", + " owf = copy.deepcopy(iwf)\n", + " hidden_cqubit = random.randint(1, self.n_ -1)\n", + " \n", + " for st in range(self.max_State+1):\n", + " if np.absolute(iwf[st]) > 0.0:\n", + " if func(st, hidden_cqubit):\n", + " newState = st ^ 1\n", + " owf[st] = iwf[newState]\n", + " owf[newState] = iwf[st]\n", + " if self.verbose:\n", + " print(\"wavefunction after Uf gate:\\n\", iwf, \"\\n\")\n", + " return owf\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 104, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# functions for DJ:\n", + "\n", + "def f_balanced(st, qbit):\n", + " return ( (st & (1<<qbit))>>qbit )\n", + "\n", + "def f_const0(st, qbit):\n", + " return 0\n", + "\n", + "def f_const1(st, qbit):\n", + " return 1\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 105, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "State of Alices's |psi> qubit (id=0) is 0.942809041582 |0> + 0.333333333333 |1> \n", + "\n", + "State of the Bob's (id=2) qubit is 0.942809041582 |0> + 0.333333333333 |1> \n", + "\n" + ] + } + ], + "source": [ + "# teleportation:\n", + "verbose = False\n", + "n = 3\n", + "teleport = QSimulator(n, verbose)\n", + "iwf = np.array( [0.0] * (teleport.max_State+1), dtype = complex)\n", + "\n", + "iwf[0] = 2*sqrt(2)/3\n", + "iwf[1] = 1/3.\n", + "\n", + "#print(\"Initial wave function:\",iwf, \"\\n\")\n", + "print(\"State of Alices's |psi> qubit (id=0) is\", sqrt(teleport.tomography(iwf,0)[0]) ,\"|0> + \",sqrt(teleport.tomography(iwf,0)[1]),\"|1> \\n\" )\n", + "\n", + "# generate Bell state:\n", + "owf = teleport.H(iwf,1)\n", + "iwf, owf = copy.deepcopy(owf), copy.deepcopy(iwf) # swap\n", + "owf = teleport.CNOT(iwf, 1,2)\n", + "iwf, owf = copy.deepcopy(owf), copy.deepcopy(iwf) # swap\n", + "\n", + "\n", + "# do teleportation:\n", + "owf = teleport.CNOT(iwf, 0,1)\n", + "iwf, owf = copy.deepcopy(owf), copy.deepcopy(iwf) # swap\n", + "owf = teleport.H(iwf,0)\n", + "iwf, owf = copy.deepcopy(owf), copy.deepcopy(iwf) # swap\n", + "\n", + "\n", + "#measure:\n", + "M0 = teleport.measureQubit(iwf,0)\n", + "M1 = teleport.measureQubit(iwf,1)\n", + "\n", + "if M0[1]:\n", + " owf = teleport.Z(iwf,2)\n", + " iwf, owf = copy.deepcopy(owf), copy.deepcopy(iwf) # swap\n", + "if M1[1]:\n", + " owf= teleport.X(iwf,2)\n", + " iwf, owf = copy.deepcopy(owf), copy.deepcopy(iwf) # swap\n", + " \n", + " \n", + "prob_pair = teleport.tomography(iwf,2)\n", + "\n", + "print(\"State of the Bob's (id=2) qubit is \", sqrt(prob_pair[0]), \"|0> + \" , sqrt(prob_pair[1]), \"|1> \\n\")" + ] + }, + { + "cell_type": "code", + "execution_count": 129, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "before UF [ 0.5+0.j -0.5+0.j 0.5+0.j -0.5+0.j]\n", + "after Uf [-0.5+0.j 0.5+0.j -0.5+0.j 0.5+0.j]\n", + "f is constant\n" + ] + } + ], + "source": [ + "# Deutsch-Jozsa:\n", + "\n", + "DJN = 2 # DJN >= 2\n", + "f = f_const1 # f_balanced, f_const0, f_const1\n", + "\n", + "DJ = QSimulator(DJN, False)\n", + "\n", + "dj_iwf = np.array( [0.0] * (DJ.max_State+1), dtype = complex)\n", + "initial_state = 1\n", + "dj_iwf[initial_state] = complex(1,0)\n", + "\n", + "for qid in range(DJN):\n", + " dj_owf = DJ.H(dj_iwf, qid)\n", + " dj_iwf, dj_owf = copy.deepcopy(dj_owf), copy.deepcopy(dj_iwf) # swap\n", + " \n", + "print(\"before UF\", dj_iwf)\n", + " \n", + "dj_owf = DJ.Uf(dj_iwf, f)\n", + "dj_iwf, dj_owf = copy.deepcopy(dj_owf), copy.deepcopy(dj_iwf) # swap\n", + "\n", + "print(\"after Uf\", dj_iwf)\n", + "\n", + "for qid in range(1,DJN):\n", + " dj_owf = DJ.H(dj_iwf, qid)\n", + " dj_iwf, dj_owf = copy.deepcopy(dj_owf), copy.deepcopy(dj_iwf) # swap\n", + "\n", + " \n", + "allZeros = True\n", + "\n", + "for qid in range(1,DJN):\n", + " allZeros = allZeros and (not DJ.measureQubit(dj_iwf,qid)[1])\n", + " if not allZeros:\n", + " break\n", + " \n", + "#print( \"f is constant\") if allZeros else print(\" f is balanced \")\n", + "if allZeros:\n", + " print(\"f is constant\")\n", + "else:\n", + " print(\" f is balanced \")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.4.4" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +}