    Next: util.h Up: Code for my Simulation Previous: qureg.C   Contents

## shor.C

```#include <iostream>
#include <math.h>
#include <stdlib.h>
#include <time.h>

#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;
}
```    Next: util.h Up: Code for my Simulation Previous: qureg.C   Contents
Matthew Hayward - Quantum Computing and Shor's Algorithm GitHub Repository