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;
}