shorsim/0000755000175000017500000000000012455435401012077 5ustar viblioviblioshorsim/complex.C0000644000175000017500000000267012455427735013672 0ustar viblioviblio#include
#include "complex.h"
//Complex constructor, initializes to 0 + i0.
Complex::Complex(): real( 0 ), imaginary( 0 ) {}
//Argument constructor.
Complex::Complex( double r, double i ): real( r ), imaginary( i ) {}
//Complex destructor.
Complex::~Complex() {}
//Overloaded = operator.
Complex & Complex::operator=(const Complex &c) {
if (&c != this) {
real = c.Real();
imaginary = c.Imaginary();
return *this;
}
}
//Overloaded + operator.
Complex & Complex::operator+(const Complex &c) {
real += c.Real();
imaginary += c.Imaginary();
return *this;
}
//Overloaded * operator.
Complex & Complex::operator*(const Complex &c) {
real = real * c.Real() - imaginary * c.Imaginary();
imaginary = real * c.Imaginary() + imaginary * c.Real();
return *this;
}
//Overloaded == operator. Small error tolerances.
bool Complex::operator==(const Complex &c) const {
//This is to take care of round off errors.
if (fabs(c.Real() - real) > pow(10,-14)) {
return false;
}
if (fabs(c.Imaginary()- imaginary) > pow(10,-14)) {
return false;
}
return true;
}
//Sets private data members.
void Complex::Set(double new_real, double new_imaginary) {
real = new_real;
imaginary = new_imaginary;
return;
}
//Returns the real part of the complex number.
double Complex::Real() const {
return real;
}
//Returns the imaginary part of the complex number.
double Complex::Imaginary() const {
return imaginary;
}
shorsim/qubit.C0000644000175000017500000000714312455427616013345 0ustar viblioviblio#include
#include
#include
#include "complex.h"
using namespace std;
class Qubit {
public:
Qubit(); //Default constructor.
virtual ~Qubit(); //Default destructor.
int Measure(); //Returns zero_state = 0 or one_state = 1 in accordance
//with the probabilities of zero_state and one_state.
void Dump() const; //Prints our zero_state, and one_state, without
//disturbing anything, this operation has no physical
//realization, it is only for information and debugging.
//It should never be used in an algorithm for
//information.
void SetState(const Complex& zero_prob, const Complex& one_prob);
// Takes two complex numbers and sets the states to
//those values.
void SetAverage(); //Sets the state to 2^(1/2) zero_state, 2^(1/2)
//one_state. No imaginary/phase component.
double MCC(int state) const; //Multiply the zero or one state by its complex
//conjugate, and return the value. This value
//will always be a real number, with no imaginary
//component.
private:
Complex zero_state;
//The probability of finding the Qubit in the zero or all imaginary
//state. I currently use only the real portion.
Complex one_state;
//The probability of finding the Qubit in the one or all real state.
//I currently use only the real portion.
//|zero_state|^2 + |one_state|^2 should always be 1.
//This notation means z_s * ComplexConjugate(z_s) + o_s *
//ComplexConjugate(o_s) = 1.
};
//Qubit constructor, initially sets things in the zero state.
Qubit::Qubit() {
zero_state = Complex(1,0);
one_state = Complex(0,0);
srand(time(NULL));
}
//Returns _state * ComplexConjugate(_state).
double Qubit::MCC(int state) const {
if (state == 0) {
return (pow(zero_state.Real(), 2) + pow(zero_state.Imaginary(), 2));
} else {
return (pow(one_state.Real(), 2) + pow(one_state.Imaginary(), 2));
}
}
//Measurement operator. Destructively collapses superpositions.
int Qubit::Measure() {
double rand1 = rand()/(double)RAND_MAX;
//This assumes that the sum of the squares of the amplitudes are = 1
if (MCC(0) > rand1) {
zero_state.Set(1,0);
one_state.Set(0,0);
return 0;
} else {
zero_state.Set(0,0);
one_state.Set(1,0);
return 1;
}
}
//Outputs state info for our qubit. For debugging purposes.
void Qubit::Dump() const{
cout << "Amplitude of the zero state is "
<< zero_state.Real() << " +i" << zero_state.Imaginary() << endl
<< flush;
cout << "Amplitude of the one state is "
<< one_state.Real() << " +i" << one_state.Imaginary() << endl
<< flush;
}
//Sets the zero and one states to arbitrary amplitudes. Outputs
//an error message if the two values MCC'ed != 1 + 0i.
void Qubit::SetState(const Complex& zero_prob, const Complex& one_prob) {
zero_state = zero_prob;
one_state = one_prob;
//Determine if |zero_state|^2 + |one_state|^2 == 1, if not we
//are not in a valid state.
double probab;
probab = MCC(0) + MCC(1);
if (fabs(probab - 1) > pow(10, -10)) {
//This funny expression
//allows us some rounding errors.
cout << "Warning, total probability for in SetState is different from 1." << endl << flush;
}
}
//Sets the qubit 1/2 way between the 0 state and the 1 state. No phase.
void Qubit::SetAverage() {
zero_state.Set(pow(2, -.5), 0);
one_state.Set(pow(2, -.5), 0);
}
shorsim/qureg.h0000644000175000017500000000465212455430160013377 0ustar viblioviblio#include "complex.h"
#ifndef QUREG_H
#define QUREG_H
using namespace std;
class QuReg {
public:
//Default constructor. Size is the size in bits of our register.
//In our implementation of Shor's algorithm we will need size bits
//to represent our value for "q" which is a number we have chosen
//with small prime factors which is between 2n^2 and 3n^2 inclusive
//where n is the number we are trying to factor. We envision our
//the description of our register of size "S" as 2^S complex
//numbers, representing the probability of finding the register on
//one of or 2^S base states. Thus we use an array of size 2^S, of
//Complex numbers. Thus if the size of our register is 3 bits
//array[7] is the probability amplitude of the state |1,1,1>, and
//array[7] * Complex Conjugate(array[7]) = probability of choosing
//that state. We use normalized state vectors thought the
//simulation, thus the sum of all possible states times their
//complex conjugates is = 1.
QuReg(unsigned long long int size);
QuReg(); //Default Constructor
QuReg(const QuReg &); //Copy constructor
virtual ~QuReg(); //Default destructor.
//Measures our quantum register, and returns the decimal
//interpretation of the bit string measured.
unsigned long long int DecMeasure();
//Dumps all the information about the quantum register. This
//function would not be possible for an actual quantum register, it
//is only there for debugging. When verbose != 0 we return every
//value, when verbose = 0 we return only probability amplitudes
//which differ from 0.
void Dump(int verbose) const;
//Sets state of the qubits using the arrays of complex amplitudes.
void SetState(Complex *new_state);
//Sets the state to an equal superposition of all possible states
//between 0 and number inclusive.
void SetAverage(unsigned long long int number);
//Normalize the state amplitudes.
void Norm();
//Get the probability of a given state. This is used in the
//discrete Fourier transformation. In a real quantum computer such
//an operation would not be possible, on the flip side, it would
//also not be necessary as you could simply build a DFT gate, and
//run your superposition through it to get the right answer.
Complex GetProb(unsigned long long int state) const;
//Return the size of the register.
int Size() const;
private:
unsigned long long int reg_size;
Complex *State;
};
#endif
shorsim/util.h0000644000175000017500000000344012455430100013215 0ustar viblioviblio#include "qureg.h"
#ifndef UTIL_H
#define UTIL_H
//This function takes an integer input and returns 1 if it is a prime
//number, and 0 otherwise.
unsigned long long int TestPrime(unsigned long long int n);
//This function takes an integer input and returns 1 if it is equal to
//a prime number raised to an integer power, and 0 otherwise.
unsigned long long int TestPrimePower(unsigned long long int n);
//This function computes the greatest common denominator of two integers.
//Since the modulus of a number mod 0 is not defined, we return a -1 as
//an error code if we ever would try to take the modulus of something and
//zero.
unsigned long long int GCD(unsigned long long int a, unsigned long long int b);
//This function takes and integer argument, and returns the size in bits
//needed to represent that integer.
unsigned long long int RegSize(unsigned long long int a);
//q is the power of two such that n^2 <= q < 2n^2.
unsigned long long int GetQ(unsigned long long int n);
//This function takes three integers, x, a, and n, and returns x^a mod
//n. This algorithm is known as the "Russian peasant method," I
//believe, and avoids overflow by never calculating x^a directly.
unsigned long long int modexp(unsigned long long int x, unsigned long long int a, unsigned long long int n);
// This function finds the denominator q of the best rational
// denominator q for approximating p / q for c with q < qmax.
unsigned long long int denominator(double c, unsigned long long int qmax);
//This function takes two integer arguments and returns the greater of
//the two.
unsigned long long int max(unsigned long long int a, unsigned long long int b);
//This function computes the discrete Fourier transformation on a register's
//0 -> q - 1 entries.
void DFT(QuReg * reg, unsigned long long int q);
#endif
shorsim/qureg.C0000644000175000017500000000706112455430125013330 0ustar viblioviblio#include
#include
#include
#include
#include "complex.h"
#include "qureg.h"
using namespace std;
QuReg::QuReg(): reg_size( 0 ), State( NULL ) {}
//Constructor.
QuReg::QuReg(unsigned long long int size) {
reg_size = size;
State = new Complex[(unsigned long long int)pow(2, reg_size)];
srand(time(NULL));
}
//Copy Constructor
QuReg::QuReg(const QuReg & old) {
reg_size = old.reg_size;
unsigned long long int reg_length = (unsigned long long int) pow(2, reg_size);
State = new Complex[reg_length];
for (unsigned int i = 0 ; i < reg_length ; i++) {
State[i] = old.State[i];
}
}
//Destructor.
QuReg::~QuReg() {
if ( State ) {
delete [] State;
}
}
//Return the probability amplitude of the state'th state.
Complex QuReg::GetProb(unsigned long long int state) const {
if (state >= pow(2, reg_size)) {
cout << "Invalid state index " << state << " requested!"
<< endl << flush;
throw -1;
} else {
return(State[state]);
}
}
//Normalize the probability amplitude, this ensures that the sum of
//the sum of the squares of all the real and imaginary components is
//equal to one.
void QuReg::Norm() {
double b = 0;
double f, g;
for (unsigned long long int i = 0; i < pow(2, reg_size) ; i++) {
b += pow(State[i].Real(), 2) + pow(State[i].Imaginary(), 2);
}
b = pow(b, -.5);
for (unsigned long long int i = 0; i < pow(2, reg_size) ; i++) {
f = State[i].Real() * b;
g = State[i].Imaginary() * b;
State[i].Set(f, g);
}
}
//Returns the size of the register.
int QuReg::Size() const {
return reg_size;
}
//Measure a state, and return the decimal value measured. Collapse
//the state so that the probability of measuring the measured value in
//the future is 1, and the probability of measuring any other state is
//0.
unsigned long long int QuReg::DecMeasure() {
int done = 0;
int DecVal = -1; //-1 is an error, we did not measure anything.
double rand1, a, b;
rand1 = rand()/(double)RAND_MAX;
a = b = 0;
for (unsigned long long int i = 0 ; i < pow(2, reg_size) ;i++) {
if (!done ){
b += pow(State[i].Real(), 2) + pow(State[i].Imaginary(), 2);
if (b > rand1 && rand1 > a) {
//We have just measured the i state.
for (unsigned long long int j = 0; j < pow(2, reg_size) ; j++) {
State[j].Set(0,0);
}
State[i].Set(1,0);
DecVal = i;
done = 1;
}
a += pow(State[i].Real(), 2) + pow(State[i].Imaginary(), 2);
}
}
return DecVal;
}
//For debugging, output information about the register.
void QuReg::Dump(int verbose) const {
for (unsigned long long int i = 0 ; i < pow(2, reg_size) ; i++) {
if (verbose || fabs(State[i].Real()) > pow(10,-14)
|| fabs(State[i].Imaginary()) > pow(10,-14)) {
cout << "State " << i << " has probability amplitude "
<< State[i].Real() << " + i" << State[i].Imaginary()
<< endl << flush;
}
}
}
//Set the states to those given in the new_state array.
void QuReg::SetState(Complex *new_state) {
//Beware, this function will cause segfaults if new_state is too
//small.
for (unsigned long long int i = 0 ; i < pow(2, reg_size) ; i++) {
State[i].Set(new_state[i].Real(), new_state[i].Imaginary());
}
}
//Set the State to an equal superposition of the integers 0 -> number
//- 1
void QuReg::SetAverage(unsigned long long int number) {
if (number >= pow(2, reg_size)) {
cout << "Error, initializing past end of array in qureg::SetAverage.\n";
} else {
double prob = pow(number, -.5);
for (unsigned long long int i = 0 ; i <= number ; i++) {
State[i].Set(prob, 0);
}
}
}
shorsim/README.txt0000664000175000017500000000102112455421313013566 0ustar viblioviblio1. The Shor's algorithm simulator can be built with:
make all
Which will produce the shor-simulator binary, which can be run with:
./shor-simulator
2. The simulator is defined in these files:
complex.h, complex.C - Simple complex number implementation.
qureg.h, qureg.C - Simple quantum register implementation.
util.h, util.C - Utility functions used by the simulator.
shor.C - Main logic of the simulator
3. There is also included a simple qubit.C program which defines a
single qubit, it is not used bit the simulation.
shorsim/complex.h0000644000175000017500000000133112455416671013725 0ustar viblioviblio#ifndef COMPLEX_H
#define COMPLEX_H
class Complex {
public:
Complex(); //Default constructor
Complex( double, double ); // Argument constructor
virtual ~Complex(); //Default destructor.
void Set(double new_real, double new_imaginary); //Set data members.
double Real() const; //Return the real part.
double Imaginary() const; //Return the imaginary part.
Complex & operator+(const Complex&); //Overloaded + operator
Complex & operator*(const Complex&); //Overloaded * operator
Complex & operator=(const Complex&); //Overloaded = operator
bool operator==(const Complex&) const; //Overloaded == operator
private:
double real;
double imaginary;
};
#endif
shorsim/util.C0000644000175000017500000001171712455427711013174 0ustar viblioviblio#include
#include
#include "qureg.h"
using namespace std;
//This function takes an integer input and returns 1 if it is a prime
//number, and 0 otherwise.
//
// Not optimized at all.
unsigned long long int TestPrime(unsigned long long int n) {
unsigned long long int i;
for (i = 2 ; i <= floor(sqrt(n)) ; i++) {
if (n % i == 0) {
return(0);
}
}
return(1);
}
//This function takes an integer input and returns 1 if it is equal to
//a prime number raised to an integer power, and 0 otherwise.
unsigned long long int TestPrimePower(unsigned long long int n) {
unsigned long long int i,j;
j = 0;
i = 2;
while ((i <= floor(pow(n, .5))) && (j == 0)) {
if((n % i) == 0) {
j = i;
}
i++;
}
for (unsigned long long int i = 2 ; i <= (floor(log(n) / log(j)) + 1) ; i++) {
if(pow(j , i) == n) {
return(1);
}
}
return(0);
}
//This function computes the greatest common denominator of two integers.
//Since the modulus of a number mod 0 is not defined, we return a -1 as
//an error code if we ever would try to take the modulus of something and
//zero.
unsigned long long int GCD(unsigned long long int a, unsigned long long int b) {
unsigned long long int d;
if (b != 0) {
while (a % b != 0) {
d = a % b;
a = b;
b = d;
}
} else {
return -1;
}
return(b);
}
//This function takes and integer argument, and returns the size in bits
//needed to represent that integer.
unsigned long long int RegSize(unsigned long long int a) {
unsigned long long int size = 0;
while(a != 0) {
a = a>>1;
size++;
}
return(size);
}
//q is the power of two such that n^2 <= q < 2n^2.
unsigned long long int GetQ(unsigned long long int n) {
unsigned long long int power = 8; //256 is the smallest q ever is.
while (pow(2,power) < pow(n,2)) {
power = power + 1;
}
return 1 << power;
}
//This function takes three integers, x, a, and n, and returns x^a mod
//n. This algorithm is known as the "Russian peasant method," I
//believe, and avoids overflow by never calculating x^a directly.
unsigned long long int modexp(unsigned long long int x, unsigned long long int a, unsigned long long int n) {
unsigned long long int value = 1;
unsigned long long int tmp;
tmp = x % n;
while (a > 0) {
if (a & 1) {
value = (value * tmp) % n;
}
tmp = tmp * tmp % n;
a = a>>1;
}
return value;
}
// This function finds the denominator q of the best rational
// denominator q for approximating p / q for c with q < qmax.
unsigned long long int denominator(double c, unsigned long long int qmax) {
double y = c;
double z;
unsigned long long int q0 = 0;
unsigned long long int q1 = 1;
unsigned long long int q2 = 0;
while (1) {
z = y - floor(y);
if (z < 0.5 / pow(qmax,2)) {
return(q1);
}
if (z != 0) {
//Can't divide by 0.
y = 1 / z;
} else {
//Warning this is broken if q1 == 0, but that should never happen.
return(q1);
}
q2 = (unsigned long long int)floor(y) * q1 + q0;
if (q2 >= qmax) {
return(q1);
}
q0 = q1;
q1 = q2;
}
}
//This function takes two integer arguments and returns the greater of
//the two.
unsigned long long int max(unsigned long long int a, unsigned long long int b) {
if (a > b) {
return(a);
}
return(b);
}
//This function computes the discrete Fourier transformation on a register's
//0 -> q - 1 entries.
void DFT(QuReg * reg, unsigned long long int q) {
//The Fourier transform maps functions in the time domain to
//functions in the frequency domain. Frequency is 1/period, thus
//this Fourier transform will take our periodic register, and peak it
//at multiples of the inverse period. Our Fourier transformation on
//the state a takes it to the state: q^(-.5) * Sum[c = 0 -> c = q - 1,
//c * e^(2*Pi*i*a*c / q)]. Remember, e^ix = cos x + i*sin x.
Complex * init = new Complex[q];
Complex tmpcomp( 0, 0 );
//Here we do things that a real quantum computer couldn't do, such
//as look as individual values without collapsing state. The good
//news is that in a real quantum computer you could build a gate
//which would what this out all in one step.
unsigned long long int count = 0;
double tmpreal = 0;
double tmpimag = 0;
double tmpprob = 0;
for (unsigned long long int a = 0 ; a < q ; a++) {
//This if statement helps prevent previous round off errors from
//propagating further.
if ((pow(reg->GetProb(a).Real(),2) +
pow(reg->GetProb(a).Imaginary(),2)) > pow(10,-14)) {
for (unsigned long long int c = 0 ; c < q ; c++) {
tmpcomp.Set(pow(q,-.5) * cos(2*M_PI*a*c/q),
pow(q,-.5) * sin(2*M_PI*a*c/q));
init[c] = init[c] + (reg->GetProb(a) * tmpcomp);
}
}
count++;
if (count == 100) {
cout << "Making progress in Fourier transform, "
<< 100*((double)a / (double)(q - 1)) << "% done!"
<< endl << flush;
count = 0;
}
}
reg->SetState(init);
reg->Norm();
delete [] init;
}
shorsim/shor.C0000644000175000017500000003064412455430037013165 0ustar viblioviblio#include
#include
#include
#include
#include "complex.h"
#include "qureg.h"
#include "util.h"
using namespace std;
int main() {
//Establish a random seed.
srand(time(NULL));
//Output standard greeting.
cout << "Welcome to the simulation of Shor's algorithm." << endl
<< "There are four restrictions for Shor's algorithm:" << endl
<< "1) The number to be factored (n) must be >= 15." << endl
<< "2) The number to be factored must be odd." << endl
<< "3) The number must not be prime." << endl
<< "4) The number must not be a prime power." << endl
<< endl << "There are efficient classical methods of factoring "
<< "any of the above numbers, or determining that they are prime."
<< endl << endl << "Input the number you wish to factor." << endl
<< flush;
//n is the number we are going to factor, get n.
unsigned long long int n;
cin >> n;
cout << "Step 1 starting." << endl << flush;
//Test to see if n is factorable by Shor's algorithm.
//Exit if the number is even.
if (n%2 == 0) {
cout << "\tError, the number must be odd!" << endl << flush;
exit(0);
}
//Exit if the number is prime.
if (TestPrime(n)) {
cout << "\tError, the number must not be prime!" << endl << flush;
exit(0);
}
//Prime powers are prime numbers raised to integral powers.
//Exit if the number is a prime power.
if (TestPrimePower(n)) {
cout << "\tError, the number must not be a prime power!" << endl << flush;
exit(0);
}
cout << "Step 1 complete." << endl << flush;
cout << "Step 2 starting." << endl << flush;
//Now we must figure out how big a quantum register we need for our
//input, n. We must establish a quantum register big enough to hold
//an equal superposition of all integers 0 through q - 1 where q is
//the power of two such that n^2 <= q < 2n^2.
unsigned long long int q;
cout << "\tSearching for q, the smallest power of 2 greater than or equal to n^2." << endl << flush;
q = GetQ(n);
cout << "\tFound q to be " << q << "." << endl << flush;
cout << "Step 2 complete." << endl << flush;
cout << "Step 3 starting." << endl << flush;
//Now we must pick a random integer x, coprime to n. Numbers are
//coprime when their greatest common denominator is one. One is not
//a useful number for the algorithm.
unsigned long long int x = 0;
cout << "\tSearching for x, a random integer coprime to n." << endl << flush;
x = 1+ (unsigned long long int)((n-1)*(double)rand()/(double)RAND_MAX);
while (GCD(n,x) != 1 || x == 1) {
x = 1 + (unsigned long long int)((n-1)*(double)rand()/(double)RAND_MAX);
}
cout << "\tFound x to be " << x << "." << endl << flush;
cout << "Step 3 complete." << endl << flush;
//Create the register.
cout << "Step 4 starting." << endl << flush;
QuReg * reg1 = new QuReg(RegSize(q) - 1);
cout << "\tMade register 1 with register size = " << RegSize(q) << endl
<< flush;
//This array will remember what values of q produced for x^q mod n.
//It is necessary to retain these values for use when we collapse
//register one after measuring register two. In a real quantum
//computer these registers would be entangled, and thus this extra
//bookkeeping would not be needed at all. The laws of quantum
//mechanics dictate that register one would collapse as well, and
//into a state consistent with the measured value in resister two.
unsigned long long int * modex = new unsigned long long int[q];
//This array holds the probability amplitudes of the collapsed state
//of register one, after register two has been measured it is used
//to put register one in a state consistent with that measured in
//register two.
Complex * collapse = new Complex[q];
//This is a temporary value.
Complex tmp;
//This is a new array of probability amplitudes for our second
//quantum register, that populated by the results of x^a mod n.
Complex * mdx = new Complex[(unsigned long long int)pow(2,RegSize(n))];
// This is the second register. It needs to be big enough to hold
// the superposition of numbers ranging from 0 -> n - 1.
QuReg *reg2 = new QuReg(RegSize(n));
cout << "\tCreated register 2 of size " << RegSize(n) << endl << flush;
cout << "Step 4 complete." << endl << flush;
//This is a temporary value.
unsigned long long int tmpval;
//This is a temporary value.
unsigned long long int value;
//c is some multiple lambda of q/r, where q is q in this program,
//and r is the period we are trying to find to factor n. m is the
//value we measure from register one after the Fourier
//transformation.
double c, m;
//This is used to store the denominator of the fraction p / den where
//p / den is the best approximation to c with den <= q.
unsigned long long int den;
//This is used to store the numerator of the fraction p / den where
//p / den is the best approximation to c with den <= q.
unsigned long long int p;
//The integers e, a, and b are used in the end of the program when
//we attempts to calculate the factors of n given the period it
//measured.
//Factor is the factor that we find.
unsigned long long int e, a, b, factor;
//Shor's algorithm can sometimes fail, in which case you do it
//again. The done variable is set to 0 when the algorithm has
//failed. Only try a maximum number of tries.
unsigned int done = 0;
unsigned int tries = 0;
while (!done) {
if (tries >= 5) {
cout << "\tThere have been five failures, giving up." << endl << flush;
exit(0);
}
cout << "Step 5 starting attempt: " << tries+1 << endl << flush;
//Now populate register one in an even superposition of the
//integers 0 -> q - 1.
reg1->SetAverage(q - 1);
cout << "Step 5 complete." << endl << flush;
cout << "Step 6 starting attempt: " << tries+1 << endl << flush;
//Now we preform a modular exponentiation on the superposed
//elements of reg 1. That is, perform x^a mod n, but exploiting
//quantum parallelism a quantum computer could do this in one
//step, whereas we must calculate it once for each possible
//measurable value in register one. We store the result in a new
//register, reg2, which is entangled with the first register.
//This means that when one is measured, and collapses into a base
//state, the other register must collapse into a superposition of
//states consistent with the measured value in the other.. The
//size of the result modular exponentiation will be at most n, so
//the number of bits we will need is therefore less than or equal
//to log2 of n. At this point we also maintain a array of what
//each state produced when modularly exponised, this is because
//these registers would actually be entangled in a real quantum
//computer, this information is needed when collapsing the first
//register later.
//This counter variable is used to increase our probability amplitude.
tmp.Set(1,0);
//This for loop ranges over q, and puts the value of x^a mod n in
//modex[a]. It also increases the probability amplitude of the value
//of mdx[x^a mod n] in our array of complex probabilities.
for (unsigned long long int i = 0 ; i < q ; i++) {
//We must use this version of modexp instead of c++ builtins as
//they overflow when x^i is large.
tmpval = modexp(x,i,n);
modex[i] = tmpval;
mdx[tmpval] = mdx[tmpval] + tmp;
}
//Set the state of register two to what we calculated it should be.
reg2->SetState(mdx);
//Normalize register two, so that the probability of measuring a
//state is given by summing the squares of its probability
//amplitude.
reg2->Norm();
cout << "Step 6 complete." << endl << flush;
cout << "Step 7 starting attempt: " << tries+1 << endl << flush;
//Now we measure reg2.
value = reg2->DecMeasure();
//Now we must using the information in the array modex collapse
//the state of register one into a state consistent with the value
//we measured in register two.
for (unsigned long long int i = 0 ; i < q ; i++) {
if (modex[i] == value) {
collapse[i].Set(1,0);
} else {
collapse[i].Set(0,0);
}
}
//Now we set the state of register one to be consistent with what
//we measured in state two, and normalize the probability
//amplitudes.
reg1->SetState(collapse);
reg1->Norm();
cout << "Step 7 complete." << endl << flush;
//Here we do our Fourier transformation.
cout << "Step 8 starting attempt: " << tries+1 << endl << flush;
DFT(reg1, q);
cout << "Step 8 complete." << endl << flush;
cout << "Step 9 starting attempt: " << tries+1 << endl << flush;
//Next we measure register one, due to the Fourier transform the
//number we measure, m will be some multiple of lambda/r, where
//lambda is an integer and r is the desired period.
m = reg1->DecMeasure();
cout << "\tValue of m measured as: " << m << endl << flush;
cout << "Step 9 complete." << endl << flush;
//If nothing goes wrong from here on out we are done.
done = 1;
//If we measured zero, we have gained no new information about the
//period, we must try again.
if (m == 0) {
cout << "\tMeasured, 0 this trial a failure!" << endl << flush;
done = 0;
}
//The DecMeasure subroutine will return -1 as an error code, due
//to rounding errors it will occasionally fail to measure a state.
if (m == -1) {
cout << "\tWe failed to measure anything, this trial a failure!" << endl << flush;
done = 0;
}
//If nothing has gone wrong, try to determine the period of our
//function, and get factors of n.
if (done) {
//Now c =~ lambda / r for some integer lambda. Borrowed with
//modifications from Berhnard Ohpner.
c = (double)m / (double)q;
cout << "Steps 10 and 11 starting attempt: " << tries+1 << endl << flush;
//Calculate the denominator of the best rational approximation
//to c with den < q. Since c is lambda / r for some integer
//lambda, this will provide us with our guess for r, our period.
den = denominator(c, q);
//Calculate the numerator from the denominator.
p = (unsigned long long int)floor(den * c + 0.5);
//Give user information.
cout << "\tMeasured m: " << m << ", rational approximation for m/q=" << c << " is: "
<< p << " / " << den << endl << flush;
//The denominator is our period, and an odd period is not
//useful as a result of Shor's algorithm. If the denominator
//times two is still less than q we can use that.
if (den % 2 == 1 && 2 * den < q ){
cout << "\tOdd candidate for r found, expanding by 2\n";
p = 2 * p;
den = 2 * den;
}
//Initialize helper variables.
e = a = b = factor = 0;
// Failed if odd denominator.
if (den % 2 == 1) {
cout << "\tOdd period found. This trial failed."
<< " \t Trying again." << endl << flush;
done = 0;
} else {
//Calculate candidates for possible common factors with n.
cout << "\tCandidate period is " << den << endl << flush;
e = modexp(x, den / 2, n);
a = (e + 1) % n;
b = (e + n - 1) % n;
cout << "\t" << x << "^" << den / 2 << " + 1 mod " << n << " = " << a
<< "," << endl
<< "\t" << x << "^" << den / 2 << " - 1 mod " << n << " = " << b
<< endl << flush;
factor = max(GCD(n,a),GCD(n,b));
}
}
//GCD will return a -1 if it tried to calculate the GCD of two
//numbers where at some point it tries to take the modulus of a
//number and 0.
if (factor == -1) {
cout << "\tError, tried to calculate n mod 0 for some n. Trying again."
<< endl << flush;
done = 0;
}
if ((factor == n || factor == 1) && done == 1) {
cout << "\tFound trivial factors 1 and " << n
<< ". Trying again." << endl << flush;
done = 0;
}
//If nothing else has gone wrong, and we got a factor we are
//finished. Otherwise start over.
if (factor != 0 && done == 1) {
cout << "\t" << n << " = " << factor << " * " << n / factor << endl << flush;
} else if (done == 1) {
cout << "\tFound factor to be 0, error. Trying again." << endl
<< flush;
done = 0;
}
cout << "Steps 10 and 11 complete." << endl << flush;
tries++;
}
delete reg1;
delete reg2;
delete [] modex;
delete [] collapse;
delete [] mdx;
return 1;
}
shorsim/Makefile0000664000175000017500000000031512455427324013545 0ustar viblioviblioFILES=\
complex.C\
qureg.C\
util.C\
shor.C
OBJS=$(FILES:.C=.o)
BINARY=shor-simulator
%.o : %.cpp
g++ -O3 -c $< -o $@
all : $(OBJS)
g++ -O3 -o $(BINARY) $(OBJS)
clean :
rm -f $(OBJS) $(BINARY)