First working implementation. Can now factor 15*3.

Fixed all known normalization errors.
Debugging code (and two reg dumps) still active.

quantum_gates.c

@@ -25,16 +25,15 @@ int quda_quantum_hadamard_gate(int target, quantum_reg* qreg) {
 		// Copy amplitude to created state
 		qreg->states[qreg->num_states+i].amplitude = qreg->states[i].amplitude;
 
-		// TODO: Look at overhead of conditional rmul vs negation
 		if(qreg->states[i].state & mask) {
 			qreg->states[i].amplitude = quda_complex_neg(qreg->states[i].amplitude);
 		}
-		/*// DEBUG BLOCK
+		// DEBUG BLOCK
 		if(quda_complex_abs_square(qreg->states[i].amplitude) > 1.0f) {
 			printf("PROBLEM (hadamard): state[%d].amplitude = (%f,%f)\n",i,
 					qreg->states[i].amplitude.real,qreg->states[i].amplitude.imag);
 		}
-		*/// END DEBUG BLOCK
+		// END DEBUG BLOCK
 	}
 
 	printf("HADAMARD: increasing # states from %d ",qreg->num_states); // DEBUG
@@ -226,6 +225,7 @@ void quda_quantum_controlled_phase_gate(int control, int target, quantum_reg* qr
 }
 
 void quda_quantum_controlled_rotate_k_gate(int control, int target, quantum_reg* qreg, int k) {
+	if(control-target != k-1) printf("C-r-k Violation!\n"); // DEBUG
 	float temp = QUDA_PI / (1 << (k-1));
 	complex_t c = { .real = cos(temp), .imag = sin(temp) };
 	uint64_t mask = 1 << control;

quantum_reg.c

@@ -130,17 +130,20 @@ int quda_quantum_reg_measure_and_collapse(quantum_reg* qreg, uint64_t* retval) {
 	return -1;
 }
 
+// TODO: Write optimized version that avoids unnecessary duplicate operations
 int quda_quantum_range_measure_and_collapse(int start, int end, quantum_reg* qreg, uint64_t* retval) {
 	int i;
 	uint64_t res = 0;
-	for(i=0;i<end;i++) {
+	for(i=start;i<end;i++) {
 		res |= quda_quantum_bit_measure_and_collapse(i,qreg) << i;
 	}
 
 	if(retval) {
 		*retval = res;
 	}
 
+	// Renormalize already handled within each bit collapse
+
 	return 0;
 }
 /*	// old range_measure_and_collapse - appears to be invalid
@@ -387,9 +390,73 @@ int quda_amplitude_coalesce(complex_t* dest, complex_t* toadd) {
 		if(dest->real == 0.0f) {
 			dest->real = toadd->real;
 		} else {
-			if((dest->real > 0 && toadd->real < 0) || (dest->real < 0 && toadd->real > 0)) {
+			if(dest->real > 0 && toadd->real < 0) {
+				float temp = dest->real*dest->real - toadd->real*toadd->real;
+				if(temp < 0) {
+					dest->real = -sqrt(temp);
+				} else {
+					dest->real = sqrt(temp);
+				}
+				renorm = 1;
+			} else if(dest->real < 0 && toadd->real > 0) {
+				float temp = toadd->real*toadd->real - dest->real*dest->real;
+				if(temp < 0) {
+					dest->real = -sqrt(temp);
+				} else {
+					dest->real = sqrt(temp);
+				}
+				renorm = 1;
+			} else if(dest->real > 0 && toadd->real > 0) {
+				dest->real = sqrt(dest->real * dest->real + toadd->real * toadd->real);
+			} else {
+				dest->real = -sqrt(dest->real * dest->real + toadd->real * toadd->real);
+			}
+		}
+		toadd->real = 0.0f;
+	}
+
+	if(toadd->imag != 0.0f) {
+		if(dest->imag == 0.0f) {
+			dest->imag = toadd->imag;
+		} else {
+			if(dest->imag > 0 && toadd->imag < 0) {
+				float temp = dest->imag*dest->imag - toadd->imag*toadd->imag;
+				if(temp < 0) {
+					dest->imag = -sqrt(temp);
+				} else {
+					dest->imag = sqrt(temp);
+				}
 				renorm = 1;
+			} else if(dest->imag < 0 && toadd->imag > 0) {
+				float temp = toadd->imag*toadd->imag - dest->imag*dest->imag;
+				if(temp < 0) {
+					dest->imag = -sqrt(temp);
+				} else {
+					dest->imag = sqrt(temp);
+				}
+				renorm = 1;
+			} else if(dest->imag > 0 && toadd->imag > 0) {
+				dest->imag = sqrt(dest->imag * dest->imag + toadd->imag * toadd->imag);
+			} else {
+				dest->imag = -sqrt(dest->imag * dest->imag + toadd->imag * toadd->imag);
 			}
+		}
+		toadd->imag = 0.0f;
+	}
+
+	return renorm;
+}
+
+/* Old (wrong) implementation - did not account for cancellation
+int quda_amplitude_coalesce(complex_t* dest, complex_t* toadd) {
+	int renorm = 0;
+	if(toadd->real != 0.0f) {
+		if(dest->real == 0.0f) {
+			dest->real = toadd->real;
+		} else {
+			if((dest->real > 0 && toadd->real < 0) || (dest->real < 0 && toadd->real > 0)) {
+				renorm = 1;
+			} // FIXME: negatives DO NOT CANCEL
 			dest->real = sqrt(dest->real * dest->real + toadd->real * toadd->real);
 		}
 		toadd->real = 0.0f;
@@ -409,6 +476,7 @@ int quda_amplitude_coalesce(complex_t* dest, complex_t* toadd) {
 
 	return renorm;
 }
+*/
 
 float quda_rand_float() {
 	return rand()/(float)RAND_MAX;

quantum_stdlib.c

@@ -49,10 +49,13 @@ int quda_weak_check_amplitudes(quantum_reg* qreg) {
 
 void quda_quantum_reg_dump(quantum_reg* qreg, char* tag) {
 	int i;
+	uint64_t mask = (1 << qreg->qubits)-1;
+	uint64_t smask = ((1 << qreg->scratch)-1) << qreg->qubits;
 	printf("QREG_DUMP: %d states\n",qreg->num_states);
 	for(i=0;i<qreg->num_states;i++) {
 		if(tag) printf("%s: ",tag);
-		printf("qreg->states[%d].state = %lu\n",i,qreg->states[i].state);
+		printf("qreg->states[%d].state = %lu (bits,scratch)=(%lu,%lu)\n",i,qreg->states[i].state,
+				qreg->states[i].state & mask,(qreg->states[i].state & smask) >> qreg->qubits);
 		if(tag) printf("%s: ",tag);
 		printf("qreg->states[%d].amplitude = (%f,%f)\n",i,qreg->states[i].amplitude.real,
 				qreg->states[i].amplitude.imag);
@@ -79,12 +82,13 @@ int quda_quantum_hadamard_all(quantum_reg* qreg) {
 	return quda_quantum_hadamard_range(0,qreg->qubits,qreg);
 }
 
+// Original implementation
 void quda_quantum_fourier_transform(quantum_reg* qreg) {
 	int q = qreg->qubits-1;
 	int i,j;
 	for(i=q;i>=0;i--) {
 		for(j=q;j>i;j--) {
-			printf("Performing c-R_k(%d,%d)\n",j,i); // DEBUG
+			printf("Performing c-R_%d (PI/%lu) on (%d,%d)\n",j-i+1,(uint64_t)1 << (j-i),j,i); // DEBUG
 			//quda_quantum_controlled_phase_gate(j,i,qreg);
 			quda_quantum_controlled_rotate_k_gate(j,i,qreg,j-i+1);
 			// DEBUG BLOCK
@@ -103,6 +107,7 @@ void quda_quantum_fourier_transform(quantum_reg* qreg) {
 			printf("violation following hadamard gate\n");
 			exit(1);
 		}
+		//quda_quantum_reg_dump(qreg,"H:");
 		// END DEBUG BLOCK
 	}
 
@@ -114,6 +119,42 @@ void quda_quantum_fourier_transform(quantum_reg* qreg) {
 	}
 }
 
+/*// TESTING (endianness)
+void quda_quantum_fourier_transform(quantum_reg* qreg) {
+	int q = qreg->qubits-1;
+	int i,j;
+	for(i=0;i<=q;i++) {
+		for(j=0;j<i;j++) {
+			printf("Performing c-R_%d (PI/%lu) on (%d,%d)\n",i-j+1,(uint64_t)1 << (i-j),j,i); // DEBUG
+			quda_quantum_controlled_rotate_k_gate(j,i,qreg,i-j+1);
+			// DEBUG BLOCK
+			int err2 = quda_weak_check_amplitudes(qreg);
+			if(err2) {
+				printf("violation following c-phase gate\n");
+				exit(2);
+			}
+			// END DEBUG BLOCK
+		}
+		printf("Performing hadamard(bit %d)\n",i); // DEBUG
+		quda_quantum_hadamard_gate(i,qreg);
+		// DEBUG BLOCK
+		int err = quda_weak_check_amplitudes(qreg);
+		if(err) {
+			printf("violation following hadamard gate\n");
+			exit(1);
+		}
+		// END DEBUG BLOCK
+	}
+
+	// TODO: Consider using SWAP gate here instead
+	for(i=0;i<qreg->qubits/2;i++) {
+		quda_quantum_controlled_not_gate(i,q-i,qreg);
+		quda_quantum_controlled_not_gate(q-i,i,qreg);
+		quda_quantum_controlled_not_gate(i,q-i,qreg);
+	}
+}
+*/
+
 // Classical functions
 
 void quda_classical_exp_mod_n(int x, int n, quantum_reg* qreg) {

shor.c

@@ -1,12 +1,11 @@
 /* shor.c: implementation of Shor's quantum algorithm for factorization
  */
 
-#include "quantum_stdlib.h"
-#include "quantum_reg.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <time.h>
 #include <math.h>
+#include "quantum_stdlib.h"
 #include "shor.h"
 
 /* TODO: Change input parsing and/or accepted parameters
@@ -27,18 +26,23 @@ int main(int argc, char** argv) {
 	}
 
 	// Find a number x relatively prime to n
+	//srand(13); // TESTING - guarantees |4> as output reg state
 	srand(time(NULL));
 	int x;
 	do {
 		x = rand() % N;
 	} while(quda_gcd_div(N,x) > 1 || x < 2);
 
-	x = 8; // DEBUG
+	//x = 8; // DEBUG
+	//x = 7; // TESTING
+	x = 13; // TESTING2
 	printf("Random seed: %i\n", x);
 
 	int L = qubits_required(N);
-	int width = qubits_required(N*N);
+	//int width = qubits_required(N*N);
 	//int width = 2*L+2;
+	//int width = 11; // TESTING
+	int width = 4; // TESTING2
 
 	printf("N = %i, %i qubits required\n", N, width+L);
 
@@ -67,18 +71,29 @@ int main(int argc, char** argv) {
 	printf("After quda_quantum_hadamard_all()\n");
 	/* END DEBUG */
 
-	quda_quantum_add_scratch(3*L+2,&qr1); // Full scratch space not needed for classical algorithm
-	//quda_quantum_add_scratch(L,&qr1); // effectively creates 'output' subregister for exp_mod_n()
+	//quda_quantum_add_scratch(3*L+2,&qr1); // Full scratch may not be needed for classical exp_mod_n
+	quda_quantum_add_scratch(L,&qr1); // effectively creates 'output' subregister for exp_mod_n() (+TESTING)
 	quda_classical_exp_mod_n(x,N,&qr1);
 
+	/*// TESTING - Verifies exactly for x = 7 (width = 11)
+	int outputs = 7;
+	dump_mod_exp_results(&qr1,&outputs);
+	exit(0);
+	*/// END TESTING
+
+	/*// TESTING 2 - Verfies exactly for x = 13, width = 4
+	dump_mod_exp_results(&qr1,NULL); // TESTING 2 (x=13,width=4)
+	exit(0);
+	*/// END TESTING2
+
 	/* DEBUG */
 	err = quda_check_normalization(&qr1);
 	err |= quda_weak_check_amplitudes(&qr1);
 	if(err != 0) {
 		printf("ERROR: ");
 	}
 	printf("After quda_classical_exp_mod_n()\n");
-	printf("Ok going into quantum_collapse_scratch()\n");
+	printf("Going into quantum_collapse_scratch()\n");
 	/* END DEBUG */
 
 	/* libquantum measures all of its scratch bits here for some reason.
@@ -88,8 +103,8 @@ int main(int argc, char** argv) {
 	 * significant register bits to store scratch (with no differentiation from regular bits).
 	 * Comment the next line for the opposite effect (no coalescing).
 	 */
-	//quda_quantum_collapse_scratch(&qr1);
-	quda_quantum_clear_scratch(&qr1);
+	// TESTING - this function is now correct
+	quda_quantum_collapse_scratch(&qr1); // Explain implicit measurement; stupid not to call this
 
 	/* DEBUG */
 	err = quda_check_normalization(&qr1);
@@ -99,7 +114,9 @@ int main(int argc, char** argv) {
 	}
 	printf("After quantum_collapse_scratch()\n");
 	printf("Going into quantum_fourier_transform()\n");
-	//quda_quantum_reg_dump(&qr1,"BEFORE_FOURIER");
+	//qsort(qr1.states,qr1.num_states,sizeof(quantum_state_t),qstate_compare); // TESTING
+	// TESTING - dump verified as correct for (x,width) in {(7,11),(13,4)}
+	quda_quantum_reg_dump(&qr1,"BEFORE_FOURIER");
 	/* END DEBUG */
 
 	quda_quantum_fourier_transform(&qr1);
@@ -112,7 +129,7 @@ int main(int argc, char** argv) {
 	}
 	printf("After quda_quantum_fourier_transform()\n");
 	printf("Going into reg_measure_and_collapse\n");
-	//quda_quantum_reg_dump(&qr1,"AFTER_FOURIER");
+	quda_quantum_reg_dump(&qr1,"AFTER_FOURIER");  // TESTING - results appear incorrect
 	/* END DEBUG */
 
 	uint64_t result;
@@ -121,10 +138,13 @@ int main(int argc, char** argv) {
 		printf("Invalid result (normalization error).\n");
 		return -1;
 	} else if(result == 0) {
+		// NOTE: This is actually a valid result for 15 with (x=7,width=11) ~.25 prob
+		// Creates fraction 0/1, expands to 0/2, 2 is a valid period
 		printf("Measured zero.\n");
 		return 0;
 	}
 
+	// FIXME - Fractional expansion may be incorrect for our implementation
 	uint64_t denom = 1 << width;
 	quda_classical_continued_fraction_expansion(&result,&denom);
 
@@ -162,13 +182,34 @@ int main(int argc, char** argv) {
 	return 0;
 }
 
+// Testing functions
+void dump_mod_exp_results(quantum_reg* qreg, int* num_outputs) {
+	int outputs;
+	if(num_outputs) {
+		outputs = *num_outputs;
+	} else {
+		outputs = qreg->num_states;
+	}
+
+	uint64_t mask = (1 << qreg->qubits)-1;
+	uint64_t smask = ((1 << qreg->scratch)-1) << qreg->qubits;
+
+	uint64_t val,sval;
+	int i;
+	for(i=0;i<outputs;i++) {
+		val = qreg->states[i].state & mask;
+		sval = (qreg->states[i].state & smask) >> qreg->qubits;
+		printf("state[%d]: r_input = %lu, r_output =  %lu\n",i,val,sval);
+	}
+}
+
 // Classical functions
 
 int qubits_required(int num) {
-	int j = 1;
 	int i;
-	for(i=1;j < num;i++) {
-		j <<= 1;
+	num >>= 1;
+	for(i=1;num > 0;i++) {
+		num >>= 1;
 	}
 	return i;
 }

shor.h

@@ -1,6 +1,15 @@
 /* shor.c: header for implementation of Shor's quantum algorithm for factorization
  */
 
+#include "quantum_reg.h"
+
+// Testing functions
+
+/* Prints each qreg state and its corresponding modular exponentiation (stored in scratch).
+ * Only prints 'num_outputs' entries (or all if NULL).
+ */
+void dump_mod_exp_results(quantum_reg* qreg, int* num_outputs);
+
 // Classical functions (perform purely non-quantum computation)
 
 /* Determines number of qubits required to store the given number */