www.pudn.com > shor.rar > shorDlg.cpp, change:2009-02-23,size:37018b


// shorDlg.cpp : implementation file 
// 
 
#include "stdafx.h" 
#include "shor.h" 
#include "shorDlg.h" 
 
#ifdef _DEBUG 
#define new DEBUG_NEW 
#undef THIS_FILE 
static char THIS_FILE[] = __FILE__; 
#endif 
#include <iostream.h> 
#include <math.h> 
#include <stdlib.h> 
#include <time.h> 
#define PI 3.14159265359 
CString s; 
 
 
 
 
class Complex { 
public: 
	Complex();                  //Default constructor 
	~Complex();                 //Default destructor. 
	void Set(double new_real, double new_imaginary); //Set data members. 
	double Real();              //Return the real part. 
	double Imaginary();         //Return the imaginary part. 
	Complex operator+(Complex); //Overloaded + operator 
	Complex operator*(Complex); //Overloaded * operator 
	Complex operator=(Complex); //Overloaded = operator 
	int operator==(Complex);    //Overloaded == operator 
private: 
	double real; 
	double imaginary; 
}; 
 
//Complex constructor, initialises to 0 + i0. 
Complex::Complex() { 
	real = 0; 
	imaginary = 0; 
} 
 
//Complex destructor. 
Complex::~Complex() {}; 
 
//Overloaded = operator. 
Complex Complex::operator=(Complex c) { 
	if (&c != this) { 
		real = c.Real(); 
		imaginary = c.Imaginary(); 
		return *this; 
	}else 
		exit(0); 
} 
 
//Overloaded + operator. 
Complex Complex::operator+(Complex c) {  
	Complex tmp;  
	double new_real, new_imaginary;  
	new_real = real + c.Real(); 
	new_imaginary = imaginary + c.Imaginary();  
	tmp.Set(new_real,new_imaginary);  
	return tmp; 
} 
 
//Overloaded * operator. 
Complex Complex::operator*(Complex c) {  
	Complex tmp;  
	double new_real, new_imaginary;  
	new_real = real * c.Real() - imaginary * c.Imaginary(); 
	new_imaginary = real * c.Imaginary() + imaginary * c.Real(); 
	tmp.Set(new_real,new_imaginary);  
	return tmp; 
} 
 
//Overloaded == operator.  Small error tolerances. 
int Complex::operator==(Complex c) { 
	//This is to take care of round off errors. 
	if (fabs(c.Real() - real) > pow(10,-14)) { 
		return 0; 
	} 
	if (fabs(c.Imaginary()- imaginary) > pow(10,-14)) { 
		return 0; 
	} 
	return 1; 
} 
 
//Sets private data members. 
void Complex::Set(double new_real, double new_imaginary) { 
	real = new_real; 
	imaginary = new_imaginary; 
} 
 
//Returns the real part of the complex number. 
double Complex::Real() { 
	return real; 
} 
 
//Returns the imaginary part of the complex number. 
double Complex::Imaginary() { 
	return imaginary; 
} 
 
 
 
class Qubit { 
public: 
	Qubit();            //Default constructor. 
	~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();        //Prints our zero_state, and one_state, without  
	//disturbing anything, this operation has no physical 
	//realization, it is only for information and debuging.  
	//It should never be used in an algorithm for  
	//information. 
	void SetState(Complex zero_prob, 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);  //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 imarinary  
	//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.Set(1,0);  
	one_state.Set(0,0); 
	srand(time(NULL)); 
} 
 
//Returns <state>_state * ComplexConjugate(<state>_state). 
double Qubit::MCC(int state) { 
	if (state == 0) { 
		return (pow(zero_state.Real(), 2) + pow(zero_state.Imaginary(), 2)); 
	} 
	return (pow(one_state.Real(), 2) + pow(one_state.Imaginary(), 2)); 
} 
 
//Measurement operator.  Destructively collapses superpositions. 
int Qubit::Measure() { 
	double rand1; 
	rand1 = rand()/(double)RAND_MAX; 
	//Assumes that the sum of the squares of the amplitudes are = 1 
	 
	if (MCC(0) > rand1) { 
		//Should I be collapsing phases as well? 
		zero_state.Set(1,0); 
		one_state.Set(0,0); 
		return 0; 
	} else { 
		//Should I be collapsing phases as well? 
		zero_state.Set(0,0); 
		one_state.Set(1,0); 
		return 1; 
	}  
} 
 
//Outputs state info for our qubit.  For debugging purposes. 
void Qubit::Dump() { 
	cout << "Amplitude of the zero state is " 
		<< zero_state.Real() << " +i" << zero_state.Imaginary() << endl  
		<< flush; 
	cout << "Amplutide 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(Complex zero_prob, 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); 
} 
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 number, 
	//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 normalised state vectors thought the simulation, thus the 
	//sum of all possible states times their complex conjugates is = 1. 
	QuReg(int size);  
	 
	QuReg(); //Default Constructor 
	 
	QuReg(const QuReg &); //Copy constructor 
	 
	~QuReg(); //Default destructor. 
	 
	int DecMeasure(); //Measures our quantum register, and returns the 
	//decimal interpretation of the bitstring measured.   
	 
	//Dumps all the information about the quantum register.  This has no 
	//physical reality, 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); 
	 
	//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(int number); 
	 
	//Normalise 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(int state); 
	 
	//Return the size of the register. 
	int Size(); 
	 
private: 
	int reg_size; 
	Complex *State; 
}; 
 
QuReg::QuReg() { 
	reg_size = 0; 
	State = 0; 
} 
 
//Constructor. 
QuReg::QuReg(int size) { 
	reg_size = size; 
	State = new Complex[(int)pow(2, reg_size)]; 
	srand(time(NULL)); 
} 
 
//Copy Constructor 
QuReg::QuReg(const QuReg & old) { 
	reg_size = old.reg_size; 
	int reg_length = (int) pow(2, reg_size); 
	State = new Complex[reg_length]; 
	for (int i = 0 ; i < reg_length ; i++) { 
		State[i] = old.State[i]; 
	} 
} 
 
//Destructor. 
QuReg::~QuReg() { 
	delete [] State; 
} 
 
//Return the probability amplitude of the state'th state. 
Complex QuReg::GetProb(int state) { 
	if (state >= pow(2, reg_size)) { 
		cout << "You are trying to measure past the end of an array in qureg::GetProb!"  
			<< endl << flush; 
		exit(0); 
	}  
	else 
    { 
		return(State[state]); 
	} 
} 
 
//Normalise 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; 
	double f, g; 
	b = 0; 
	for (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 ( 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() { 
	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. 
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 (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 (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) { 
	for (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 (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(int number) { 
	if (number >= pow(2, reg_size)) { 
		cout << "Error, initialising past end of array in qureg::SetAverage.\n"; 
	} else { 
		double prob; 
		prob = pow(number, -.5); 
		for (int i = 0 ; i <= number ; i++) { 
			State[i].Set(prob, 0); 
		} 
	} 
} 
//This function takes an integer input and returns 1 if it is a prime 
//number, and 0 otherwise. 
int TestPrime(int n) { 
	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. 
int TestPrimePower(int n) { 
	int i,j; 
	j = 0; 
	i = 2; 
	while ((i <= floor(pow(n, .5))) && (j == 0)) { 
         if((n % i) == 0) { 
			j = i; 
		} 
		i++; 
	} 
	for ( 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. 
int GCD(int a, int b) { 
	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. 
int RegSize(int a) { 
	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. 
int GetQ(int n) { 
	int power = 8; //265 is the smallest q ever is. 
	while (pow(2,power) < pow(n,2)) { 
		power = power + 1; 
	} 
	return((int)pow(2,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 belive. 
int modexp(int x, int a, int n) { 
	int value = 1; 
	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. 
int denominator(double c, int qmax) { 
	double y = c; 
	double z; 
	int q0 = 0; 
	int q1 = 1; 
	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 = (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. 
int max1(int a, 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, 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; 
	tmpcomp.Set(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. 
	 
	int count = 0; 
	double tmpreal = 0; 
	double tmpimag = 0; 
	double tmpprob = 0; 
	for (int a = 0 ; a < q ; a++) { 
		//This if statement helps prevent previos round off errors from 
		//propogating further. 
		if ((pow(reg->GetProb(a).Real(),2) +  
			pow(reg->GetProb(a).Imaginary(),2)) > pow(10,-14)) { 
			for (int c = 0 ; c < q ; c++) { 
				tmpcomp.Set(pow(q,-.5) * cos(2*PI*a*c/q), 
					pow(q,-.5) * sin(2*PI*a*c/q)); 
				init[c] = init[c] + (reg->GetProb(a) * tmpcomp); 
			} 
		} 
		count++; 
		if (count == 100) { 
			s.Format("傅立叶变换进行中...,完成%d",100*((double)a / (double)(q - 1))); 
			count = 0; 
		} 
	} 
	reg->SetState(init); 
	reg->Norm(); 
	delete [] init; 
} 
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 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. 
	int n; 
	cin >> n; 
	 
	//Test to see if n is factorable by Shor's algorithm. 
	//Exit if the number is even. 
	if (n%2 == 0) { 
		cout << "Error, the number must be odd!" << endl << flush; 
		exit(0); 
	}  
	//Exit if the number is prime. 
	if (TestPrime(n)) { 
		cout << "Error, 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 << "Error, the number must not be a prime power!" << endl << flush; 
		exit(0); 
	} 
	 
	//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. 
	int x = 0; 
	x = 1+ (int)((n-1)*(double)rand()/(double)RAND_MAX); 
	while (GCD(n,x) != 1 || x == 1) { 
		x = 1 + (int)((n-1)*(double)rand()/(double)RAND_MAX); 
	} 
	cout << "Found x to be " << x << "." << 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. 
	int q; 
	q = GetQ(n); 
	cout << "Found q to be " << q << "." << endl << flush; 
	 
	//Create the register. 
	QuReg * reg1 = new QuReg(RegSize(q) - 1); 
	cout << "Made 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. 
	int * modex = new 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[(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 << "Created register 2 of size " << RegSize(n) << endl << flush; 
	 
	//This is a temporary value. 
	int tmpval; 
	 
	//This is a temporary value. 
	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. 
	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. 
	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. 
	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. 
	int done = 0; 
	int tries = 0; 
	while (!done) { 
		if (tries >= 5) { 
			cout << "There have been five failures, giving up." << endl << flush; 
			exit(0); 
		} 
		//Now populate register one in an even superposition of the 
		//integers 0 -> q - 1. 
		reg1->SetAverage(q - 1); 
		 
		//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 (int i = 0 ; i < q ; i++) { 
			//We must use this version of modexp instead of c++ builtins as 
			//they overflow when x^i > 2^31. 
			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); 
		 
		//Normalise register two, so that the probability of measuring a 
		//state is given by summing the squares of its probability 
		//amplitude. 
		reg2->Norm(); 
		 
		//Now we measure reg1.  
		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 ( 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 normalise the probability 
		//amplitudes. 
		reg1->SetState(collapse); 
		reg1->Norm(); 
		 
		//Here we do our Fourier transformation.   
		cout << "Begin Discrete Fourier Transformation!" << endl << flush; 
		DFT(reg1, q); 
		 
		//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(); 
		 
		//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 << "Measured, 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 << "We failed to measure anything, this trial a failure!"  
				<< "  Trying again." << 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; 
			 
			//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 = (int)floor(den * c + 0.5); 
			 
			//Give user information. 
			cout << "measured " << m <<", approximation for " << 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 << "Odd denominator, expanding by 2\n"; 
				p = 2 * p;  
				den = 2 * den; 
			} 
			 
			//Initialise helper variables. 
			e = a = b = factor = 0; 
			 
			// Failed if odd denominator. 
			if (den % 2 == 1)  {           
				cout <<  "Odd period found.  This trial failed."  
					<< "  Trying again." << endl << flush; 
				done = 0; 
			} else { 
				//Calculate candidates for possible common factors with n. 
				cout <<  "possible period is " << den << endl << flush; 
				e = modexp(x, den / 2, n);     
				a = (e + 1) % n;      
				b = (e + n - 1) % n;           
				cout << x << "^" << den / 2 << " + 1 mod " << n << " = " << a  
					<< "," << endl  
					<< x << "^" << den / 2 << " - 1 mod " << n << " = " << b  
					<< endl << flush; 
				factor = max1(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 << "Error, tried to calculate n mod 0 for some n.  Trying again." 
				<< endl << flush; 
			done = 0; 
		} 
		 
		if ((factor == n || factor == 1) && done == 1) { 
			cout << "Found 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 << n << " = " << factor << " * " << n / factor << endl << flush; 
		} else if (done == 1) { 
			cout << "Found factor to be 0, error.  Trying again." << endl  
				<< flush; 
			done = 0; 
		} 
		tries++; 
  } 
  delete reg1; 
  delete reg2; 
  delete [] modex; 
  delete [] collapse; 
  delete [] mdx; 
  return 1; 
} 
 
 
 
///////////////////////////////////////////////////////////////////////////// 
// CAboutDlg dialog used for App About 
 
class CAboutDlg : public CDialog 
{ 
public: 
	CAboutDlg(); 
	 
	// Dialog Data 
	//{{AFX_DATA(CAboutDlg) 
	enum { IDD = IDD_ABOUTBOX }; 
	//}}AFX_DATA 
	 
	// ClassWizard generated virtual function overrides 
	//{{AFX_VIRTUAL(CAboutDlg) 
protected: 
	virtual void DoDataExchange(CDataExchange* pDX);    // DDX/DDV support 
	//}}AFX_VIRTUAL 
	 
	// Implementation 
protected: 
	//{{AFX_MSG(CAboutDlg) 
	//}}AFX_MSG 
	DECLARE_MESSAGE_MAP() 
}; 
 
CAboutDlg::CAboutDlg() : CDialog(CAboutDlg::IDD) 
{ 
	//{{AFX_DATA_INIT(CAboutDlg) 
	//}}AFX_DATA_INIT 
} 
 
void CAboutDlg::DoDataExchange(CDataExchange* pDX) 
{ 
	CDialog::DoDataExchange(pDX); 
	//{{AFX_DATA_MAP(CAboutDlg) 
	//}}AFX_DATA_MAP 
} 
 
BEGIN_MESSAGE_MAP(CAboutDlg, CDialog) 
//{{AFX_MSG_MAP(CAboutDlg) 
// No message handlers 
//}}AFX_MSG_MAP 
END_MESSAGE_MAP() 
 
///////////////////////////////////////////////////////////////////////////// 
// CShorDlg dialog 
 
CShorDlg::CShorDlg(CWnd* pParent /*=NULL*/) 
: CDialog(CShorDlg::IDD, pParent) 
{ 
	//{{AFX_DATA_INIT(CShorDlg) 
	m_edit1 = 0; 
	m_edit2 = _T(""); 
	m_edit3 = _T(""); 
	m_edit4 = _T(""); 
	m_edit5 = _T(""); 
	m_edit6 = _T(""); 
	m_edit7 = _T(""); 
	m_edit8 = _T(""); 
	m_edit9 = _T(""); 
	m_edit10 = _T(""); 
	//}}AFX_DATA_INIT 
	// Note that LoadIcon does not require a subsequent DestroyIcon in Win32 
	m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME); 
} 
 
void CShorDlg::DoDataExchange(CDataExchange* pDX) 
{ 
	CDialog::DoDataExchange(pDX); 
	//{{AFX_DATA_MAP(CShorDlg) 
	DDX_Text(pDX, IDC_EDIT1, m_edit1); 
	DDX_Text(pDX, IDC_EDIT2, m_edit2); 
	DDX_Text(pDX, IDC_EDIT3, m_edit3); 
	DDX_Text(pDX, IDC_EDIT4, m_edit4); 
	DDX_Text(pDX, IDC_EDIT5, m_edit5); 
	DDX_Text(pDX, IDC_EDIT6, m_edit6); 
	DDX_Text(pDX, IDC_EDIT7, m_edit7); 
	DDX_Text(pDX, IDC_EDIT8, m_edit8); 
	DDX_Text(pDX, IDC_EDIT9, m_edit9); 
	DDX_Text(pDX, IDC_EDIT10, m_edit10); 
	//}}AFX_DATA_MAP 
} 
 
BEGIN_MESSAGE_MAP(CShorDlg, CDialog) 
//{{AFX_MSG_MAP(CShorDlg) 
ON_WM_SYSCOMMAND() 
ON_WM_PAINT() 
ON_WM_QUERYDRAGICON() 
	ON_BN_CLICKED(IDC_BUTTON1, OnButton1) 
	ON_BN_CLICKED(IDC_BUTTON2, OnButton2) 
	//}}AFX_MSG_MAP 
END_MESSAGE_MAP() 
 
///////////////////////////////////////////////////////////////////////////// 
// CShorDlg message handlers 
 
BOOL CShorDlg::OnInitDialog() 
{ 
	CDialog::OnInitDialog(); 
	 
	// Add "About..." menu item to system menu. 
	 
	// IDM_ABOUTBOX must be in the system command range. 
	ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX); 
	ASSERT(IDM_ABOUTBOX < 0xF000); 
	 
	CMenu* pSysMenu = GetSystemMenu(FALSE); 
	if (pSysMenu != NULL) 
	{ 
		CString strAboutMenu; 
		strAboutMenu.LoadString(IDS_ABOUTBOX); 
		if (!strAboutMenu.IsEmpty()) 
		{ 
			pSysMenu->AppendMenu(MF_SEPARATOR); 
			pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu); 
		} 
	} 
	 
	// Set the icon for this dialog.  The framework does this automatically 
	//  when the application's main window is not a dialog 
	SetIcon(m_hIcon, TRUE);			// Set big icon 
	SetIcon(m_hIcon, FALSE); 
	 
	// Set small icon 
 
	// TODO: Add extra initialization here 
	 
	return TRUE;  // return TRUE  unless you set the focus to a control 
} 
 
void CShorDlg::OnSysCommand(UINT nID, LPARAM lParam) 
{ 
	if ((nID & 0xFFF0) == IDM_ABOUTBOX) 
	{ 
		CAboutDlg dlgAbout; 
		dlgAbout.DoModal(); 
	} 
	else 
	{ 
		CDialog::OnSysCommand(nID, lParam); 
	} 
} 
 
// If you add a minimize button to your dialog, you will need the code below 
//  to draw the icon.  For MFC applications using the document/view model, 
//  this is automatically done for you by the framework. 
 
void CShorDlg::OnPaint()  
{ 
	if (IsIconic()) 
	{ 
		CPaintDC dc(this); // device context for painting 
		 
		SendMessage(WM_ICONERASEBKGND, (WPARAM) dc.GetSafeHdc(), 0); 
		 
		// Center icon in client rectangle 
		int cxIcon = GetSystemMetrics(SM_CXICON); 
		int cyIcon = GetSystemMetrics(SM_CYICON); 
		CRect rect; 
		GetClientRect(&rect); 
		int x = (rect.Width() - cxIcon + 1) / 2; 
		int y = (rect.Height() - cyIcon + 1) / 2; 
		 
		// Draw the icon 
		dc.DrawIcon(x, y, m_hIcon); 
	} 
	else 
	{ 
		CDialog::OnPaint(); 
	} 
} 
 
// The system calls this to obtain the cursor to display while the user drags 
//  the minimized window. 
HCURSOR CShorDlg::OnQueryDragIcon() 
{ 
	return (HCURSOR) m_hIcon; 
} 
 
void CShorDlg::OnButton1()  
{    
	 
	CButton *p1; 
	p1=(CButton*)GetDlgItem(IDC_BUTTON1); 
	 
	srand(time(NULL)); 
	UpdateData(true); 
	int n=m_edit1; 
	if (n%2 == 0){  
		MessageBox("所输入的数必须是奇数!"); 
		m_edit1=0; 
		UpdateData(false); 
		 
    }  
	//Exit if the number is prime. 
	 else if (TestPrime(n)) { 
		MessageBox("所输入的数不能是质数!"); 
		m_edit1=0; 
		UpdateData(false); 
		 
	} 
	//Prime powers are prime numbers raised to integral powers. 
	//Exit if the number is a prime power. 
	else if (TestPrimePower(n)) { 
		MessageBox("所输入的数不能是质数的幂!"); 
		m_edit1=0; 
		UpdateData(false); 
		 
	} 
	else if(n<15){ 
		MessageBox("所输入的数必须大于或等于15!"); 
		m_edit1=0; 
		UpdateData(false); 
		 
	} 
	else 
	{ 
	p1->SetWindowText("执行中..."); 
	p1->EnableWindow(false); 
	int x = 0; 
	x = 1+ (int)((n-1)*(double)(rand())/(double)RAND_MAX); 
	while (GCD(n,x) != 1 || x == 1) { 
		x = 1 + (int)((n-1)*(double)(rand())/(double)RAND_MAX); 
	} 
	m_edit2.Format("随机数x的值为:%d",x); 
	UpdateData(false); 
	int q; 
	q = GetQ(n); 
	m_edit3.Format("q为满足大于等于N^2,小于2N^2的数,取值为:%d",q); 
	UpdateData(false); 
	QuReg * reg1 = new QuReg(RegSize(q) - 1); 
	m_edit4.Format("寄存器1的大小为:%d",RegSize(q)); 
	UpdateData(false); 
	int * modex = new int[q];      
	Complex *collapse = new Complex[q]; 
	Complex tmp; 
	Complex *mdx = new Complex[(int)pow(2,RegSize(n))]; 
	QuReg *reg2 = new QuReg(RegSize(n));  
	m_edit5.Format("建立寄存器2,大小为%d",RegSize(n)); 
    UpdateData(false); 
	int tmpval; 
	int value; 
	double c,m;  
	int den; 
	int p; 
	int e,a,b, factor; 
	int done = 0; 
	int tries = 0; 
	/*while (!done) { 
		if (tries >= 5) { 
			MessageBox( "There have been five failures, giving up." ); 
			exit(0); 
		} 
    }*/ 
	reg1->SetAverage(q - 1); 
		tmp.Set(1,0); 
		for (int i = 0 ; i < q ; i++) { 
			tmpval = modexp(x,i,n); 
			modex[i] = tmpval; 
			mdx[tmpval] = mdx[tmpval] + tmp; 
		} 
		reg2->SetState(mdx); 
		reg2->Norm(); 
		value = reg2->DecMeasure(); 
		for ( i = 0 ; i < q ; i++) { 
			if (modex[i] == value) { 
				collapse[i].Set(1,0); 
			} else { 
				collapse[i].Set(0,0); 
			} 
		} 
		reg1->SetState(collapse); 
		reg1->Norm();  
		m_edit6= "进行傅立叶变换!" ; 
		UpdateData(false); 
		DFT(reg1, q); 
	 
		m = reg1->DecMeasure(); 
		done = 1; 
		if (m == 7) { 
		    MessageBox("测量结束,试验失败!"); 
				m_edit1=0; 
	m_edit2=""; 
	m_edit3=""; 
	m_edit4=""; 
	m_edit5=""; 
	m_edit6=""; 
	m_edit7=""; 
	m_edit8=""; 
	m_edit9=""; 
	m_edit10=""; 
	CButton *p; 
	p=(CButton*)GetDlgItem(IDC_BUTTON1); 
	p->SetWindowText("执行"); 
	p->EnableWindow(true); 
	UpdateData(false); 
			done = 0; 
			UpdateData(false); 
		} 
		if (m == -1) { 
			MessageBox("没有测量到任何东西,这次试验失败!\n再试一次"); 
				m_edit1=0; 
	m_edit2=""; 
	m_edit3=""; 
	m_edit4=""; 
	m_edit5=""; 
	m_edit6=""; 
	m_edit7=""; 
	m_edit8=""; 
	m_edit9=""; 
	m_edit10=""; 
	CButton *p; 
	p=(CButton*)GetDlgItem(IDC_BUTTON1); 
	p->SetWindowText("执行"); 
	p->EnableWindow(true); 
	UpdateData(false); 
			done = 0; 
		} 
		if (done) { 
			c = (double)m  / (double)q; 
			den = denominator(c, q); 
			p = (int)floor(den * c + 0.5); 
            
			UpdateData(false); 
			cout << "measured " << m <<", approximation for " << c << " is " 
				<< p << " / " << den << endl << flush; 
			if (den % 2 == 1 && 2 * den < q ){ 
				 
				p = 2 * p;  
				den = 2 * den; 
			} 
			e = a = b = factor = 0; 
			if (den % 2 == 1)  {  
				MessageBox("得到的周期是奇数,这次试验失败!\n再试一次"); 
					m_edit1=0; 
	m_edit2=""; 
	m_edit3=""; 
	m_edit4=""; 
	m_edit5=""; 
	m_edit6=""; 
	m_edit7=""; 
	m_edit8=""; 
	m_edit9=""; 
	m_edit10=""; 
	CButton *p; 
	p=(CButton*)GetDlgItem(IDC_BUTTON1); 
	p->SetWindowText("执行"); 
	p->EnableWindow(true); 
	UpdateData(false); 
				cout <<  "Odd period found.  This trial failed."  
					<< "  Trying again." << endl << flush; 
				done = 0; 
			} else { 
				m_edit7.Format( "可能的周期是 %d", den); 
				e = modexp(x, den / 2, n);     
				a = (e + 1) % n;      
				b = (e + n - 1) % n;  
				m_edit8.Format("(%d^%d+1)mod %d = %d   (%d^%d-1)mod %d = %d",x,den/2,n,a,x,den/2,n,b); 
				cout << x << "^" << den / 2 << " + 1 mod " << n << " = " << a  
					<< "," << endl  
					<< x << "^" << den / 2 << " - 1 mod " << n << " = " << b  
					<< endl << flush; 
				factor = max1(GCD(n,a),GCD(n,b)); 
			} 
		} 
		if (factor == -1) { 
			MessageBox("错误, 零不能做分母!再试一次!"); 
				m_edit1=0; 
	m_edit2=""; 
	m_edit3=""; 
	m_edit4=""; 
	m_edit5=""; 
	m_edit6=""; 
	m_edit7=""; 
	m_edit8=""; 
	m_edit9=""; 
	m_edit10=""; 
	CButton *p; 
	p=(CButton*)GetDlgItem(IDC_BUTTON1); 
	p->SetWindowText("执行"); 
	p->EnableWindow(true); 
	UpdateData(false); 
			cout << "Error, tried to calculate n mod 0 for some n.  Trying again." 
				<< endl << flush; 
			done = 0; 
		} 
		 
		if ((factor == n || factor == 1) && done == 1) { 
			MessageBox("找到无用因子1,再试一次!"); 
				m_edit1=0; 
	m_edit2=""; 
	m_edit3=""; 
	m_edit4=""; 
	m_edit5=""; 
	m_edit6=""; 
	m_edit7=""; 
	m_edit8=""; 
	m_edit9=""; 
	m_edit10=""; 
	CButton *p; 
	p=(CButton*)GetDlgItem(IDC_BUTTON1); 
	p->SetWindowText("执行"); 
	p->EnableWindow(true); 
	UpdateData(false); 
			cout << "Found trivial factors 1 and " << n  
				<< ".  Trying again." << endl << flush; 
			done = 0; 
		} 
		if (factor != 0 && done == 1) { 
			p1->SetWindowText("结果"); 
			m_edit9.Format("%d=%d*%d",n,factor,n/factor); 
			UpdateData(false); 
			cout << n << " = " << factor << " * " << n / factor << endl << flush; 
		} else if (done == 1) { 
			MessageBox("错误,找到因子0,再试一次!"); 
				m_edit1=0; 
	m_edit2=""; 
	m_edit3=""; 
	m_edit4=""; 
	m_edit5=""; 
	m_edit6=""; 
	m_edit7=""; 
	m_edit8=""; 
	m_edit9=""; 
	m_edit10=""; 
	CButton *p; 
	p=(CButton*)GetDlgItem(IDC_BUTTON1); 
	p->SetWindowText("执行"); 
	p->EnableWindow(true); 
	UpdateData(false); 
			cout << "Found factor to be 0, error.  Trying again." << endl  
				<< flush; 
			done = 0; 
		} 
		tries++; 
		 
		} 
		 
   
   
} 
 
void CShorDlg::OnButton2()  
{ 
	m_edit1=0; 
	m_edit2=""; 
	m_edit3=""; 
	m_edit4=""; 
	m_edit5=""; 
	m_edit6=""; 
	m_edit7=""; 
	m_edit8=""; 
	m_edit9=""; 
	m_edit10=""; 
	CButton *p; 
	p=(CButton*)GetDlgItem(IDC_BUTTON1); 
	p->SetWindowText("执行"); 
	p->EnableWindow(true); 
	UpdateData(false); 
}