qgtypes.cpp

/***************************************************************************
 *   QGame++                                                               *
 *   =======                                                               *
 *   A Quantum Gate And Measurement Emulator.                              *
 *                                                                         *
 *   Copyright (C) 2003/04 by Manuel Nickschas                             *
 *   <qgame-devel@nickschas.de>                                            *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "qgame++.h"
#include "qgtypes.h"
#include "qgerror.h"

using namespace qgame;
using namespace std;

/***************************************************************************
 * Class Complex                                                           *
 ***************************************************************************/

Complex::Complex(double r, double i) {
  re = r; im = i;
}

double Complex::Re() {
  return re;
}

double Complex::Im() {
  return im;
}

Complex Complex::operator+(const Complex &c) {
  double r = re + c.re;
  double i = im + c.im;
  return Complex(r,i);
}

Complex Complex::operator+=(const Complex &c) {
  re += c.re;
  im += c.im;
  return *this;
}

Complex Complex::operator-(const Complex &c) {
  double r = re - c.re;
  double i = im - c.im;
  return Complex(r,i);
}

Complex Complex::operator*(const Complex &c) {
  double r = re*c.re - im*c.im;
  double i = re*c.im + im*c.re;
  return Complex(r,i);
}

Complex Complex::operator*=(const Complex &c) {
  double r = re*c.re - im*c.im;
  double i = re*c.im + im*c.re;
  re = r; im = i;
  return *this;
}

std::ostream& qgame::operator<<(std::ostream &o, const Complex &c) {
  o << '(' << c.re << ',' << c.im << ')';
  return o;
}

double Complex::sqabs() {
  return re*re + im*im;
}

/* (namespace) global function for convenience */
double qgame::sqabs(const Complex &c) {
  return c.re*c.re + c.im*c.im;
}

/***************************************************************************
 * Class QubitList                                                         *
 ***************************************************************************/

unsigned int QubitList::size() const {
  return qbl.size();
}

void QubitList::append(int x) {
  qbl.push_back(x);
}

int QubitList::operator[](const unsigned int i) const throw(Error) {
  if(i>=qbl.size()) throw Error(E_OUTOFBOUNDS);
  return qbl[i];
}

ostream & qgame::operator<<(ostream &o, const QubitList &l) {
  o << "{ ";
  for(int i=l.size()-1; i>=0; i--) o << l[i] << " ";
  o << "}";
}

/***************************************************************************
 * Class QuSubReg                                                          *
 ***************************************************************************/

QuSubReg::QuSubReg(int s) {
  num = s;
  sub = vector<Complex*>(1<<s, (Complex*)NULL);
  //for(int i=0; i<1<<s; i++) { sub.push_back(NULL); }
}

QuSubReg::~QuSubReg() {

}

unsigned int QuSubReg::numAmps() const {
  return 1<<num;
}

unsigned int QuSubReg::numQb() const {
  return num;
}

Complex & QuSubReg::operator[](const unsigned int i) throw(Error) {
  if(i>=sub.size()) throw Error(E_OUTOFBOUNDS);
  return *sub[i];
}

/***************************************************************************
 * Class Qureg                                                             *
 ***************************************************************************/

Qureg::Qureg(int s) {
  num = s;
  amp = vector<Complex>(1<<s, Complex(0,0));
  amp[0] = Complex(1,0);  // |00...00>
}

vector<QuSubReg> Qureg::getSubRegs(const QubitList &qubits) {
  vector<QuSubReg> result;
  int numqb = qubits.size();
  int numot = num - numqb;
  QubitList otbits; // remaining qubits
  for(int i=0; i<num; i++) {
    bool flg = false;
    for(int j=0; j<numqb; j++) {
      if(qubits[j]==i) { flg = true; break; }
    }
    if(!flg) otbits.append(i);
  }
  for(int sr=0; sr<(1<<numot); sr++) {
    int index2 = 0; int dum = sr;
    for(int i = 0; i<numot; i++) {
      if(dum%2) index2 += (1<<otbits[i]);
      dum = dum>>1;
    }
    QuSubReg sreg(numqb);
    for(int sqb=0; sqb<(1<<numqb); sqb++) {
      int index = index2; dum = sqb;
      for(int i = 0; i<numqb; i++) {
        if(dum%2) index += (1<<qubits[i]);
        dum = dum>>1;
      }
      sreg.sub[sqb] = &amp[index];
    }
    result.push_back(sreg);
  }
  return result;
}

/* TODO: this needs some serious cleaning... */
vector<double> Qureg::getProbabilities(const QubitList &qubits) {
  vector<double> result = vector<double>(1<<qubits.size(),0);
  vector<QuSubReg> subregs = getSubRegs(qubits);
  double t = 0; // DEBUG // 
  //cerr << "\nabs: ";
  for(int s = 0; s < subregs.size(); s++) {
    for(int a = 0; a < subregs[s].numAmps(); a++) {
      result[a] += sqabs(subregs[s][a]);
      //cerr << sqabs(subregs[s][a])<< " ";
      t += sqabs(subregs[s][a]); // DEBUG
    }
  }
  /* sanity check works only for normalized register! */
  if(fabs(1-t)>.0000001) cerr << "\n\nWARNING: Non-unitary transformation? p="<<t;
  //if(t>1.000000001) cerr << "\n\nWARNING: Something wrong with probs in register! (p="<<t<<")";
  //cerr << "\n probs: "; for(int i=0; i< result.size(); i++) cerr << result[i]<<" "; 
  return result;
}

/* return false if the system completely collapses! */
void Qureg::collapse(const QubitList &qubits, int value) {
  vector<QuSubReg> subregs = getSubRegs(qubits);
  for(int s = 0; s < subregs.size(); s++) {
    for(int a = 0; a < subregs[s].numAmps(); a++) {
      if(a != value) subregs[s][a] = 0;
    }
  }
  normalize();
}

void Qureg::normalize() {
  double sum = 0;
  for(int i=0; i<amp.size(); i++) sum += sqabs(amp[i]);
  if(sum == 0) cerr << "\nWARNING: Cannot happen!"; 
  sum = 1.0/sqrt(sum);
  for(int i=0; i<amp.size(); i++) amp[i]*=sum; 
  //return true;
}

unsigned int Qureg::size() const {
  return num;
}

Complex & Qureg::operator[](const unsigned int i) throw(Error) {
  if(i>=amp.size()) throw Error(E_OUTOFBOUNDS);
  return amp[i];
};

#if 0
/* old, more compact representation */
std::ostream& qgame::operator<<(std::ostream &o, const Qureg &r) {
  o << "[ ";
  for(int i=0; i<r.amp.size();i++) o << r.amp[i] << ' ';
  o << "]";
  return o;
}
#else
/* this prints everything */
std::ostream& qgame::operator<<(std::ostream &o, Qureg r) {
  for(int i = 0; i<r.amp.size(); i++) {
    o << "\n   |";
    int dum = i;
    for(int j = r.num-1; j>=0; j--) {
      if(dum >= 1<<j) { dum -= 1<<j; o << "1"; }
      else o << "0";
    }
    o << ">: " << r[i] << " (p = " << sqabs(r.amp[i]) << ")";
  }
  return o;
}

#endif

/***************************************************************************
 * TestCase                                                                *
 ***************************************************************************/
 
/* parse case and set desiredRes to -1 if invalid */
TestCase::TestCase(const string &casestr) {
  char buf[256];
  desiredRes = -1;
  if(casestr.length()<3) return;
  
  int k;
  for(k=0; k<casestr.length()-2; k++) {
    if(casestr[k]=='-') break;
    if(casestr[k]!='0' && casestr[k]!='1') return;
    buf[k] = casestr[k];
  }
  buf[k] = '\0';
  if(casestr[k]!='-') return;
  char *endpnt;
  desiredRes = strtol(casestr.c_str()+k+1,&endpnt,10);
  if(*endpnt != '\0') desiredRes = -1;
  oracleTT = buf; 
}

---